Finite element formulation of plane stress plasticity

In [1], we stated the basic FEM procedures to deal with the elastic-plastic problem applicable to the general three-dimensional state of stress, plane strain, and an axisymmetric body. These problems are those in which all three components of the normal stress are free variables. However, in plane stress, the component perpendicular to the midplane of the body is equal to zero, but the component of the stress deviator in this direction is not zero, and these relationships cannot be used for this case without modification.

7.1 Plane stress linear elasticity

The plane stress assumption was introduced in the analysis of bodies in which one dimension was significantly smaller than the other dimensions. Consider a loaded body in equilibrium (Fig. 1), where the plane (x, y) is the plane of symmetry (midplane $S$) of the body. Let all external loads (including reactions) be distributed symmetrically to the relatively small thickness $t$ (x, y) of the body. Thus, the load resultants acting in the midplane, that is, the continuous line load $\mathbf{p}(x,y)$ and body load $\mathbf{b}(x,y)$, are independent of $z$. Under these conditions, it is enough to consider only the displacement functions of the points on the midplane of the body

The non-zero components of the stress tensor placed in the column vector $\boldsymbol{\sigma}$ are only the components that lie in the plane (x, y)

Hence, we say that the body (wall, membrane) is in a plane stress state.

In the column vector $\mathbf{\varepsilon}$ of the strain components, we will also consider only three components (In fact material reacts to the plane stress state by the transverse deformation $\varepsilon_z$, but this can be determined from the deformations in the midplane.)

After simplifications

the strains in the plane (x,y) are given by

From the physical (material) equations we have

and from the condition $\sigma_z = 0$

Inverting (7.6) we get the stress components

Now the strain $\varepsilon_z$ can be expressed by substitution (7.8) to (7.7)

The physical (constitutive) relationships in the matrix form are given by

In the equations above $E$ is the Young modulus, $\nu$ denotes the Poisson ratio and

is the shear modulus of the material.

7.2 Von Mises plane stress isotropic small strain plasticity with the linear hardening

In this case, we again consider only those strain components that lie in the midplane because $\varepsilon_z$ can be expressed by using these components, and for the plastic strain $\varepsilon_z^p$, we consider an incompressible material at plastic yielding ($\nu =0.5$). Thus, we have

where the splitting of $\varepsilon_z$ into elastic and plastic parts is introduced based on the assumption of small strain.

We need the von Mises equivalent stress $\bar{\sigma}$ for the plasticity condition

which at the plane stress state has form

where $\bar{\sigma}_k(\bar{\varepsilon}^p)$ is the uniaxial hardened yield stress, which is a function of the accumulated plastic strain $\bar{\varepsilon}^p$.

Then the vector of plastic yielding is given by

7.3 Stress increment specification

For the elastic test (trial) stress in the numerical incremental integration step from $n$ to $n+1$ we have

and the equations of stress components after the stress point projection onto the updated yield surface are

where $\Delta\lambda$ is the unknown incremental plastic multiplier.

After an arrangement of equations (7.17), we get

where

The unknown plastic multiplier $\Delta\bar{\lambda}$ has to satisfy the plasticity condition (7.13)

where for a linear hardening material the hardened yield stress $\bar{\sigma}_k$ is

depending on the initial yield stress $\sigma_k$, the accumulated plastic strain $\bar{\varepsilon}_p$, and the hardening modulus

where $E_t$ is the elastoplastic tangent modulus.

The relationship for the $\Delta\lambda$ determination has a more simple form when one uses the square of the condition (7.22)

Squaring the von Mises stress $\bar{\sigma}$ (7.14) we get

and after substituting (7.18)-(7.20)

Then the auxiliary function $q$ (7.25) is given by

The nonlinear equation (7.28), in fact, the other form of the consistency condition (7.22), enables us to solve $\Delta\bar{\lambda}$ using the Newton-Raphson iteration. Decomposition of the function $q$ in the reduced Taylor series gives the iterative relationship

where $q'$ is the derivative of function $q$ by $\Delta\bar{\lambda}$ and which form according to (7.28) and (7.23) is as follows

From (7.27) after extraction and differentiation by $\Delta\bar{\lambda}$ we get the last needed relationship for the iterative procedure

The calculation of $\Delta\bar{\lambda}$ and the resulting stress $\mathsf{\sigma}_{n+1}$ are shown in the following example (the numerical results were calculated using the Mathematica 7 program, as shown in Fig. 7.2).

Example 7.1

Consider a body made of steel material and loaded with plane stress. Its properties were approximated using the von Mises loading function and linear hardening model. Consider a material particle, which previous state is without stress, deformed by strains $\varepsilon_x = 0,002$, $\varepsilon_y = -\,0,001$ and $\gamma_{xy} = 2\varepsilon_{xy} = 0,002$ and determine from the explicit relationships given above the final values of the stress components in the next integration step, that is, $\mathsf{\sigma}_{n+1}$. Given is E = 200000 MPa, $\sigma_k$ = 200 MPa, the linear elastoplastic tangent modulus $E_t = 100000$ MPa, $\nu = 0,3$. The required accuracy was $f_{n+1} = \bar{\sigma}_{n+1} - \bar{\sigma}_k < 0,001$ MPa.

Solution

The test stress is

Test of plasticity condition (7.14) gives

Because the test elastic stress $\bar{\sigma}_0$ is greater than the initial yield stress $\sigma_k$ of the material point in the investigation, it follows that an elastic-plastic deformation will occur. The resulting stress in such case must be determined from backward integration relationships.

The initial values for the calculation of stress at the elastic-plastic loading were $\Delta\bar{\lambda}_0 = 0$ and $\bar{\sigma}_{k0} = \sigma_k = 200$. Then, according to (7.24) - (7.28), the iterative calculation was performed by the program (the used relations do not contain matrices, and so we have chosen program Mathematica 7), as shown in Fig. 7.2.

Fig. 7.2 Mathematica 7 program and results of the example 7.1 solution (Výsledky = Results)

The components of the resulting stress at the investigated point of the body are as follows:

The resulting (hardened) yield stress of the material and resulting the equivalent plastic strain are as follows:

Example 7.2

The previous example was solved using ANSYS Mechanical APDL.

Solution

The material particle of the body can be replaced by a single plane stress finite element of unit thickness with given material properties. The problem of the $\gamma_{xy}$ input can be solved by placing the element edges in the directions of the main strains $\varepsilon_1$ and $\varepsilon_2$ for which we have

and for the direction of $\varepsilon_1$ we receive

The example can be solved in the ANSYS interactive environment using the following steps:

1. Job name
Utility Menu>File>Change Jobname..., /FILNAM = PlaneStress1, OK;

2. Element type
Main Menu>Preprocessor>Element Type>Add/Edit/Delete, Add..., Solid Quad 8 node 183, OK, Close;

3. Material properties
Linear
Preprocessor>Material Props>Material Models, Structural, Linear, Elastic, Isotropic, EX = 2E5, PRXY = 0.3, OK,
Nonlinear
Nonlinear, Inelastic, Rate Independent, Isotropic Hardening Plasticity, Mises Plasticity, Bilinear, Yield Strss = 200, Tang Mod = 1E5, OK, Material, Exit;

4. Element keypoints (numbering is automatic)
Utility Menu>Work Plane>Offset WP by Increments..., Degrees = 16.845, 0, 0, OK;
Utility Menu>Work Plane>Local Coordinate Systems>Create Local CS>At WP Origin..., OK;
Preprocessor>Modeling>Create>Keypoints>In Active CS: X = 0, Y = 0, Apply, X = 1, Y = 0, Apply, X = 1, Y = 1, Apply, X = 0, Y = 1, OK;

5. Element area
Preprocessor>Modeling>Create>Areas>Arbitrary>Trough KPs: KP1,KP2,KP3,KP4, OK;

6. Element creating
Preprocessor>Meshing>Mesh Tool: Size Controls: Lines, Set: Pick All, NDIV = 1, OK; (Preprocessor>Meshing>Mesh Tool): Mesh, Pick All, Close;

7. $\varepsilon_1$ and $\varepsilon_2$ input
Preprocessor>Modeling>Move/Modify>Rotate Node CS>To Active CS, Pick All;
Main Menu>Solution>Define Loads>Apply>Structural>Displacement>On Nodes
List of Items: 1, 8, 6, OK, UX, Value = 0, OK, On Nodes 1, 3, 2, OK, UY, Value = 0, OK, On Nodes 2, 5, 4, OK, UX, Value = 0.0023027, OK, On Nodes 6, 7, 4, OK, UY, Value = -0.0013027, OK;

8. Solution
Main Menu>Solution>Analysis Type>Sol'n Controls...,Time at end of loadstep = 1, Number of substeps = 1, OK; Solution>Solve>Current LS, Solve Current Load Step, OK;

9. Stress results
Main Menu>General Postproc>List Results> Nodal Solution..., Stress, X-Component of stress, OK;

10. End of job
Ansys Toolbar>Quit>Save Geom+Loads, OK;

The results were identical to those of the analytical solution.

7.4 Matrix formulation

The analytical determination of the stress state in the elastic-plastic load step is indeed an elegant and simple procedure; however, if we are interested in creating an FEM program for solving elastic-plastic plane stress problems, it is necessary to establish basic relations and save the necessary values for the next load step in the matrix form. In such case, the state variables are written to the column matrices (vectors)

where at the deformation variables we consider tensor components for shear strain ($\gamma_{xy}=2\varepsilon_{xy}$) and for the deviatoric stress component we will implement

where

Now we can express the von Mises equivalent stress as follows

The plasticity condition is given by

and its preferred squared form is

The easiest way to get the increment of the plastic deformation vector is by using the formula for deriving the scalar quadratic form

where for the vector of plasticity in this case we have

The increment of plastic deformation is

In an increment step from time $t_n$ to time $t_{n+1}$, the values of the elastic predictor are

where $\Delta\mathsf{\varepsilon}$ is the known linear strain increment calculated from the load increment. In these equations all values are known, and we can check the plasticity condition (7.37). First, we calculate

If $f_2^{test}\le0$, then the test values are valid for this (elastic) load step; otherwise, these values must be returned to the changed plasticity area using the plastic corrector.

7.5 Equations for the return of the stress point on the yield area (plastic corrector)

The first equation, which the unknown stress $\mathsf{\sigma}_{n+1}$ have to fulfill at a plastic load step, is the condition (7.37) at the time $t_{n+1}$

The equivalent plastic deformation in the load step will change by the increment (7.40) to

and the last equation in this system of nonlinear equations is the implicit relationship of the backward Euler method for calculation of the elastic deformation

After a suitable arrangement the system of equations to can be changed to only one nonlinear equation with only one unknown $\Delta\lambda$. Starting with the elastic constitutive equations [1]

gradually we get

With the known $\Delta\lambda$ it is possible, according to (7.47), to calculate the stress components from

or, because $\mathsf{\sigma}_{n+1}^{test}=\mathbf{D}^e(\mathsf{\varepsilon}_{n+1}-\mathsf{\varepsilon}_n^p)$, from

where

7.6 Specifying the plastic multiplier from the consistency condition

There is only one unknown in the above equations, namely the plastic multiplier $\Delta\lambda$. Its value must fulfill the consistency condition (7.43). Direct substitution of the matrix relations for $\mathsf{\sigma}_{n+1}$ (7.49) under this condition leads to an unpleasant matrix form of a nonlinear equation; therefore, it is better to use an explicit form, the derivation of which is presented in this section. In the case of isotropic elastic material, the matrices $\mathbf{D}^e$ and $\mathbf{P}$ have the same eigenvectors, and it is possible to change them to diagonal by their orthogonal transformation in the plasticity condition, and thus obtain an explicit equation form suitable for calculation. The transformation relations are

where

and $\mathbf{D}^e$ is given by (7.10). From we get

where

Substituting (7.53) into the consistency condition (7.43) gives the simplified explicit equation

where

with

For the yield stress at the end of the load step we have

The last problem is to calculate $\Delta\lambda$ from the nonlinear equation (7.55).

7.7 Determination of the plastic multiplier by the Newton-Raphson method

We are searching for the value of the plastic multiplier at which the equation (7.55) will be fulfilled with sufficient accuracy. In other words at the Newton-Raphson method the value of $\Delta\lambda$ has to change in such way that function

converges to zero. The classical iterative formula for an equation with one unknown in our case gives

where

The partial derivatives needed for $q'$ can be determined from (7.56) and (7.58)

7.8 Matrix relations in Fortran

In the FEM programs above relationships are written in the form of a subroutine which is called in the cycle for each integration element point of the body model. To conclude this section, we show such a subroutine in the form of an executable program. The program is based on the above relations and should be comprehensible for people with elementary knowledge of FORTRAN language. The program solves the same example as that solved above with other programs (examples 1 and 2).

Results (Napätia = Stresses, Medza sklzu = Yield stress):

7.9 Consistent tangent material modulus

The global tangential stiffness matrix of an FE model is created by an organized summation of its element tangential stiffness matrices

where, in a material nonlinear problem, the matrix $\mathbf{D}$ is the so-called consistent tangential material modulus of the element. The term consistent expresses a requirement that its value must be in every incremental step consistent with the integration procedure for the solution of the stress $\mathsf{\sigma}_{n+1}$. The definition of the tangential material modulus is

The derivations equality in (7.65) follows from validity

The consistent linearization of the local material element matrices produces the global tangent stiffness matrix in such a form, which guarantees the quadratic convergence of the global Newton-Raphson iteration in the load step.

The derivation of the $\mathbf{D}$ matrix begins with differentiation of expression for stress (7.46)

and after similar arrangements as at (7.46) we get

It remains to express $d\Delta\lambda$ in (7.68) from the consistency condition. Differentiation of the consistency condition (7.43) gives (note: $d(\mathbf{x}^T\mathbf{x})=2\mathbf{x}^Td\mathbf{x}$)

From (7.44) we have

where

By substituting (7.70) in (7.69) we have

and after substitute a we get a simplified equation

where

Substitution (7.68) into (7.73) gives the following relationship

and from this equation the wanted differential is given by

The use of $d\Delta\lambda$ in (7.68) gives the resulting expression for the stress increment

with the consistent tangential modulus

Thus, if an FEM program for solving the elastic-plastic plane stress problem uses the incremental calculation of the stress according to equations (7.48) or (7.49) , then to ensure the quadratic convergence of the Newton-Raphson method at the solution of the global nonlinear system, it is necessary to use the material modulus (7.78) at the tangential element matrices (7.64) calculation. Use of the classical tangential material modulus (so-called elastoplastic continuum tangent operator)

leads to a rapid decline in convergence.

7.10 Consistent material modulus computation in FORTRAN.

According to the relationships in the previous section, we processed the program MATMODUL2D for a consistent tangential material modulus with linear isotropic hardening in FORTRAN. In a real FEM program, such a part acts as a subroutine called for every integral point from the subroutine, which calculates the tangential element stiffness matrices. The input variables come to such subroutine as parameters of the calling command, but in this program illustration, we have input them directly in the program part noted as Input values.

The calculated consistent tangential modulus is as follows:

Table 1 presents the results of the iterative procedure of the test program, which illustrates the effectiveness of using the tangential material modulus consistent with the integration procedure (third column), where Fnorma is the norm of the residual forces and ITE is the number of iteration steps. The first two columns in the table indicate these values using the elastic material modulus $\mathbf{D}^e$ and the classical modulus of continuous integration $\mathbf{D}_{cont}$. During the integration procedures, the convergence criterion Fnorma should fall below 0.001.

It is seen that only the consistent tangent modulus preserves the quadratic convergence of the global equilibrium iteration.

[1] Benča Š.: Aplikovaná nelineárna mechanika kontinua, Vydavateľstvo 1000knih.sk, 2013
        (Applied nonlinear continuum mechanics – in Slovak).