Implementation of a large elasto-plastic strain formalism in an industrial finite element code
Romain Carré
Royal Institute of Technology, 100 44 Stockholm, Sweden
Master thesis 2017
Abstract : This master thesis gives a brief sum up of the implementation of a large strain formulation in an existing industrial finite element code. The logarithmic strain tensor additively decomposed in an elastic and a plastic contribution was selected as a basis of this formulation. The main advantage of this formalism is that the existing plastic resolution algorithms written in the small strain formalism can be re-used. Only a pre-processing phase (construction of the logarithmic strain tensor) and a post-processing of the stress tensor is needed. This formulation, detailed in this paper was implemented in FORTRAN and integrated to the existing code. The final implementation was then validated thanks to some benchmarks with Abaqus.
1. Introduction
Analytical solutions do not exist for every structure;
some numerical solutions must be calculated instead.
The finite element method (FEM) is a powerful numerical tool to analyze mechanical structures. The linear FEM is usually enough to deal with the majority of the structures. Nevertheless, some advanced finite element codes must sometimes be used to solve the problems involving some mechanical nonlinearities.
Those can appear because of the material behavior (plasticity), can be created by some contacts or can be geometrical (significant changes in the shape of the structure). The aim of this paper is to implement a formalism dealing with the large strain nonlinearities and integrate it in an existing industrial finite element code. A simple example of application is the analysis of a metallic tensile specimen. For a high enough loading, some plastic areas will appear and the section of the specimen will shrink significantly. Those phenomena will generate some nonlinearities that cannot be treated with a linear code.
The existing finite element code takes into account the material nonlinearities so the general procedure to solve a problem involving some nonlinearities is already implemented in the code. For metallic material, a large strain analysis must include an elasto-plastic analysis. One of the objectives of this study is to take advantage as much as possible of what is already implemented and validated in the code. The general nonlinear resolution procedure (already implemented) is presented in Appendix A.
2. Selection of the formalism 1) Hypotheses of the formalism
In this paper, the main basic hypotheses involved in the finite element formulation will be used.
Displacements and strains can be large but only the nonlinear membrane effects will be considered and a decoupling between the nonlinear membrane terms and the linear bending terms is assumed. The hypothesis is valid if the finite elements are small enough and if the rotations remain small.
The code is supposed to be used with some metallic materials whose hardening behavior can be isotropic, kinematic or both. Some isotropic or anisotropic plasticity criteria can be used. In this paper, the formalism is presented for the case of an isotropic hardening behavior and a Von Mises plasticity criterion.
2) Choice of the formulation
a. Kinematic description
In order to take into account the effects of a large transformation involving some large displacements and some large strains, a distinction between the initial state or initial configuration and the deformed configuration must be made. Thus, the path of a material point belonging to the solid must be traced along the transformation.
The Lagrangian kinematic description was chosen because it involves the notion of trajectories well adapted to the deformable solids.
Among the sub-families of the Lagrangian kinematic description, the total Lagrangian description was selected. In this kinematic description, the initial configuration is used as the reference configuration.
This simple kinematic description brings a limitation.
Indeed, this description is more adapted to moderate strains (20-30%, it is the case for elasto-plastic materials) and would be inadequate for the “flow- like” analyses occurring in some manufacturing process. This kinematic description is also unreliable to model some topology changes. Nevertheless, this finite element code is not used to treat those kinds of applications.
According to this description, assuming a given transformation, a material point P is considered. The position of P in the reference configuration is called 𝑿. The position of P in the deformed configuration corresponding to the instant t can be followed with
𝒙 = 𝒙(𝑿 , 𝑡) (1)
The displacement field is defined as the difference between the deformed and initial positions:
𝒖 = 𝒙 − 𝑿 (2)
If the relative motion between the material points P and P’ located in the position 𝑿 and 𝑿 + 𝒅𝑿in the reference configuration is studied, one can write
𝑑𝒙 =𝜕𝒙
𝜕𝑿 ∙ 𝑑𝑿 = 𝑭 ∙ 𝑑𝑿 (3)
This leads to the definition of the transformation gradient tensor 𝑭 =𝜕𝑿𝜕𝒙 which can be rewritten as
𝑭 = [
𝜕𝑥
𝜕𝑋𝜕𝑦
𝜕𝑋𝜕𝑧
𝜕𝑋
𝜕𝑥
𝜕𝑌𝜕𝑦
𝜕𝑌𝜕𝑧
𝜕𝑌
𝜕𝑥
𝜕𝑍𝜕𝑦
𝜕𝑍𝜕𝑧
𝜕𝑍]
(4)
b. General formulation of the large strain formalism
The formulation described in this paper was selected for several reasons. The main advantage of this formulation is that it can allow the use of the algorithms in charge of the plastic resolution written in the small strain formalism. Indeed, the small strain formulation relies on some properties of the linearized strain tensor. This formalism is developed
so that the logarithmic strain tensor keeps those properties :
Multiplicative decomposition of the transformation gradient tensor (see the definition in Equation 3) in an elastic part 𝑭𝒆 and a plastic part 𝑭𝒑 so that 𝑭 = 𝑭𝒆∙ 𝑭𝒑
Additive decomposition of the strain tensor in an elastic part 𝜺𝒆𝒍 and a plastic part defined as 𝜺𝒑𝒍=1
2∙ ln (𝑭𝒑𝑻∙ 𝑭𝒑) so that the total deformation 𝜺 = 𝜺𝒆𝒍+ 𝜺𝒑𝒍= 12∙ ln (𝑭 𝑻∙ 𝑭)
Null trace for the plastic strain tensor
𝑭𝒑 describes the stress-free plastic configuration after the transformation. 𝑭𝒆 describes the extra elastic transformation experienced by the material to match the global configuration described by 𝑭. The additive decomposition of the strain tensor is always possible and it just depends on the right choice of 𝜺𝒆𝒍.
The last point is ensured by a mathematical property of the logarithmic function. From the expression of the logarithmic strain tensor given in Equation 8, it is possible to write
𝑡𝑟(𝜺) = ∑1 2∙
𝟑
𝒊=𝟏
ln(𝜆𝑖) (5)
Moreover, it comes from the mathematical properties of the determinant and the definition in (7):
det(𝑭) = √det (𝑪) = √𝜆1∙ 𝜆2∙ 𝜆3
ln(det(𝑭)) =12∙ (ln (𝜆1) + ln (𝜆2) + ln (𝜆3)) thus det(𝑭) = exp ( 𝑡𝑟(𝜺))
(6)
If the last expression is used with 𝑭𝒑 and 𝛆𝐩𝐥, because the plastic transformation is isochoric, its determinant is equal to 1 and thus the trace of the plastic strain tensor is equal to 0.
This implementation was based on the formalism presented by C. MIEHE in Reference 1. This formulation takes advantage of the additive decomposition of the logarithmic strain tensor in elasto-plasticity to develop an efficient resolution algorithm.
3. Mathematical formulation
The finite element formulation is presented in this paper for a 3D iso-parametric element.
Some mathematical operators used in the formulation are also presented in Appendix E.
The different steps of the resolution algorithm are presented in this section. A summarized version of the resolution procedure is available in Appendix A.
1) Strain tensor
The logarithmic strain tensor is derived from the transformation gradient tensor 𝑭 for a given displacement field. After the calculation of 𝑭, the symmetric right Green-Cauchy tensor 𝑪 is defined as
𝑪 = 𝑭𝑇∙ 𝑭 (7)
The theoretical formulation of the logarithmic strain tensor is
𝜺 =1
2∙ ln (𝑪) (8)
In order to actually implement the logarithmic strain tensor, an efficient way is to diagonalize C and to apply the logaritmic function to the eigenvalues.
So in a second step, an eigenvalue problem must be solved to obtain the eigenvalues λi and the eigenvectors Ni of C.
Then the logarithmic strain tensor is implemented as
𝜺 = ∑3 𝑒𝑖∙ 𝑵𝒊
𝑖=1 ⨂ 𝑵𝒊
Where 𝒆𝒊=12∙ ln (𝜆𝑖) for i=1:3
(9)
One can notice that the logarithmic strain tensor is a member of the Seth-Hill strain tensor family (m tends towards 0 in the formula given in Equation 10 with 𝑰 the identity matrix).
𝜺 = 1
𝑚∙ (𝑪𝑚2− 𝑰) (10)
This family was developed in order to respect the requirements to elaborate a strain tensor:
Symmetric second order tensor;
Null for a rigid body motion;
Can be expressed as monotonic and continuous function of the displacement gradients and in a unique way;
Must degenerate in the small strain formalism when the strains are small.
2) Stress tensor
In order to treat correctly the plastic resolution, the mathematical formulation of the energy involved in the problem must be performed carefully. To ensure this, the work conjugate tensors must be used. The elastic energy in the system will be the same
regardless of the choice of the strain tensor if its dual stress tensor is used.
Two equivalent conventions can be used to model the energy in the system. The difference between those depends on the choice of the volume of integration.
For a given strain tensor, the integration can be performed on the deformed volume (current convention) or on the initial volume (Lagrangian convention). For the same strain tensor, the stress tensor is then adapted to ensure that the energy involved is the same independently of the choice of the convention.
In the particular case of the logarithmic strain tensor, the elastic power can be written as
𝑃 = ∫ 𝑻′: 𝜺̇
𝑽 ∙ 𝑑𝑉 = ∫ 𝑑𝑒𝑡(𝑭) ∙ 𝑻′: 𝜺̇
𝑽𝒐 ∙ 𝑑𝑉
= ∫ 𝑻 : 𝜺̇
𝑽𝒐 ∙ 𝑑𝑉 (11)
𝑻′ is then the work conjugate stress tensor of the logarithmic strain tensor 𝜺 in the current convention and 𝑻 is the work conjugate tensor of 𝜺 in the Lagrangian convention.
From Equation 11, one can derive a simple relationship between 𝑻 to 𝑻′:
𝑻 = 𝑑𝑒𝑡(𝑭) ∙ 𝑻′ (12)
3) Material behavior law
The material behavior law must be well defined between the two conjugated tensors. In pure traction, 𝑻′ degenerates in the stress “ 𝐹𝑆 ”, the true stress (Cauchy stress).
If the case of an elasto-plasticity behavior law is considered, the material data must be used in the code with true stresses and true strains (log strain).
Because the material experimental test data are often given in terms of true stresses and true strains, those data can be directly sent as an entry in the finite element code. Thus, using 𝑻′ instead of 𝑻 makes more physical sense for the material behavior law and the plastic resolution.
The elastic material behavior matrix 𝑫 is then defined as the link between the 𝑻′stresses and the logarithmic strain so that
𝑻′ = 𝑫 ∙ 𝑬𝑷𝑺
(13) Thanks to Equation 13, the stress tensor corresponding to a given strain state can be calculated.
4) Plastic resolution
After this pre-processing phase, the plastic resolution algorithms written in the small strain formalism can be used.
Those algorithms give the plastic increment corresponding to the given displacement field. Then, it can be used to calculate the hardening and the plastic strain. Knowing the plastic strain and the total strain, it is easy to derive the elastic strain and thus the stresses thanks to the behavior matrix 𝑫.
Before the plastic resolution, the accumulated plastic strain updated at the previous equilibrium 𝜺𝒑𝒐 is known. 𝜺𝒑𝒐 must be removed from the total strain tensor constructed thanks to Equation 9 :
𝑬𝑷𝑺 = 𝜺 − 𝜺𝒑𝒐 (14)
The matrix 𝑫 is used to calculate the elastic prediction stresses 𝑻𝒆𝒑′ which will be used as an input of the plastic resolution algorithm. The elastic prediction stress tensor is defined as
𝑻′𝒆𝒑= 𝑫 ∙ 𝑬𝑷𝑺
(15) 𝜺𝒑𝒍, 𝜺𝒆𝒍 and 𝑻′ are the outputs of the plastic resolution.
5) Post processing of the stress tensor
The Piola-Kirchhoff 2 stress tensor, 𝕊, is widely used in the Lagrangian nonlinear mechanics. This stress tensor is defined as the work conjugate of the well- known Green-Lagrange strain tensor 𝜸 (specially used in large displacement and small strain analyses). This strain tensor is a member of the Seth-Hill family with m=2. Its expression is then given by the formula
𝜸 =1
2(𝑪 − 𝑰) (16)
𝑻′ is not classic in continuum mechanics and is more difficult to use. This tensor will be converted into 𝕊 during a post-processing phase after the plastic resolution. This is performed in order to calculate the internal forces and the tangent stiffness matrix in a classical way.
A relation between 𝑻 and 𝕊 is derived thanks to the Principle of the Virtual Powers (PVP).
𝑃 = ∫ 𝑻: 𝜺̇
𝑉𝑜 ∙ 𝑑𝑉 = ∫ 𝕊: 𝜸̇ ∙ 𝑑𝑉
𝑉𝑜
(17)
𝑃 = ∫ 𝑻:𝜕𝜺
𝜕𝑪: 𝑪̇ ∙ 𝑑𝑉
𝑉𝑜
= ∫ 𝕊:1 2∙ 𝑪̇ ∙ 𝑑𝑉
𝑉𝑜
(18)
𝑻:𝜕𝜺
𝜕𝑪= 𝕊:1 2 𝕊 = 𝑻: (2 ∙𝜕𝜺
𝜕𝑪) = 𝑻: 𝓟
(19)
One of the main difficulties of the formulation is the construction and the computation of the 4th-order operator 𝓟. The derivation of 𝓟, presented in Appendix B, leads to the following expression :
𝓟 = ∑1
𝜆𝑖∙ 𝑴𝒊⨂𝑴𝒊+ ∑ ∑ 𝜗𝑖𝑗 3
𝑗≠𝑖
∙ 𝑮𝒊𝒋 3
𝑖=1 3
𝑖=1
with 𝜗𝑖𝑗=𝑒𝜆1−𝑒2
𝑖−𝜆𝑗, 𝑴𝒊= 𝑵𝒊⨂ 𝑵𝒊 and 𝐺𝑖𝑗𝑎𝑏𝑑𝑐= 𝑀𝑖𝑎𝑐∙ 𝑀𝑗𝑏𝑑+ 𝑀𝑖𝑎𝑑∙ 𝑀𝑗𝑏𝑐
(20)
Since the operator 𝓟 is defined between 𝑻 and 𝕊 , the stresses 𝑻′ must be converted into 𝑻 thanks to Equation 12. To do so, the determinant of the transformation gradient 𝑭 must be calculated.
In 3D, deriving the determinant is easy since 𝑭 takes into account the transformation in the 3 directions. In 1D or 2D, the calculation of the determinant is more complex. Indeed the transformation is only considered in one or two directions.
The multiplicative decomposition of the transformation gradient is used to calculate the complete determinant. The use of one of the mathematical properties of the determinant gives
detF = detF𝑒∙ detF𝑝 (21) The plastic transformation is isochoric, thanks to the property of this formalism detailed in Section 2.2.b, detFp= 1. And so
detF = detF𝑒= exp (𝜀𝑥𝑥𝑒+ 𝜀𝑦𝑦𝑒+ 𝜀𝑧𝑧𝑒)
(22) In 2D under the plane stress assumptions, σzz= 0 means 𝜀𝑧𝑧𝑒= −1−𝜈𝜈 ∙ (𝜀𝑥𝑥𝑒+ 𝜀𝑦𝑦𝑒).
The missing strain 𝜀𝑧𝑧𝑒can then be used to calculate the complete determinant of the transformation. A similar procedure can be performed for the plane strain assumption.
6) Derivation of the internal forces
Once the conversion is achieved, one can calculate the internal forces in a classical way with the well-known Green-Lagrange and PK2 tensors. The Voigt notation is used to simplify the formulation.
The internal forces are defined as the derivative of the elastic energy with respect to the nodal displacement vector 𝑼.
The elastic energy is defined as
𝑈𝑒𝑙=1 2∙ ∫ 𝜸𝑻∙
𝑉𝑜 𝕊 ∙ 𝑑𝑉
(23) The interpolation matrix 𝑩 links the Green-Lagrange strain tensor 𝜸 to 𝑼 so that
𝑑𝜸 = 𝑩 ∙ 𝑑𝑼 with 𝑩 = 𝑩𝑳+ 𝑩𝑵𝑳
The expressions of 𝑩𝑳 and 𝑩𝑵𝑳 are available in Appendix C.
Once the interpolation matrix B is obtained, the internal forces can be calculated as
𝑭𝒊𝒏𝒕 =𝜕𝑈𝑒𝑙
𝜕𝑼
(24)
The mathematical formulation of the derivation of a scalar with respect to a vector is given in Appendix E.
𝑭𝒊𝒏𝒕 =𝝏𝑼𝒆𝒍
𝝏𝑼
=𝟏
𝟐∙ ∫ (𝝏𝕊
𝝏𝑼)
𝑻
∙
𝑽𝒐 𝜸 ∙ 𝒅𝑽 +𝟏
𝟐∙ ∫ (𝝏𝜸
𝝏𝑼)
𝑻
∙ 𝕊 ∙ 𝒅𝑽
𝑽𝒐
=𝟏
𝟐∙ ∫ (𝝏𝕊
𝝏𝜸∙𝝏𝜸
𝝏𝑼)
𝑻
∙
𝑽𝒐 𝜸 ∙ 𝒅𝑽 +𝟏
𝟐∙ ∫ 𝑩𝑻∙ 𝕊 ∙ 𝒅𝑽
𝑽𝒐
=𝟏
𝟐∙ ∫ (𝝏𝜸
𝝏𝑼)
𝑻
∙ (𝝏𝕊
𝝏𝜸)
𝑻
∙
𝑽𝒐 𝜸 ∙ 𝒅𝑽 +𝟏
𝟐
∙ ∫ 𝑩𝑻∙ 𝕊 ∙ 𝒅𝑽
𝑽𝒐
= ∫ 𝑩𝑻∙ 𝕊 ∙ 𝒅𝑽
𝑽𝒐
(25)
7) Tangent Stiffness matrix
The tangent stiffness matrix was also studied in order to reduce the number of iterations needed to converge towards the solution. Indeed, as it is presented in Appendix A, this matrix is used to derive the displacement increment that will be considered in the next iteration. The material nonlinear behavior is not fully taken into account in this matrix.
The tangent stiffness matrix is defined as the derivatives of the residuals (difference between the nodal internal and external forces) with respect to the nodal displacement vector 𝑼. Here again, the Green- Lagrange and PK2 tensors are used to calculate the tangent stiffness matrix in the classical way.
In the simple case where the external nodal forces are independent of 𝑼, it comes that
𝑲𝒕𝒈=𝜕𝑭𝒊𝒏𝒕
𝜕𝑼 = 𝜕
𝜕𝑼(∫ 𝑩𝑻∙ 𝕊 ∙ 𝑑𝑉
𝑉𝑜 )
(26)
The complete derivation of the expression of the stiffness matrix is too long to be written in this paper, but it is done in the same way of the derivation of the internal forces. Using the sixth order operator Ł defined as = 4 ∙𝜕𝜕𝑪2𝜺2 .
The final expression of the tangent stiffness matrix is 𝑲𝒕𝒈 = ∫ 𝑩𝑻
𝑉𝑜
∙ 𝓟𝑇: 𝑫: 𝓟 ∙ 𝑩 ∙ det(𝑭) ∙ 𝑑𝑉
+ ∫ 𝑩𝑻
𝑉𝑜
∙ 𝑻: Ł ∙ 𝑩 ∙ 𝑑𝑉
+ ∫ 𝑮𝑻
𝑉𝑜 ∙ 𝑴 ∙ 𝑮 ∙ 𝑑𝑉
(27)
The terms 𝑻: Ł, 𝑮 and 𝑴 are explicated in Appendix D.
4. Benchmark validation
The described formalism was first implemented in a Matlab prototype to validate the resolution process for simple elements. A formulation for a 2-node-bar element, a 4-node-membrane element and a 8-node- cubic element was implemented. The validation was made thanks to some benchmarks with Abaqus and Ansys with the corresponding elements. The simulations were launched on the same models with the same simulation parameters such as the convergence tolerance or the number of steps. The benchmark cases involved some pure traction, some combinations of shear and traction in elasticity and in elasto-plasticity.
Once the main aspects of the formalism were validated with the Matlab prototype, an implementation in FORTRAN directly in the finite element code was performed. The same simple benchmarks were performed again to ensure a correct FORTRAN implementation. Then, some more complex benchmarks validated definitely the integration of the implementation in the code.
Two “real-life” cases were studied.
1) 2D-case: Holed plate in traction
A plate with a hole in its center is pulled at its 2 ends.
More than 21 000 elements (triangles and quadrangles) are involved in the model. The displacement in the longitudinal direction of the plate is measured at one end on the location identified by a
blue cross. All the degrees of freedom (DOF’s) are blocked in the normal direction. At the extremities, the DOF’s in the X-direction are blocked. The DOF’s in the Y-direction are blocked for the 2 nodes in the middle of the plate identified by 2 black triangles. Two external forces are distributed on the black lines of nodes. Some kinematic constraints are applied to keep those lines straight.
The material used is Aluminium and the behavior law used comes from some physical testing. The simulations are launched on the same mesh with the same parameters with Abaqus and the large strain code. The results in displacements are presented on Figure 1 and 2.
Figure 1: Force-displacement plot from the benchmark with Abaqus on the holed plate
Less than 0.11% of difference exists between the 2 codes in displacement.
Figure 2: Displacement plot given by the larget strain formalism
2) 3D-case: Cylinder in compression
This 3D-simulation is performed on a cylinder (diameter 15 mm and a height of 25 mm). A 10 mm- vertical displacement is imposed on the nodes of the higher face of the cylinder whereas the nodes of the
lower face are blocked in the vertical direction. The DOF’s of the node in the center of the cylinder are blocked in the X and Y directions. More than 12 000 elements (cubic and prismatic elements) are involved in the simulation. The material used is Aluminium. The reaction forces on the blocked face are calculated and compared between Abaqus and the large strain code.
The results in reaction forces are presented in Figure 3.
Figure 3: Displacement-force plot from the benchmark with Abaqus on cylinder
The difference is smaller than 0.3% in reaction forces.
3) Tangent stiffness matrix
It was difficult to validate the implementation of the tangent stiffness matrix thanks to some benchmarks with Abaqus. Indeed, Abaqus works like a black box and knowing the technics actually used by Abaqus is difficult (B.F.G.S, linesearch). The material nonlinearity may also be taken into account in Abaqus.
The convergence was tested on the elementary benchmarks used to validate the calculation of the internal forces. Without any material nonlinearity, the convergence is significantly faster when the stiffness matrix is updated every 10 iterations than when it is only factorized at the beginning of the analysis. The number of iterations is approximately divided by two for simple traction case on one element not experiencing large movements.
If the plasticity is considered, the effect of the re- factorization is less significant since the nonlinearities from the material are more important than the geometrical nonlinearities and the material ones are not treated in the stiffness matrix. The results are presented in Figure 4.
Figure 4: Convergence study on a 3D elementary model in traction
For all the particular cases investigated, the first term of Equation 27 seems more significant. These elementary cases in pure traction and shear mixed with traction include only some large strain nonlinearities, the gain in iterations would be more important if some large displacements were also included (the tangent matrix would re-orientate the stiffness in the right directions improving the convergence). To ensure a valid implementation of the tangent stiffness matrix, some benchmarks (buckling of shells) are in progress.
5. Conclusion
The large strain formalism was, at first, implemented in a Maltab prototype and then in FORTRAN in the existing finite element code. An additive formulation of the logarithmic strain tensor was chosen to be able to use the algorithms for the plastic resolution developed in small strain. Some benchmarks were made to compare the large strain finite element code to Abaqus to ensure the implementation of the formalism is valid. The development of a well- constructed tangent stiffness matrix decreases the number of iterations needed to converge towards the solution. Another way to speed up the convergence would be to adapt the linesearch technics for the formalism.
With this implementation, a distinction between the initial configuration and the deformed configuration is made and the variation of the section is taken into account. This makes the simulations more realistic and closer to the physical testing. One drawback is the complexity of the formalism and the heavier computation process.
6. References
[1] C. Miehe, N. Apel, M. Lambrecht, Anisotropic additive plasticity in the logarithmic strain space : modular kinematic formulation and implementation based on incremental minimization principles for standard materials, Institut für Mechanik (Bauwesen) Lehrstuhl I, Universität Stuttgart, 2002.
[2] C. Miehe, M. Lambrecht, Algorithms for computation of stresses and elasticity moduli in terms of Seth–Hill’s family of generalized strain tensors, Institut für Mechanik (Bauwesen) Lehrstuhl I, Universität Stuttgart, 2001.
[3] R. Bargellini, Documentation Code Aster: R.03.24, Modèles de grandes déformations GDEF_LOG et GDEF_HYPO_ELAS, 2013.
[4] M. Latorre, F.J. Montáns, On the interpretation of the logarithmic strain tensor in an arbitrary system of representation, International Journal of Solids and Structures 51, 2014.
[5] R. de Borst, M. A. Crisfield, J.J.C. Remmers, C.V.
Verhoosel, Non-linear finite element analysis of solids and structures, Second edition, 2012.
[6] T. Helfer, C. Ling, ÉCRITURE DE LOIS DE COMPORTEMENT MÉCANIQUE EN GRANDES TRANSFORMATIONS AVEC LE GÉNÉRATEUR DE CODE MFRONT,2014.
[7] K.J. Bathe, Finite element procedures, Prentice Hall, 1996.
Appendix A: Nonlinear resolution procedure
General nonlinear resolution procedure in the small strain formalism
General nonlinear resolution procedure in the large strain formalism
Fext imposed ith iteration, U(i)
Internal force (Fint) Plastic resolution
Construction of the Strain tensor
Residuals=Fext-Fint <Tol ?
Factorization of the tangent stiffness matrix (if needed) Elastic prediction stress tensor
Fext imposed ith iteration, U(i)
U(i+1)=U(i)+V
Construction of 𝑭 then 𝑪
Elastic prediction stress tensor
Plastic resolution
Conversion of 𝑻′ into 𝑻
Conversion of 𝑻 into PK2 (with 𝓟)
Internal force (Fint)
Residuals=Fext-Fint <Tol ?
𝑽 = 𝑲𝒕𝒈−1∙ 𝑹𝒆𝒔𝒊𝒅𝒖𝒂𝒍𝒔
U(i+1)=U(i)+V NO
YES
NO
Factorization of the tangent stiffness matrix (if needed)
YES
Construction of the Strain tensor
𝑽 = 𝑲𝒕𝒈−1∙ 𝑹𝒆𝒔𝒊𝒅𝒖𝒂𝒍𝒔
Appendix B : Derivation of the operator P
The fourth order operator 𝓟 is defined as 𝓟 = 2 ∙𝜕𝑪𝜕𝜺 . In the case where C has 3 different eigenvalues:
𝓟 = 2 ∙𝜕𝜺
𝜕𝑪= 2 ∙ 𝜕
𝜕𝑪(∑ 𝑒𝑖∙ 𝑴𝒊
3
𝑖=1
) = 2 ∙ ∑ (𝜕𝑒𝑖
𝜕𝑪⨂ 𝑴𝒊+ 𝑒𝑖⨂ 𝜕𝑴𝒊
𝜕𝑪 )
3
𝑖=1
= 2 ∙ ∑ (𝜕𝑒𝑖
𝜕𝜆𝑖∙𝜕𝜆𝑖
𝜕𝑪⨂ 𝑴𝒊+ 𝑒𝑖⨂ 𝜕𝑴𝒊
𝜕𝑪 )
3
𝑖=1
With the classical matrix derivation mathematical formulas, it comes for i varying from 1 to 3,
𝜕𝜆𝑖
𝜕𝑪= 𝑵𝒊⨂ 𝑵𝒊= 𝑴𝒊 , 2 ∙𝜕𝑵𝜕𝑪𝒊= ∑ 𝜆1
𝑖−𝜆𝑗
3𝑗≠𝑖 ∙ 𝑵𝒋⨂ (𝑵𝑖⨂ 𝑵𝒋+ 𝑵𝒋⨂ 𝑵𝒊)
and 2 ∙𝜕𝑴𝜕𝑪𝒊= ∑ 𝜆1
𝑖−𝜆𝑗
3𝑗≠𝑖 ∙ (𝑮𝒊𝒋+ 𝑮𝒋𝒊)
with 𝐺𝑖𝑗𝑎𝑏𝑑𝑐= 𝑀𝑖𝑎𝑐∙ 𝑀𝑗𝑏𝑑+ 𝑀𝑖𝑎𝑑∙ 𝑀𝑗𝑏𝑐
with a, b , c and d some integers varying from 1 to 6.
Using those formulas in the expression of 𝓟 gives the final expression presented in Equation 20. In the case of equal eigenvalues, the limits of the expressions are considered.
Appendix C : Interpolation matrix B
𝑩 is defined so that 𝑑𝜸 = 𝑩 ∙ 𝑑𝑼 with 𝑩 = 𝑩𝑳+ 𝑩𝑵𝑳. The expressions of the (6*24) matrices 𝑩𝑳 and 𝑩𝑵𝑳 are :
𝑩𝑳=
[
𝜕𝑁1
𝜕𝑋 0 0
0 𝜕𝑁1
𝜕𝑌 0
0 0 𝜕𝑁1
𝜕𝑍
⋯ ⋯
𝜕𝑁8
𝜕𝑋 0 0 0 𝜕𝑁8
𝜕𝑌 0
0 0 𝜕𝑁8
𝜕𝑍 0 𝜕𝑁1
𝜕𝑍
𝜕𝑁1
𝜕𝑁1 𝜕𝑌
𝜕𝑍 0 𝜕𝑁1
𝜕𝑋
𝜕𝑁1
𝜕𝑌
𝜕𝑁1
𝜕𝑋 0
⋯ ⋯
0 𝜕𝑁8
𝜕𝑍
𝜕𝑁8
𝜕𝑁8 𝜕𝑌
𝜕𝑍 0 𝜕𝑁8
𝜕𝑥𝑋
𝜕𝑁8
𝜕𝑌
𝜕𝑁8
𝜕𝑋 0 ]
𝑩𝑵𝑳=
[
𝜕𝑢
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑋 𝜕𝑣
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑋 𝜕𝑤
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑢 𝜕𝑋
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑌 𝜕𝑣
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑌 𝜕𝑤
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑢 𝜕𝑌
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑍 𝜕𝑣
𝜕𝑍∙𝜕𝑁1
𝜕𝑍 𝜕𝑤
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑍
⋯ ⋯
𝜕𝑢
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑍 +𝜕𝑢
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑌
𝜕𝑣
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑍 +𝜕𝑣
𝜕𝑍∙𝜕𝑁1
𝜕𝑌
𝜕𝑤
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑍 +𝜕𝑤
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑢 𝜕𝑌
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑋 +𝜕𝑢
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑍
𝜕𝑣
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑋 + 𝜕𝑣
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑍
𝜕𝑤
𝜕𝑍 ∙𝜕𝑁1
𝜕𝑋 +𝜕𝑤
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑢 𝜕𝑍
𝜕𝑌 ∙𝜕𝑁1
𝜕𝑋 +𝜕𝑢
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑌
𝜕𝑣
𝜕𝑌∙𝜕𝑁1
𝜕𝑋 + 𝜕𝑣
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑌
𝜕𝑤
𝜕𝑌∙𝜕𝑁1
𝜕𝑋 +𝜕𝑤
𝜕𝑋 ∙𝜕𝑁1
𝜕𝑌
⋯ ⋯ ]
Appendix D : Tangent stiffness matrix
Some terms appearing in the expressions of the tangent stiffness matrix are explicated here:
The matrix 𝑮 is obtained thanks to the The last term involved in the expression of the derivatives of the shape functions : tangent matrix is the matrix 𝑴 defined a
With 𝑰𝟑 the 3*3 identity matrix and 𝕊𝒊𝒊 the different components of the PK2 stress tensor.
The term 𝑻: Ł is a fourth order tensor. The numerical implementation takes advantage of the Voigt notations to reduce this term to a second order tensor implemented as a matrix. Its expression is given in the case where the eigenvalues of C are different :
Ł= 4 ∙𝜕2𝜺
𝜕𝑪2 𝑻:Ł= ∑ 𝑓𝑖∙ (𝑴𝒊: 𝑻) ∙ 𝑴𝒊⨂𝑴𝒊
𝟑
𝒊
+ ∑ ∑ ∑ 2 ∙ 𝜂 ∙ (𝓕𝒊(𝒋𝒌)+ 𝓕(𝒋𝒌)𝒊)
3
𝑘≠𝑖𝑘≠𝑗 3
𝑗≠𝑖 3
𝑖
+ ∑ ∑ 2 ∙ 𝜉𝑖𝑗∙ (𝓕𝒊(𝒋𝒋)+ 𝓕(𝒋𝒋)𝒊+ 𝓕𝒋(𝒊𝒋)+ 𝓕(𝒊𝒋)𝒋+ 𝓕𝒋(𝒋𝒊)+ 𝓕(𝒋𝒊)𝒋)
3
𝑗≠𝑖 3
𝑖
With 𝓕𝒂(𝒃𝒄)𝒊𝒋𝒌𝒍= 𝑴𝒂𝒊𝑘∙ (𝑴𝒃∙ 𝑻 ∙ 𝑴𝒄)𝑗𝑙+ 𝑴𝒂𝒊𝑙∙ (𝑴𝒃∙ 𝑻 ∙ 𝑴𝒄)𝒋𝒌
With 𝑓𝑖= −𝜆2
𝑖2 , 𝜈𝑖𝑗=𝑒𝜆𝑖−𝑒𝑗
𝑖−𝜆𝑗 , 𝜉𝑖𝑗=𝜈𝑖𝑗−0.5∙
1 𝜆𝑖 𝜆𝑖−𝜆𝑗 et
𝜂 = ∑ ∑ ∑ 𝑒𝑖
2 ∙ (𝜆𝑖− 𝜆𝑗) ∙ (𝜆𝑖− 𝜆𝑘)
3 𝑘≠𝑖𝑘≠𝑗 3
𝑗≠𝑖 3
𝑖
If the eigenvalues are equal, the limit value is considered.
Appendix E: Definition of the mathematical operators used
Different mathematical operators are used to implement the large strain formalism. Some of them are briefly presented in this section.
Tensor product
The notation for this mathematical operator is « ⨂ ». The tensor resulting from this operation with a tensor 𝑨 (order p) and a tensor 𝑩 (order q) is a p+q order tensor. In particular:
If 𝑨 and 𝑩 are two first order tensors (vectors), the resulting tensor 𝑪 = 𝑨 ∗ 𝑩𝑇 will be a second order tensor.
If 𝑨 and 𝑩 are two second order tensors (matrices), the resulting tensor 𝑪 = 𝑨 ⨂ 𝑩 will be a fourth order tensor defined as
𝑪(𝑖, 𝑗, 𝑘, 𝑙) = 𝑨(𝑖, 𝑘) ∗ 𝑩(𝑗, 𝑙)
with the integers I, j , k and l varying in the rows and columns of the matrices 𝑨 and 𝑩.
Doubly contracted product
The notation for this mathematical operator is « : ». The tensor resulting from this operation with a tensor 𝑨 (order p≥2) and a tensor 𝑩 (order q≥2) is a p+q-4 order tensor. In particular:
If 𝑨 and 𝑩 are two second order tensors (matrices), the resulting tensor 𝑪 = 𝑨: 𝑩 will be a zero order tensor (scalar) defined as 𝐶 = ∑ ∑ 𝐴𝑖 𝑗 𝑖𝑗∙ 𝐵𝑗𝑖
If 𝑨 and 𝑩 are two tensors (order 4 et 2), the resulting tensor 𝑪 = 𝑨: 𝑩 will be a second order tensor (matrix) defined as 𝐶𝑖𝑗 = ∑ ∑ 𝐴𝑘 𝑙 𝑖𝑗𝑘𝑙∙ 𝐵𝑘𝑙
If 𝑨 and 𝑩 are two fourth order tensors, the resulting tensor 𝑪 = 𝑨: 𝑩 will be a fourth order tensor defined as 𝐶𝑖𝑗𝑚𝑛 = ∑ ∑ 𝐴𝑘 𝑙 𝑖𝑗𝑘𝑙∙ 𝐵𝑘𝑙𝑚𝑛
Derivation of a scalar with respect to a vector
Let u, v and x be two column vectors. The derivation of the scalar 𝒖𝑇∙ 𝒗 is given by the following formula:
𝜕
𝜕𝒙(𝒖𝑇∙ 𝒗 ) = (𝜕𝒗
𝜕𝒙)
𝑇
∙ 𝒖 + (𝜕𝒖
𝜕𝒙)
𝑇
∙ 𝒗
