Second gradient models as a particular case of microstructured models: a large strain finite elements analysis (original) (raw)

Second gradient models as a particular case of microstructured models: a large strain finite elements analysis

Takashi MATSUSHIMA a { }^{\text {a }}, René CHAMBON b { }^{\text {b }}, Denis CAILLERIE b { }^{\text {b }}
a{ }^{a} University of Tsukuba, Tsukuba, Japan formerly visiting researcher at Laboratoire 3S
b { }^{\text {b }} Laboratoire 3S Grenoble, Université Joseph-Fourier, Institut national polytechnique, CNRS UMR 5521, BP 53, 38041 Grenoble cedex 9, France

(Reçu le 4 mai 1999, accepté après révision le 27 mai 1999)

Abstract

Second gradient models are usually non local. Recently local second gradient models have been developed. They can be seen as a particular case of microstructured media which obey a kinematic constraint. It is then easy to build up a large strain finite element method useful for 1D as well as for 2D or 3D problems. A one-dimensional example is finally given. © 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS finite elements / second gradient / large strains / localization / microstructure

Matériaux avec microstructures, le cas du second gradient: une méthode d’éléments finis en grande déformation

Résumé. Les modèles de second gradient les plus étudiés sont non locaux. Récemment nous avons développé des modèles de second gradient locaux. Il est intéressant de les considérer comme des milieux avec microstructure cinématiquement contraints. Ce point de vue permet de construire naturellement une méthode d’éléments finis en grande déformation adaptée. On présente ainsi une application unidimensionnelle de cette méthode d’éléments finis facilement généralisable à deux dimensions. © 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
eléments finis / second gradient / grandes déformations / localisation / microstructure

Version française abrégée

1. Introduction : matériaux avec microstructures

En reprenant les travaux de Germain [1,2], on peut définir la cinématique d’un milieu avec microstructure avec le champ de déplacement usuel noté uiu_{i} et un champ de tenseur microstructural noté fijf_{i j}. Si l’on note par * des quantités virtuelles et si l’on définit Fij⋆=∂ui⋆∂xjF_{i j}^{\star}=\frac{\partial u_{i}^{\star}}{\partial x_{j}} et Dij⋆=12(Fij⋆+Fji⋆)D_{i j}^{\star}=\frac{1}{2}\left(F_{i j}^{\star}+F_{j i}^{\star}\right), le principe des puissances virtuelles s’écrit comme dans l’équation (1), où σij\sigma_{i j} est appelée contrainte macroscopique, τij\tau_{i j} qui n’est pas nécessairement symétrique est associée à la microstructure, et χijk\chi_{i j k} est appelée double contrainte. We⋆W_{e}^{\star}, le travail virtuel des forces extérieures n’est pas détaillé ici (voir [1,2]).

Note présentée par Évariste SANCHEZ-PALENCIA.

2. Matériaux avec microstructure contraints : milieux du second gradient, application unidimensionnelle

Supposons l’égalité du tenseur microstructural fijf_{i j} et du gradient des déplacements usuel (cf. équation (2)), dans ce cas le principe des travaux virtuels s’écrit comme dans l’équation (3). On peut l’utiliser directement pour développer une méthode d’éléments finis (voir [3,4]). Cette méthode utilise des éléments d’Hermite, ce qui ne peut être fait que pour des problèmes unidimensionnels. Une autre façon de procéder est d’écrire la contrainte cinématique (2) sous forme faible et d’utiliser des multiplicateurs de Lagrange. Le principe des puissances virtuelles est alors écrit dans les équations (4) et (5) où λij\lambda_{i j} est le champ de multiplicateur. Dans ce cadre, on peut construire des modèles de comportement locaux dont l’énergie associée se définit de la même manière que les travaux virtuels. C’est une différence importante avec la plupart des modèles de second gradient qui, excepté les modèles que nous avons développés dans [3,4] et les modèles de Fleck et Hutchinson [5], sont non locaux [6].

On se restreint désormais à des problèmes unidimensionnels comme dans [3] et [4]. Les équations (4) et (5) se particularisent alors en les équations (6) et (7) valides pour le domaine correspondant à la position courante du milieu.

3. Analyse numérique

On discrétise le temps et à la fin de chaque pas de temps les équations de l’équilibre ((6) et (7)) doivent être satisfaites. On fait une hypothèse sur l’histoire cinématique entre le début du pas et la fin du pas de temps, puis on entame une procédure itérative. Soit CnC_{n} la nième configuration finale supposée. L’équation (8) nous fournit un résidu Rn∗R_{n}^{*}. En linéarisant au voisinage de CnC_{n} on obtient l’équation (10) dont les inconnues sont les différences entre deux itérations successives notées par un Δ.ΔNn\Delta . \Delta N_{n} et ΔMn\Delta M_{n} sont calculés par linéarisation consistante au sens de Simo et Taylor [7] exprimée en fonction de dΔue dxo\frac{\mathrm{d} \Delta u_{e}}{\mathrm{~d} x_{o}} et dΔee dxn\frac{\mathrm{d} \Delta e_{e}}{\mathrm{~d} x_{n}}. On obtient alors le système linéaire auxiliaire utilisé dans le processus itératif. Notons [D][D] la matrice utilisée dans l’équation (12).

4. Eléments finis

On utilise un élément parent pour calculer la contribution élémentaire de la version discrétisée de l’équation (12). Cet élément est défini pour ξ∈[−1,+1]\xi \in[-1,+1]. Les fonctions de formes pour uu sont : ϕ1(ξ)=\phi_{1}(\xi)= −ξ2+ξ22,ϕ2(ξ)=1−ξ2-\frac{\xi}{2}+\frac{\xi^{2}}{2}, \phi_{2}(\xi)=1-\xi^{2} et ϕ3(ξ)=ξ2+ξ22\phi_{3}(\xi)=\frac{\xi}{2}+\frac{\xi^{2}}{2}; elles sont aussi utilisées pour faire la correspondance géométrique entre l’élément parent et l’élément réel. Pour vv les fonctions de formes sont ψ1(ξ)=12−ξ2\psi_{1}(\xi)=\frac{1}{2}-\frac{\xi}{2} et ψ2(ξ)=12+ξ2.λ\psi_{2}(\xi)=\frac{1}{2}+\frac{\xi}{2} . \lambda est supposé constant. On obtient alors l’équation (13) valide aussi bien pour les valeurs virtuelles que pour les différences réelles. On appelle comme d’habitude [B][B] la matrice présente dans cette équation. Finalement l’expression de la matrice de rigidité élémentaire est écrite dans l’équation (14) et la partie du second membre non dépendant des conditions aux limites est détaillée dans l’équation (15)).

5. Application numérique

Soit ε\varepsilon la déformation courante, en grande déformation ε=du dx0\varepsilon=\frac{\mathrm{d} u}{\mathrm{~d} x_{0}} ce qui est différent de u′=du dxu^{\prime}=\frac{\mathrm{d} u}{\mathrm{~d} x}. En ce qui concerne le comportement, on suppose qu’en charge (u˙′⩾0)\left(\dot{u}^{\prime} \geqslant 0\right) on a N˙=Au˙′\dot{N}=A \dot{u}^{\prime} avec si ε<elim⁡,A=A1\varepsilon<e_{\lim }, A=A_{1} et si ε>elim⁡,A=A2\varepsilon>e_{\lim }, A=A_{2} (qui est négatif afin d’obtenir une possibilité de localisation). Pour la décharge (u˙′⩽0)\left(\dot{u}^{\prime} \leqslant 0\right), on suppose N˙=A1u˙′\dot{N}=A_{1} \dot{u}^{\prime}. La partie second gradient est modélisée par : M˙=Bu˙′′=B(1−v)v˙′+v˙v′(1−v)2\dot{M}=B \dot{u}^{\prime \prime}=B \frac{(1-v) \dot{v}^{\prime}+\dot{v} v^{\prime}}{(1-v)^{2}} car comme vv est identifié à u′,u˙′′u^{\prime}, \dot{u}^{\prime \prime} n’est pas égal à v˙′\dot{v}^{\prime}. La figure 1 traduit le comportement du modèle dans le cas des petites déformations. Dans ce dernier cas on connaît des solutions analytiques de problèmes aux limites [3, 4]. Une barre de 1 m de long est discrétisée en 100 éléments. Les conditions aux limites sont ua=U;ub=0u_{a}=U ; u_{b}=0 et M=0M=0 pour les deux extrémités aa et bb. On a pris A1=150;A2=−75;B=0,8;elim⁡=0,01A_{1}=150 ; A_{2}=-75 ; B=0,8 ; e_{\lim }=0,01 comme

dans les articles référencés. Analytiquement, on trouve 4 solutions post-pic différentes (au sens des courbes globales). La figure 2 montre la convergence de la valeur de la déformation à une extrémité de la barre pour des éléments de Hermite [3,4] et pour les éléments à trois champs présentés ici, dans le cas des petites déformations. La solution est la solution la plus raide (solution (4) dans [4]) et ce cas est le plus difficile à modéliser. On obtient de très bons résultats même dans ce cas difficile. La figure 3 est une comparaison entre les solutions (3) [4] en grandes déformations en utilisant les éléments d’Hermite ou les éléments à trois champs.

6. Conclusion

Les modèles locaux de second gradient représentent une solution pour régulariser les effets de radoucissement. Dans cette note, on a montré sur un exemple unidimensionnel qu’une méthode d’éléments finis à trois champs est efficace pour les modéliser. Une généralisation à deux dimensions a déjà été réalisée et donne d’aussi bons résultats.

1. Introduction: microstructured materials

Following Germain [1,2], kinematics of microstructured materials is described by a displacement field like in classical continuum denoted uiu_{i} plus a tensor field corresponding to the kinematics of the microstructure, denoted fijf_{i j}. Using the virtual work principle is the most convenient way to study such media. Indicating with a ⋆{ }^{\star} virtual quantities, defining Fij=∂ui⋆∂xjF_{i j}=\frac{\partial u_{i}^{\star}}{\partial x_{j}} and Dij⋆=12(Fij⋆+Fji⋆)D_{i j}^{\star}=\frac{1}{2}\left(F_{i j}^{\star}+F_{j i}^{\star}\right), the usual virtual strain, the virtual work principle and objectivity considerations yield [2]:

∫Ω(σijDij⋆+τij(fij⋆−Fij⋆)+χijk∂fij⋆∂xk)dv=We⋆\int_{\Omega}\left(\sigma_{i j} D_{i j}^{\star}+\tau_{i j}\left(f_{i j}^{\star}-F_{i j}^{\star}\right)+\chi_{i j k} \frac{\partial f_{i j}^{\star}}{\partial x_{k}}\right) \mathrm{d} v=W_{e}^{\star}

σij\sigma_{i j} is called here the macro stress, τij\tau_{i j} is an additional stress associated with the microstructure which is not necessarily symmetric and χijk\chi_{i j k} is called the double stress. We⋆W_{e}^{\star}, the external virtual work is not detailed here, it can be seen in [1,2][1,2].

2. Microstructured continuum with kinematics constraint: second gradient model

Let us assume that:

fij=Fij=∂ui∂xjf_{i j}=F_{i j}=\frac{\partial u_{i}}{\partial x_{j}}

In this case, assuming the same kinematic constraint for the corresponding virtual quantities, we get the virtual work principle for a second gradient model:

∫Ω(σijDij⋆+χijk∂2ui⋆∂xj∂xk)dv=We⋆\int_{\Omega}\left(\sigma_{i j} D_{i j}^{\star}+\chi_{i j k} \frac{\partial^{2} u_{i}^{\star}}{\partial x_{j} \partial x_{k}}\right) \mathrm{d} v=W_{e}^{\star}

A sraightforward way of building a corresponding finite element method is detailed in [3,4]. The drawback of such a method is well known; it implies the use of Hermitian finite elements which is easily tractable only for one dimensional problems. However, using the mathematical constraint (2) in a weak sense by mean of Lagrange multipliers is another way to consider such models. In this case the virtual work principle yields:

∫Ω[σij∂ui⋆∂xj+χijk∂fij⋆∂xk+λij(∂ui⋆∂xj−fij⋆)]dv=We⋆\int_{\Omega}\left[\sigma_{i j} \frac{\partial u_{i}^{\star}}{\partial x_{j}}+\chi_{i j k} \frac{\partial f_{i j}^{\star}}{\partial x_{k}}+\lambda_{i j}\left(\frac{\partial u_{i}^{\star}}{\partial x_{j}}-f_{i j}^{\star}\right)\right] \mathrm{d} v=W_{e}^{\star}

T. Matsushima et al.

and

∫Ωλij∗(∂ui∂xj−fij)dv=0\int_{\Omega} \lambda_{i j}^{*}\left(\frac{\partial u_{i}}{\partial x_{j}}-f_{i j}\right) \mathrm{d} v=0

where λij\lambda_{i j} is the (nonsymmetric) field of Lagrange multipliers. Clearly equation (1) gives the physical meaning of the Lagrange multipliers.

In all cases, models obtained by starting from the virtual work principle can be local models without any difficulties since the true work can be defined in the same manner as the virtual work. This is a salient difference with the most of used second gradient models (except the model studied in [3,4] and the FleckHutchinson model [5]) which are nonlocal ones [6].

3. A one-dimensional application

The previous theory is now restricted to a one-dimensional medium. xx is the only coordinate of the current configuration and uu is the corresponding displacement component. Micro deformation is then defined by vv. Ordinary stress is denoted by NN. The double stress is denoted by MM. For the sake of simplicity, the external forces are only the boundary conditions acting on the two ends of the medium a(t)a(t) and b(t)b(t). Equations (4) and (5) yield for a domain Ct=[a(t),b(t)]C_{t}=[a(t), b(t)] corresponding to the current position (at a given time tt ) of the one dimensional medium:

−∫Ct(Nt∂ut∗∂xt+Mt∂vt∗∂xt)dxt−∫Ctλt(∂ut∗∂xt−vt∗)dxt+(Nb−∂M∂xt∣b)ub∗−(Na−∂M∂xt∣a)ua∗+Mbvb∗−Mava∗=0∫Ctλt∗(∂ut∂xt−vt)dxt=0\begin{aligned} & -\int_{C_{t}}\left(N_{t} \frac{\partial u_{t}^{*}}{\partial x_{t}}+M_{t} \frac{\partial v_{t}^{*}}{\partial x_{t}}\right) \mathrm{d} x_{t}-\int_{C_{t}} \lambda_{t}\left(\frac{\partial u_{t}^{*}}{\partial x_{t}}-v_{t}^{*}\right) \mathrm{d} x_{t} \\ & +\left.\left(N_{b}-\frac{\partial M}{\partial x_{t}}\right|_{b}\right) u_{b}^{*}-\left.\left(N_{a}-\frac{\partial M}{\partial x_{t}}\right|_{a}\right) u_{a}^{*}+M_{b} v_{b}^{*}-M_{a} v_{a}^{*}=0 \\ & \int_{C_{t}} \lambda_{t}^{*}\left(\frac{\partial u_{t}}{\partial x_{t}}-v_{t}\right) \mathrm{d} x_{t}=0 \end{aligned}

These equations have to hold for any virtual kinematically admissible field. Solving a problem needs boundary conditions (i.e., the knowledge of uu or NN and vv or MM at the two ends of the domain for any time tt ).

4. Numerical analysis

A step by step finite element method is used to solve this problem. The balance equations (through the virtual work method) are written at the end of a step corresponding to time tt : see equations (6) and (7). For a given step an assumption is made for the kinematics between the beginning and the end of the step. It gives an estimation of the final configuration, and then iterations are performed to get equilibrium. Let CnC_{n} be the nnnth assumed final configuration corresponding to a given time step. Equation (8) gives us a residual Rn⋆R_{n}^{\star}.

−∫Cn(Nn∂u∗∂xn+Mn∂v∗∂xn−λn(∂u∗∂xn−v∗))dxn+∫Cnλ∗(∂un∂xn−vn)dxn+(Nb−∂M∂xt∣b)ub∗−(Na−∂M∂xt∣a)ua∗+Mbvb∗−Mava∗=Rn⋆\begin{aligned} & -\int_{C_{n}}\left(N_{n} \frac{\partial u^{*}}{\partial x_{n}}+M_{n} \frac{\partial v^{*}}{\partial x_{n}}-\lambda_{n}\left(\frac{\partial u^{*}}{\partial x_{n}}-v^{*}\right)\right) \mathrm{d} x_{n}+\int_{C_{n}} \lambda^{*}\left(\frac{\partial u_{n}}{\partial x_{n}}-v_{n}\right) \mathrm{d} x_{n} \\ & +\left.\left(N_{b}-\frac{\partial M}{\partial x_{t}}\right|_{b}\right) u_{b}^{*}-\left.\left(N_{a}-\frac{\partial M}{\partial x_{t}}\right|_{a}\right) u_{a}^{*}+M_{b} v_{b}^{*}-M_{a} v_{a}^{*}=R_{n}^{\star} \end{aligned}

A configuration Cn+1C_{n+1} is searched such that equation (9) is satisfied.

−∫Cn+1(Nn+1∂u∗∂xn+1+Mn+1∂v∗∂xn+1−λn+1(∂u∗∂xn+1−v∗))dxn+1+∫Cn+1λ∗(∂un+1∂xn+1−vn+1)dxn+1\begin{aligned} & -\int_{C_{n+1}}\left(N_{n+1} \frac{\partial u^{*}}{\partial x_{n+1}}+M_{n+1} \frac{\partial v^{*}}{\partial x_{n+1}}-\lambda_{n+1}\left(\frac{\partial u^{*}}{\partial x_{n+1}}-v^{*}\right)\right) \mathrm{d} x_{n+1} \\ & +\int_{C_{n+1}} \lambda^{*}\left(\frac{\partial u_{n+1}}{\partial x_{n+1}}-v_{n+1}\right) \mathrm{d} x_{n+1} \end{aligned}

ParseError: KaTeX parse error: Expected '\right', got 'EOF' at end of input: …}^{*}=0\right)

Substracting equation (9) from equation (8) and assuming that configuration CnC_{n} is close to configuration Cn+1C_{n+1} yields to the linearized equation (10), in which unknowns are the differences between iteration nn and iteration n+1n+1 denoted by a Δ\Delta.

∫Cn[ du⋆dxn dv⋆dxnv⋆][ΔNn−ΔλnΔMnΔλn+λn dΔun dxn]dxn+∫Cnλ⋆((1−vn)dΔun dxn−Δvn)dxn=−Rn⋆\int_{C_{n}}\left[\frac{\mathrm{~d} u^{\star}}{\mathrm{d} x_{n}} \frac{\mathrm{~d} v^{\star}}{\mathrm{d} x_{n}} v^{\star}\right]\left[\begin{array}{c} \Delta N_{n}-\Delta \lambda_{n} \\ \Delta M_{n} \\ \Delta \lambda_{n}+\lambda_{n} \frac{\mathrm{~d} \Delta u_{n}}{\mathrm{~d} x_{n}} \end{array}\right] \mathrm{d} x_{n}+\int_{C_{n}} \lambda^{\star}\left(\left(1-v_{n}\right) \frac{\mathrm{d} \Delta u_{n}}{\mathrm{~d} x_{n}}-\Delta v_{n}\right) \mathrm{d} x_{n}=-R_{n}^{\star}

ΔNn\Delta N_{n} and ΔMn\Delta M_{n} are computed by a consistent linearization in the sense of Simo and Taylor [7] of the stress point algorithm which uses the second gradient constitutive equation. This linearization is a function of dΔun dxn\frac{\mathrm{d} \Delta u_{n}}{\mathrm{~d} x_{n}} and dΔvn dxn\frac{\mathrm{d} \Delta v_{n}}{\mathrm{~d} x_{n}}.

[ΔNnΔMn]=[C11C12C21C22][dΔun dxn dΔvn dxn]\left[\begin{array}{l} \Delta N_{n} \\ \Delta M_{n} \end{array}\right]=\left[\begin{array}{ll} C_{11} & C_{12} \\ C_{21} & C_{22} \end{array}\right]\left[\begin{array}{l} \frac{\mathrm{d} \Delta u_{n}}{\mathrm{~d} x_{n}} \\ \frac{\mathrm{~d} \Delta v_{n}}{\mathrm{~d} x_{n}} \end{array}\right]

Finally we end up with the auxiliary linear system used for the iterative procedure.

∫Cn[ du⋆dxn dv⋆dxnv⋆λ⋆][C11C120−1C21C2200λn001(1−vn)0−10][dΔun dxn dΔvn dxnΔvnΔλn]dxn=−Rn⋆\int_{C_{n}}\left[\frac{\mathrm{~d} u^{\star}}{\mathrm{d} x_{n}} \frac{\mathrm{~d} v^{\star}}{\mathrm{d} x_{n}} v^{\star} \lambda^{\star}\right]\left[\begin{array}{cccc} C_{11} & C_{12} & 0 & -1 \\ C_{21} & C_{22} & 0 & 0 \\ \lambda_{n} & 0 & 0 & 1 \\ \left(1-v_{n}\right) & 0 & -1 & 0 \end{array}\right]\left[\begin{array}{c} \frac{\mathrm{d} \Delta u_{n}}{\mathrm{~d} x_{n}} \\ \frac{\mathrm{~d} \Delta v_{n}}{\mathrm{~d} x_{n}} \\ \Delta v_{n} \\ \Delta \lambda_{n} \end{array}\right] \mathrm{d} x_{n}=-R_{n}^{\star}

Let us denote [D][D] the square matrix used in the previous equation.

5. Finite elements

As usual, a parent element is used to compute equations (12) in a discretized manner. This element is defined by ξ∈[−1,+1]\xi \in[-1,+1]. In the parent element, uu is assumed to be a quadratic function, vv is assumed to be linear and λ\lambda is assumed to be constant. Nodal points for uu are the two ends of the element called node 1 and node 3 and the midpoint called node 2 ; the corresponding shape functions are ϕ1(ξ)=−ξ2+ξ22\phi_{1}(\xi)=-\frac{\xi}{2}+\frac{\xi^{2}}{2}, ϕ2(ξ)=1−ξ2\phi_{2}(\xi)=1-\xi^{2} and ϕ3(ξ)=ξ2+ξ22\phi_{3}(\xi)=\frac{\xi}{2}+\frac{\xi^{2}}{2}, which are used also for the mapping between the parent element and the current element itself. Nodal points for vv are the two ends of the element and the corresponding shape functions are ψ1(ξ)=12−12\psi_{1}(\xi)=\frac{1}{2}-\frac{1}{2} and ψ2(ξ)=12+ξ2\psi_{2}(\xi)=\frac{1}{2}+\frac{\xi}{2}. Consequently, we can get:

[dun dxn dvn dxnvnλn]=[dϕ1 dξ0 dϕ2 dξ0 dϕ3 dξ0 dϕn dξ dϕ2 dξ dϕ3 dξ0 dψ1 dξ000 dψ2 dξ dψn dξ dψ2 dξ0ψ1000ψ2000100][u∣ξ=−1v∣ξ=−1u∣ξ=0λ∣ξ=0u∣ξ=1v∣ξ=1]\left[\begin{array}{c} \frac{\mathrm{d} u_{n}}{\mathrm{~d} x_{n}} \\ \frac{\mathrm{~d} v_{n}}{\mathrm{~d} x_{n}} \\ v_{n} \\ \lambda_{n} \end{array}\right]=\left[\begin{array}{cccccc} \frac{\mathrm{d} \phi_{1}}{\mathrm{~d} \xi} & 0 & \frac{\mathrm{~d} \phi_{2}}{\mathrm{~d} \xi} & 0 & \frac{\mathrm{~d} \phi_{3}}{\mathrm{~d} \xi} & 0 \\ \frac{\mathrm{~d} \phi_{n}}{\mathrm{~d} \xi} & & \frac{\mathrm{~d} \phi_{2}}{\mathrm{~d} \xi} & & \frac{\mathrm{~d} \phi_{3}}{\mathrm{~d} \xi} & \\ 0 & \frac{\mathrm{~d} \psi_{1}}{\mathrm{~d} \xi} & 0 & 0 & 0 & \frac{\mathrm{~d} \psi_{2}}{\mathrm{~d} \xi} \\ \frac{\mathrm{~d} \psi_{n}}{\mathrm{~d} \xi} & & & & & \frac{\mathrm{~d} \psi_{2}}{\mathrm{~d} \xi} \\ 0 & \psi_{1} & 0 & 0 & 0 & \psi_{2} \\ 0 & 0 & 0 & 1 & 0 & 0 \end{array}\right]\left[\begin{array}{c} u \mid \xi=-1 \\ v \mid \xi=-1 \\ u \mid \xi=0 \\ \lambda \mid \xi=0 \\ u \mid \xi=1 \\ v \mid \xi=1 \end{array}\right]

T. Matsushima et al.

This equation is valid for the virtual quantities as well as for differences. Let us denote [B][B] the matrix used in the previous equation. It is straightforward to simplify [B][B] when small strains are assumed. Finally, the element stiffness matrix used in the auxiliary linear system for every iteration reads:

∫−11[B]T[D][B]dxn dξ dξ\int_{-1}^{1}[B]^{\mathrm{T}}[D][B] \frac{\mathrm{d} x_{n}}{\mathrm{~d} \xi} \mathrm{~d} \xi

and the part independent of the boundary conditions of the right hand side of the linear system of equation reads:

−∫−11[B]T[Nn−λnMnλn dun dxn−vn]dxn dξ dξ-\int_{-1}^{1}[B]^{\mathrm{T}}\left[\begin{array}{c} N_{n}-\lambda_{n} \\ M_{n} \\ \lambda_{n} \\ \frac{\mathrm{~d} u_{n}}{\mathrm{~d} x_{n}}-v_{n} \end{array}\right] \frac{\mathrm{d} x_{n}}{\mathrm{~d} \xi} \mathrm{~d} \xi

6. Numerical application

Let us define ε\varepsilon as the current strain; in large strain ε=du dx0\varepsilon=\frac{\mathrm{d} u}{\mathrm{~d} x_{0}} which is different from u′=du dxu^{\prime}=\frac{\mathrm{d} u}{\mathrm{~d} x}. As far as the constitutive equation is concerned, we assume that

N˙=Au˙′\dot{N}=A \dot{u}^{\prime}

for loading (if u˙′⩾0\dot{u}^{\prime} \geqslant 0 ). If ε<elim⁡,A=A1\varepsilon<e_{\lim }, A=A_{1}; and if ε>elim⁡,A=A2\varepsilon>e_{\lim }, A=A_{2} (which is negative, this model is then a softening one which is necessary to localize in 1D case). For unloading (if u˙′⩽0\dot{u}^{\prime} \leqslant 0 ) we assume:

N˙=A1u˙′\dot{N}=A_{1} \dot{u}^{\prime}

The second gradient part of the model is simply written in equation (18).

M˙=Bu˙′′=B(1−v)v˙′+v˙v′(1−v)2\dot{M}=B \dot{u}^{\prime \prime}=B \frac{(1-v) \dot{v}^{\prime}+\dot{v} v^{\prime}}{(1-v)^{2}}

Since we identify vv and u′u^{\prime} for any time t,u˙′′t, \dot{u}^{\prime \prime} is not equal to v˙′\dot{v}^{\prime}. The behavior of this model is sketched in figure 1 in the case of small strain assumption. Boundary value problems can be solved analytically with

Figure 1. Model behaviour.
Figure 1. Loi de comportement.

Figure 2. Accuracy of the solution as a function of the number of elements (small strain computation).
Figure 2. Précision de la solution fonction du nombre d’éléments (petites déformations).

Figure 3. Comparison between solutions obtained with hermitian elements and three-field elements (large strain computations).
Figure 3. Comparaison entre solutions par éléments de Hermite et par éléments à trois champs (grandes déformations).
such a model with the small strain assumption, and numerical solutions can be obtained using the present three-field elements as well as the previous Hermite finite elements (see [3,4]), either within the small strain assumption framework or without this restriction. The algorithm used is classical and can be seen for instance in [8]. A one-meter bar has been discretized in 100 elements. The boundary conditions are ua=Uu_{a}=U, ub=0u_{b}=0 and M=0M=0 for the both ends aa and bb. We assumed: A1=150,A2=−75,B=0.8,elim⁡=0.01A_{1}=150, A_{2}=-75, B=0.8, e_{\lim }=0.01 like in the referenced papers in which it is proved analytically that only four different (which means four different relationships between global force and global strain of the bar) post peak solutions are available. Figure 2 shows the rate of convergence of the strain value at one end of the bar for the Hermitian elements (presented in [3,4][3,4] ) and for the present three field elements. The (small strain) computed solution is the stiffer one (namely solution 4, see [4]) which is the most difficult to obtain. Obviously the new elements are less accurate, but even in this case the method is valid and accurate with a sufficiently fine discretization. Figure 3 shows a comparison between solution 3 (see [4]) obtained by using the three-field elements presented here and by Hermitian conforming elements for large strain computations. It is clear that the method advocated above is accurate enough to retrieve the good solutions.

7. Conclusion

Local second gradient model can be a solution to regularize softening models. However the simplest way to deal with such models from the numerical point of view can be used only for 1D problems. This paper presents a different way to compute large (or small) strain problems with a local second gradient model which is directly applicable for two or three-dimensional problems. This way has been proved to be efficient for the simple case of 1D problem. Generalization for 2D problem has already been performed and gave results as good as for the one-dimensional case.

References

[1] Germain P., La méthode des puissances virtuelle en mécanique des milieux continus: Première partie: théorie du second gradient, Journal de Mécanique 12 (1973) 235-274.
[2] Germain P., The method of virtual power in continuum mechanics, Part 2, Microstructure, SIAM J. Appl. Math. 25 (1973) 556-575.
[3] Chambon R. et al., Etude de la localisation unidimensionelle à l’aide d’un modèle de second gradient, C.R.A.S. IIb 323 (1996) 231−238231-238.

T. Matsushima et al.

[4] Chambon R. et al., One-dimensional localisation studied with a second grade model, Eur. J. Mech. A/Solids 17 (1998) 637-656.
[5] Fleck N.A., Hutchinson J.W., Strain gradient plasticity, Advances in Applied Mechanics 33 (1997) 295-361.
[6] de Borst R., Muhlhaus H., Computational strategies for gradient continuum models with a view to localisation of deformation, in: Proc. 4th Int. Conf. on Non-linear Eng. Comp., Bicanic et al. (Eds.), Pineridge Press, Swansea, 1992, pp. 239-260.
[7] Simo J.C., Taylor R.L., Consistent tangent operators for rate-independent elastoplasticity, Comput. Meth. Appl. Mech. Engrg. 33 (1985) 101-118.
[8] Chambon R., Crochepeyre S., Loss of uniqueness of the solution in a boundary value problem, in: N.U.M.O.G. V Balkema (Ed.), 1995, pp. 165-171.