Homogenization of non-linear materials

The mechanical properties of structures composed of materials with different elasto-plastic properties can be estimated based on the models of Taylor or Sachs. In this tutorial, structures composed of different materials are created and their mechanical properties are investigated and compared to the Taylor and Sachs solution. Non-linear Finite Element Analysis is performed with the pyLabFEA package, see the online documentation and the tutorial Introduction for detailed information on the functionality of the package. For linear-elastic materials, analytical solutions for the homogenzition of their properties exist, as described in the tutorial Composites. This tutorial uses the matplotlib (https://matplotlib.org/) library for the visualization of results, and NumPy (http://www.numpy.org) for mathematical operations with arrays.

Author: Alexander Hartmaier, ICAMS / Ruhr-Universität Bochum, Germany

March 2020

This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License (CC-BY-NC-SA) Creative Commons License

The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)

In [1]:
import pylabfea as FE
import matplotlib.pyplot as plt
import numpy as np

1. Create elastic-plastic materials

In the first step, the class Material in invoked to create objects for three materials, with identical elastic properties, but different yield strengths. The unit of stress is 1 MPa. All features of the class Material are described in the documentation.

In [2]:
mat1 = FE.Material(name="Soft")      # invoke class to generate object for material 1
mat1.elasticity(E=200.e3, nu=0.3)    # define material with isotropic elasticity
mat1.plasticity(sy=180., khard=15.e3) # define material with isotropic plasticity and linear work hardening
mat2 = FE.Material(name="Intermediate")  # define second material
mat2.elasticity(E=200.e3, nu=0.3)    # identic elastic properties as mat1
mat2.plasticity(sy=190., khard=15.e3) # higher yield strength than mat1, but same w.h. coefficient
mat3 = FE.Material(name="Strong")        # define third material
mat3.elasticity(E=200.e3, nu=0.3)    # identic elastic properties as mat1
mat3.plasticity(sy=230., khard=15.e3) # higher yield strength than mat1 and mat2, but same w.h. coefficient
'Calculate and plot stress strain curves of materials'
mat1.calc_properties(eps=0.01)  # invoke method to calculate properties of material 
mat2.calc_properties(eps=0.01)  # and to store stress-strain data up to total strain eps=1%
mat3.calc_properties(eps=0.01)
plt.plot(mat1.prop['sty']['eeq']*100., mat1.prop['stx']['seq'], '-b')  # plot equiv. strain (eeq) vs. equiv. stress (seq)
plt.plot(mat2.prop['sty']['eeq']*100., mat2.prop['stx']['seq'], '-m')  # for uniaxial tensile stress in vertical direction (sty)
plt.plot(mat3.prop['sty']['eeq']*100., mat3.prop['stx']['seq'], '-r')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{eq}$ (%)')
plt.ylabel(r'$\sigma_{eq}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name], loc='lower right')
plt.show()

2. Laminate structure under iso-strain condition (Taylor model)

The class Model is invoked to generate an object for FEA. The methods of the class are used to generate the geomtry, to assign materials to sections, to establish boundary conditions, and to create the mesh. After these pre-processing steps, the solver is called to determine the solution for the deformed shape of the laminate structure in mechanical equilibrium under the applied boundary conditions. In the post-processing step, the results are evaluated and visualized. See the tutorial Introduction for the basic steps; all features of the class Model are described in the documentation.

The FE model of the laminate structure is subjected to a uniaxial stress along the lamellae of the structure, leading to iso-strain conditions, in which each section of the model, i.e. each lamella, undergoes the same total strain as the entire structure. The unit of length is 1 mm. In the first step, the model is loaded up to a stress slightly higher than the yield strength of mat1, and the strains are analyzed.

In [3]:
'laminate model is generated and elastic-plastic properties are assigned to each section'
fem = FE.Model(dim=2, planestress=True) # generate an object for a 2d finite element model with plane stress conditions
fem.geom([2, 2, 2], LY=4.)      # define 3 horizontal sections with 2 micron width per section, vertical dimension is 4 micron
fem.assign([mat1, mat2, mat3])  # assign the proper material to each section
fmat1 = 1./3.                   # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress along the  direction of the lamellae'
fem.bcleft(0.)                     # fix left and bottom boundary in normal directions
fem.bcbot(0.)                      # nodes on lhs boundary can move up and down, bottom nodes can move left and right
fem.bcright(0., 'force')           # free boundary condition on right edge of plane stress model -> uniaxial stress condition
ubc = 1.05*fem.leny*mat1.sy/mat1.E # displacement on top boundary corresponding to yield strength of mat1 (+5%)
fem.bctop(ubc, 'disp') # displacement applied to top nodes
'generate mesh'
fem.mesh(NX=6, NY=2)   # create structured mesh
'solution of equations for mechanical equilibrium'
fem.solve()            # solve system of equations
'post-processing'
fem.calc_global()      # calculate global stress and strain

print('2-d Model: iso-strain, uniaxial stress')
print('--------------------------------------')
print('Global total strain: (eps_11, eps_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain2')                   # plot eps_22
fem.plot('peeq',vmin=0.,vmax=0.0015)  # plot eqiv. plastic strain

2-d Model: iso-strain, uniaxial stress
--------------------------------------
Global total strain: (eps_11, eps_22) = (-0.000300,  0.000900)
Local total strain (element solution in different sections, Voigt tensor)
[-0.0003  0.0009 -0.0003  0.      0.     -0.    ] in Section 1, Material:  Soft
[-0.00028  0.00094 -0.00028  0.       0.      -0.     ] in Section 2, Material:  Intermediate
[-0.00028  0.00094 -0.00028  0.       0.       0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = (-0.000000,  0.000000)
Local plastic strain (element solution in different sections, Voigt tensor):
[-1.5e-05  3.1e-05 -1.5e-05  0.0e+00  0.0e+00  0.0e+00]  in Section 1, Material: Soft
[0. 0. 0. 0. 0. 0.]  in Section 2, Material: Intermediate
[0. 0. 0. 0. 0. 0.]  in Section 3, Material: Strong

It is seen that the isostrain assumption is fulfilled to a very good degree, although the first section of the model, to which the soft material mat1 is assigned starts to yield plastically, while the other sections remain linear-elastic.

The analysis of the stresses reveals that the stress in section 1 remains at the level of the flow stress in mat1 (=yield strength + work hardening), while the stress in the two other sections is below the yield strength of those materials.

In [4]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)

fem.plot('seq')

Global stress (MPa): (sig_11, sig_22) = ( 0.00, 186.34)
Local stress (element solution in different sections) (MPa):
[  0.   181.11   0.     0.     0.    -0.  ]  in Section 1, Material: Soft
[  0.   188.95   0.     0.     0.    -0.  ]  in Section 2, Material: Intermediate
[  0.   188.95   0.     0.     0.     0.  ]  in Section 3, Material: Strong

In the next step the load is increased to exceed the yield strength of mat2, and the strains are analysed again. Note that only the boundary conditions need to be changed on the existing model.

In [5]:
'increase load up to yield strength of mat2'
ubc = 1.05*fem.leny*mat2.sy/mat1.E   # displacement on top boundary corresponding to yield strength of mat2 (+5%)
fem.bctop(ubc, 'disp')  # displacement applied to top nodes
fem.solve()             # solve system of equations
fem.calc_global()       # calculate global stress and strain
print('2-d Model: isostrain, uniaxial stress')
print('-------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain2')
fem.plot('peeq',vmin=0.,vmax=0.008)

2-d Model: isostrain, uniaxial stress
-------------------------------------
Global total strain: (epl_11, epl_22) = (-0.000300,  0.001000)
Local total strain (element solution in different sections, Voigt tensor)
[-0.0003  0.001  -0.0003  0.      0.      0.    ] in Section 1, Material:  Soft
[-0.00031  0.001   -0.0003   0.       0.      -0.     ] in Section 2, Material:  Intermediate
[-0.0003  0.001  -0.0003  0.      0.      0.    ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = (-0.000000,  0.000000)
Local plastic strain (element solution in different sections, Voigt tensor):
[-3.6e-05  7.2e-05 -3.6e-05  0.0e+00  0.0e+00  0.0e+00]  in Section 1, Material: Soft
[-1.2e-05  2.4e-05 -1.2e-05  0.0e+00  0.0e+00  0.0e+00]  in Section 2, Material: Intermediate
[0. 0. 0. 0. 0. 0.]  in Section 3, Material: Strong

Again, the iso-strain assumption is fulfilled very well, because all sections individually are subject to the same vertical strain as the entire structure. The plastic strain in section 1 increases, and also in section 2 plastic yielding sets in.

In [6]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)

fem.plot('seq')

Global stress (MPa): (sig_11, sig_22) = ( 0.00, 190.64)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[  0.  181.9   0.    0.    0.   -0. ]  in Section 1, Material: Soft
[  0.   190.57   0.     0.     0.    -0.  ]  in Section 2, Material: Intermediate
[  0.   199.45   0.     0.     0.     0.  ]  in Section 3, Material: Strong

The stress in sections 1 and 2 is bounded by the flow stress of the assigned materials, the stress in section 3 is still below the yield strength of mat3.

Finally, the strain applied on the model is increased to $\epsilon_{22} = 1$%.

In [7]:
'increase load up to eps_22=1%'
ubc = 0.01*fem.leny     # displacement on top boundary corresponding to 1% total strain
fem.bctop(ubc, 'disp')  # displacement applied to top nodes
fem.solve()             # solve system of equations
fem.calc_global()       # calculate global stress and strain

print('2-d Model: isostrain, uniaxial stress')
print('-------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain2')
fem.plot('peeq',vmin=0.,vmax=0.9)

2-d Model: isostrain, uniaxial stress
-------------------------------------
Global total strain: (epl_11, epl_22) = (-0.004700,  0.010000)
Local total strain (element solution in different sections, Voigt tensor)
[-0.0047  0.01   -0.0023  0.      0.     -0.    ] in Section 1, Material:  Soft
[-0.00467  0.01    -0.00228  0.       0.       0.     ] in Section 2, Material:  Intermediate
[-0.00464  0.01    -0.0023   0.       0.       0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = (-0.003000,  0.007000)
Local plastic strain (element solution in different sections, Voigt tensor):
[-0.00346   0.006919 -0.00346   0.        0.        0.      ]  in Section 1, Material: Soft
[-0.003432  0.006863 -0.003432  0.        0.        0.      ]  in Section 2, Material: Intermediate
[-0.003353  0.006706 -0.003353  0.        0.        0.      ]  in Section 3, Material: Strong

The deformation of the laminate structure is now becoming visible, and the plastic strains in the different sections are rather similar, because most of the deformation is now of plastic nature, whereas the elastic contributions is quite small.

In [8]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)

fem.plot('seq')

Global stress (MPa): (sig_11, sig_22) = (-0.00, 335.06)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[ -0.   316.94   0.     0.     0.     0.  ]  in Section 1, Material: Soft
[ -0.   325.61   0.     0.     0.     0.  ]  in Section 2, Material: Intermediate
[ -0.   362.64   0.     0.     0.    -0.  ]  in Section 3, Material: Strong

The yield strength is now exceed in all sections, due to the linear work hardening with identical rate for each material, the difference in the yield strength is still reflected in the equivalent stress in each section, which corresponds to the flow stress of the assigned materials.

In the next step, the equiv. stress vs. equiv. strain diagram for the laminate structure is plotted together with those of the three materials.

In [9]:
'plot stress-strain curves of materials and laminate structure'
plt.plot(mat1.prop['sty']['eeq']*100., mat1.prop['stx']['seq'], '-b')
plt.plot(mat2.prop['sty']['eeq']*100., mat2.prop['stx']['seq'], '-m')
plt.plot(mat3.prop['sty']['eeq']*100., mat3.prop['stx']['seq'], '-r')
plt.plot(FE.eps_eq(fem.egl)*100., FE.sig_eq_j2(fem.sgl), '-k')  # calculate and plot equivalent values of global strain and stress
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{eq}$ (%)')
plt.ylabel(r'$\sigma_{eq}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name, 'laminate structure'], loc='lower right')
plt.show()

It is seen that the global stress for the laminate structure corresponds to the average of the stresses in the three sections, given by the flow stress of the assigned materials. According to the Taylor model, the stress average is calculated at a constant total strain $\epsilon_\mathrm{tot}$, which is equal to the total strain $\epsilon^{(i)}$ in each section $i$ (iso-strain assumption), thus

\begin{equation} \epsilon^{(i)} = \epsilon_\mathrm{tot} \hspace{2em} (i=1,2,3). \end{equation}

The global stress $\sigma_0$ is obtained by averaging the stresses $\sigma^{(i)}$ of the section according to their volume fractions $f_i$ , which yields

\begin{equation} \sigma_0 = \sum\limits_{i=1}^3 f_i \, \sigma^{(i)}(\epsilon_\mathrm{tot}). \end{equation}

In the next cell, the stresses according to the Taylor model are calculated and plotted together with the FE solution. To accomplish this, the FE model is restarted an the element solution in each section is recorded for small load increments until the total strain is reached.

In [10]:
'initialize arrays for element stresses and strains'
sig1=[]
sig2=[]
sig3=[]
etot=[]
'reset model to start with zero load; geometry, sections, boundary conditions and mesh are kept'
fem.u=None # reset model
'calculate FE solution for small load increments'
for eps in [8.e-4, 1.e-3, 1.3e-3, 2.e-3, 4.e-3, 6.e-3, 8.e-3, 1.e-2]:
    ubc = eps*fem.leny      # calculate load increment
    fem.bctop(ubc, 'disp')  # apply displacement on top nodes
    fem.solve()             # solve system of equations
    sig1.append(fem.element[0].sig[1])  # store element solution for each section
    sig2.append(fem.element[4].sig[1])
    sig3.append(fem.element[8].sig[1])
    etot.append(fem.element[0].eps[1]*100)

'calculate weigthed average of stresses in different sections at constant total strains'
s_taylor = fmat1*np.array(sig1) + fmat2*np.array(sig2) + fmat3*np.array(sig3)

'plot stresses and strains'
plt.plot(etot, sig1, '-b')
plt.plot(etot, sig2, '-m')
plt.plot(etot, sig3, '-r')
plt.plot(fem.egl[:,1]*100., fem.sgl[:,1], '-k')  # plot vertical components of global strain and stress
plt.plot(etot, s_taylor, 'oc')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{22}$ (%)')
plt.ylabel(r'$\sigma_{22}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name, 'FE solution', 'Taylor model'], loc='lower right')
plt.show()

It is seen that the results of the Taylor model and the FE solution are in excellent agreement. Furthermore, the element solution for the uniaxial stress state corresponds to the stress-strain curves obtained for the individual materials.

3. Laminate model under iso-stress conditions (Sachs model)

In the next step, a new model with the identical geometry and materials is generated, but subjected to a uniaxial strain in the direction transversal to lamellae. This load case results in an iso-stress condition that can be described with the Sachs model.

In [11]:
'mechanical boundary conditions for Sachs model (iso-stress) with uniaxial strain'
fem = FE.Model(dim=2, planestress=False) # generate an object for a 2d finite element model with plane strain conditions
fem.geom([2, 2, 2], LY=4.)      # define 3 horizontal sections with 2 micron width per section, vertical dimension is 4 micron
fem.assign([mat1, mat2, mat3])  # assign the proper material to each section
fmat1 = 1./3.                   # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress in transversal direction'
fem.bcleft(0.)                  # fix left and bottom boundary in normal directions
fem.bcbot(0.)                   # nodes on lhs boundary can move up and down, bottom nodes can move left and right
ubc = 0.0015*fem.lenx           # apply uniaxial strain eps_11 = 0.15%
fem.bcright(ubc, 'disp')        # apply displacement on rhs nodes
fem.bctop(0., 'disp')           # top boundary is fixed for uniaxial strain
'generate mesh'
fem.mesh(NX=6, NY=2)            # create structured mesh
'solution of equations for mechanical equilibrium'
fem.solve()                     # solve system of equations
'post-processing'
fem.calc_global()               # calculate global stress and strain
print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('strain1')
fem.plot('strain2')
fem.plot('peeq',vmin=0., vmax=0.02)

2-d Model: iso-stress, uniaxial strain
--------------------------------------
Global total strain: (epl_11, epl_22) = ( 0.001500,  0.000000)
Local total strain (element solution in different sections, Voigt tensor)
[ 0.0016 -0.      0.      0.      0.     -0.    ] in Section 1, Material:  Soft
[0.00152 0.      0.      0.      0.      0.     ] in Section 2, Material:  Intermediate
[0.00142 0.      0.      0.      0.      0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = ( 0.000000, -0.000000)
Local plastic strain (element solution in different sections, Voigt tensor):
[ 1.75e-04 -8.70e-05 -8.70e-05  0.00e+00  0.00e+00  0.00e+00]  in Section 1, Material: Soft
[ 1.22e-04 -6.10e-05 -6.10e-05  0.00e+00  0.00e+00  0.00e+00]  in Section 2, Material: Intermediate
[0. 0. 0. 0. 0. 0.]  in Section 3, Material: Strong

It is seen that the uniaxial strain condition is applied and that sections 1 and 2 already underwent plastic deformation.

In [12]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
fem.plot('seq')

Global stress (MPa): (sig_11, sig_22) = (382.22, 183.85)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[382.22 198.3  198.3    0.     0.    -0.  ]  in Section 1, Material: Soft
[382.22 189.44 189.44   0.     0.     0.  ]  in Section 2, Material: Intermediate
[382.22 163.81 163.81   0.     0.     0.  ]  in Section 3, Material: Strong

Also the iso-stress condition in loading direction is fullfilled. Transversal stresses occur due to the restricted cross-contraction, they are different in each section, due to their different plastic deformation, which is reflected in the different equivalent stresses in each section. Note also, that the stress state has a significant triaxiality, because the stress component in x-direction is larger than the equivalent stress, which is independent of hydrostatic stress components.

Now, the strain is increased to $\epsilon_\mathrm{tot} = 1\%$, and the results are analyzed.

In [13]:
'increase load'
ubc = 0.01*fem.lenx
fem.bcright(ubc, 'disp') # displacement on rhs nodes
fem.solve()              # solve equations for mechanical equilibrium
fem.calc_global()        # calculate global stress and strain
print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('strain1')
fem.plot('peeq',vmin=0., vmax=0.6)

2-d Model: iso-stress, uniaxial strain
--------------------------------------
Global total strain: (epl_11, epl_22) = ( 0.010000,  0.000000)
Local total strain (element solution in different sections, Voigt tensor)
[ 0.0101 -0.      0.      0.      0.     -0.    ] in Section 1, Material:  Soft
[ 0.01004  0.       0.       0.       0.      -0.     ] in Section 2, Material:  Intermediate
[ 0.00989 -0.       0.       0.       0.      -0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = ( 0.005000, -0.003000)
Local plastic strain (element solution in different sections, Voigt tensor):
[ 0.005478 -0.002739 -0.002739  0.        0.        0.      ]  in Section 1, Material: Soft
[ 0.005426 -0.002713 -0.002713  0.        0.        0.      ]  in Section 2, Material: Intermediate
[ 0.005231 -0.002615 -0.002615  0.        0.        0.      ]  in Section 3, Material: Strong

Similar to the iso-strain conditions, the model is now fully plastified. Under uniaxial strain, however, the plastic deformation is smaller than for uniaxial stress, due to the pronounced stress triaxiality that is caused by the restricted cross-contraction.

Looking at the stresses reveals these strong components transversal to the loading direction.

In [14]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
fem.plot('stress2')
fem.plot('seq')

Global stress (MPa): (sig_11, sig_22) = (1858.54, 1570.73)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[1858.54 1588.63 1588.63    0.      0.      0.  ]  in Section 1, Material: Soft
[1858.54 1579.76 1579.76    0.      0.     -0.  ]  in Section 2, Material: Intermediate
[1858.54 1543.8  1543.8     0.      0.     -0.  ]  in Section 3, Material: Strong

The iso-stress condition in loading direction is fullfilled, i.e. the horizontal stress component is constant throughout the model. However, it is seen that due to the severe stress triaxiality, there are different equivalent stresses in the different sections, following the strengths of the assigned materials.

The homogenization under iso-stress conditions is covered by the Sachs model, where the stress in each section $\sigma^{(i)}$ equals the global stress $\sigma_0$, thus

\begin{equation} \sigma^{(i)} = \sigma_0 \hspace{2em} (i=1,2,3). \end{equation}

The global total strain is given as the weighted average of the total strain in each section $\epsilon_\mathrm{tot}^{(i)}$ at a constant stress $\sigma_0$ with the volume fractions $f_i$ as weights, which yields

\begin{equation} \epsilon_\mathrm{tot} = \sum\limits_{i=1}^3 f_i \, \epsilon_\mathrm{tot}^{(i)}(\sigma_0). \end{equation}

In the next cell, the strains according to the Sachs model are calculated from the element solutions in each section and plotted together with the global solution.

In [15]:
sig0=[]
seq1=[]
seq2=[]
seq3=[]
et1=[]
et2=[]
et3=[]
eq1=[]
eq2=[]
eq3=[]
fem.u=None # reset model
for eps in [1.e-3, 1.2e-3, 1.4e-3, 1.6e-3, 2.e-3, 3.e-3, 6.5e-3, 1.e-2]:
    ubc = eps*fem.lenx          # uniaxial strain increments
    fem.bcright(ubc, 'disp')    # apply displacement on rhs nodes
    fem.solve()                 # solve system of equations
    sig0.append(fem.element[0].sig[0])
    et1.append(fem.element[0].eps[0]*100)
    et2.append(fem.element[4].eps[0]*100)
    et3.append(fem.element[8].eps[0]*100)
    seq1.append(FE.sig_eq_j2(fem.element[0].sig))
    seq2.append(FE.sig_eq_j2(fem.element[4].sig))
    seq3.append(FE.sig_eq_j2(fem.element[8].sig))
    eq1.append(FE.eps_eq(fem.element[0].eps)*100)
    eq2.append(FE.eps_eq(fem.element[4].eps)*100)
    eq3.append(FE.eps_eq(fem.element[8].eps)*100)

'weigthed average of strains in different sections at same stress level'
e_sachs = fmat1*np.array(et1) + fmat2*np.array(et2) + fmat3*np.array(et3)   
'plot stresses and strains'
plt.plot(et1, sig0, '-b', label=mat1.name)
plt.plot(et2, sig0, '-m', label=mat2.name)
plt.plot(et3, sig0, '-r', label=mat3.name)
plt.plot(fem.egl[:,0]*100, fem.sgl[:,0], '-k', label='FE solution')  # plot horizontal components of global strain and stress
plt.plot(e_sachs, sig0, 'oc', label='Sachs model')
plt.xlim((0,0.3))
plt.ylim((0,700))
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{11}$ (%)')
plt.ylabel(r'$\sigma_{11}$ (MPa)')
plt.legend(loc='lower right')
plt.show()

The numerical solution and the Sachs model are in excellent agreement. Due to the large stress triaxility, the values of the flow stress are quite high, such that the differences between the different sections, caused by the comparatively small differences in the yield strength of each material, appear to be small in comparison. Thus, the evolution of equivalent stresses and equivalent strains in the elements is also evaluated.

In [16]:
'weigthed average of equiv. strains and stresses in different sections'
eeq_sachs = fmat1*np.array(eq1) +  fmat2*np.array(eq2) +  fmat3*np.array(eq3)
seq_sachs = fmat1*np.array(seq1) + fmat2*np.array(seq2) + fmat3*np.array(seq3)
'plot stresses and strains'
plt.plot(eq1, seq1, '-b', label=mat1.name)
plt.plot(eq2, seq2, '-m', label=mat2.name)
plt.plot(eq3, seq3, '-r', label=mat3.name)
plt.plot(FE.eps_eq(fem.egl)*100., FE.sig_eq_j2(fem.sgl), '-k', label='FE solution')  # plot horizontal components of global strain and stress
plt.plot(eeq_sachs, seq_sachs, 'oc', label='Sachs model')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{eq}$ (%)')
plt.ylabel(r'$\sigma_{eq}^\mathrm{J2}$ (MPa)')
plt.legend(loc='lower right')
plt.show()

The analysis of eqivalent stresses and strains makes the influence of the material properties more visible. In the next plot, Taylor and Sachs solutions are compared.

In [17]:
plt.plot(etot, s_taylor, '-b', label='Taylor model')
plt.plot(eeq_sachs, seq_sachs, '-r', label='Sachs model')
plt.title('Taylor vs. Sachs model')
plt.xlabel(r'$\epsilon_\mathrm{tot}$ (%)')
plt.ylabel(r'$\sigma_0$ (MPa)')
plt.legend(loc='lower right')
plt.ylim((0,350))
plt.show()

The Taylor model (iso-strain) results in a stronger material behavior than the Sachs model (iso-stress), in which the softer regions dominate the behavior of the entire structure.

Iso-stress approximation applied to uniaxial stress loading

In the next step, we investigate how applying a uni-axial stress transversal to the laminate structure changes the results. The number of elements is increased to accomodate the inhomogeneous transverse strains in a better way.

In [18]:
'laminate model is generated and elastic-plastic properties are assigned to each section'
fem = FE.Model(dim=2, planestress=True) # generate an object for a 2d finite element model with plane stress conditions
fem.geom([2, 2, 2], LY=4.)      # define 3 horizontal sections with 2 micron width per section, vertical dimension is 4 micron
fem.assign([mat1, mat2, mat3])  # assign the proper material to each section
fmat1 = 1./3.                   # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress in transversal direction'
fem.bcleft(0.)                     # fix left and bottom boundary in normal directions
fem.bcbot(0.)                      # nodes on lhs boundary can move up and down, bottom nodes can move left and right
fem.bctop(0., 'force')             # free boundary condition on top edge of plane stress model -> uniaxial stress condition
ubc = 1.05*fem.lenx*mat1.sy/mat1.E # displacement on rhs boundary corresponding to yield strength of mat1 (+5%)
fem.bcright(ubc, 'disp')           # displacement applied to rhs nodes
'generate mesh'
ny = 9
fem.mesh(NX=9, NY=ny)   # create structured mesh with higer number of elements in verical direction
'solution of equations for mechanical equilibrium'
fem.solve(verb=False)   # solve system of equations
'post-processing'
fem.calc_global()  # calculate global stress and strain
iel2 = 3*ny        # first element in section 2
iel3 = 6*ny        # first element in section 3

print('2-d Model: transversal uniaxial stress')
print('--------------------------------------')
print('Global total strain: (eps_11, eps_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[iel2].eps.round(decimals=5),'in Section 2, Material: ', fem.element[iel2].Mat.name)
print(fem.element[iel3].eps.round(decimals=5),'in Section 3, Material: ', fem.element[iel3].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[iel2].epl.round(decimals=6),' in Section 2, Material:',fem.element[iel2].Mat.name)
print(fem.element[iel3].epl.round(decimals=6),' in Section 3, Material:',fem.element[iel3].Mat.name)
'create graphical output of model'
fem.plot('strain1', shownodes=False)                   # plot eps_11
fem.plot('strain2', shownodes=False)                   # plot eps_22
fem.plot('peeq', shownodes=False,vmin=0.,vmax=0.0015)  # plot eqiv. plastic strain

2-d Model: transversal uniaxial stress
--------------------------------------
Global total strain: (eps_11, eps_22) = ( 0.000900, -0.000300)
Local total strain (element solution in different sections, Voigt tensor)
[ 0.001  -0.0003 -0.0003  0.      0.      0.    ] in Section 1, Material:  Soft
[ 0.00092 -0.00029 -0.00027  0.       0.       0.     ] in Section 2, Material:  Intermediate
[ 0.00091 -0.00027 -0.00027  0.       0.       0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = ( 0.000000, -0.000000)
Local plastic strain (element solution in different sections, Voigt tensor):
[ 9.2e-05 -4.6e-05 -4.7e-05  0.0e+00  0.0e+00  0.0e+00]  in Section 1, Material: Soft
[0. 0. 0. 0. 0. 0.]  in Section 2, Material: Intermediate
[0. 0. 0. 0. 0. 0.]  in Section 3, Material: Strong

Plastic yielding starts again in section 1, in which also the total strain localises. This causes also a slight difference in the cross contraction of section 1 compared to the other sections.

In [19]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[iel2].sig.round(decimals=2),' in Section 2, Material:',fem.element[iel2].Mat.name)
print(fem.element[iel3].sig.round(decimals=2),' in Section 3, Material:',fem.element[iel3].Mat.name)
fem.plot('stress1', shownodes=False)
fem.plot('stress2', shownodes=False)

Global stress (MPa): (sig_11, sig_22) = (182.87,  0.00)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[1.8389e+02 3.1700e+00 0.0000e+00 0.0000e+00 0.0000e+00 7.0000e-02]  in Section 1, Material: Soft
[ 1.8243e+02 -4.1800e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.7000e-01]  in Section 2, Material: Intermediate
[ 1.8162e+02 -4.9000e-01  0.0000e+00  0.0000e+00  0.0000e+00  3.0000e-02]  in Section 3, Material: Strong

The different cross-contractions of the sections cause significant stress components transversal to the loading direction. hence, altough the applied boundary condition correspond to a uniaxial stress, within the structure, the stress state deviates from the applied one. Note that the stresses are highest at the x-axis, because there the the displacements in y-direction are constrained to zero.

In [20]:
'increase load above the yield strength of mat2'
ubc = 1.3*fem.lenx*mat2.sy/mat2.E
fem.bcright(ubc, 'disp') # displacement on rhs nodes
fem.solve()              # solve equations for mechanical equilibrium
fem.calc_global()        # calculate global stress and strain

print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[iel2].eps.round(decimals=5),'in Section 2, Material: ', fem.element[iel2].Mat.name)
print(fem.element[iel3].eps.round(decimals=5),'in Section 3, Material: ', fem.element[iel3].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[iel2].epl.round(decimals=6),' in Section 2, Material:',fem.element[iel2].Mat.name)
print(fem.element[iel3].epl.round(decimals=6),' in Section 3, Material:',fem.element[iel3].Mat.name)
fem.plot('peeq', shownodes=False, vmin=0., vmax=0.06)
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[iel2].sig.round(decimals=2),' in Section 2, Material:',fem.element[iel2].Mat.name)
print(fem.element[iel3].sig.round(decimals=2),' in Section 3, Material:',fem.element[iel3].Mat.name)
fem.plot('stress1', shownodes=False)

2-d Model: iso-stress, uniaxial strain
--------------------------------------
Global total strain: (epl_11, epl_22) = ( 0.001200, -0.000400)
Local total strain (element solution in different sections, Voigt tensor)
[ 0.0018 -0.0006 -0.0005  0.      0.      0.    ] in Section 1, Material:  Soft
[ 1.23e-03 -4.50e-04 -3.30e-04  0.00e+00  0.00e+00  1.00e-05] in Section 2, Material:  Intermediate
[ 0.00095 -0.00031 -0.00027  0.       0.       0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = ( 0.000000, -0.000000)
Local plastic strain (element solution in different sections, Voigt tensor):
[ 0.00066  -0.000312 -0.000348  0.        0.        0.      ]  in Section 1, Material: Soft
[ 0.000213 -0.000113 -0.0001    0.        0.        0.      ]  in Section 2, Material: Intermediate
[0. 0. 0. 0. 0. 0.]  in Section 3, Material: Strong

Global stress (MPa): (sig_11, sig_22) = (192.84,  0.00)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[1.9761e+02 8.7300e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.9000e-01]  in Section 1, Material: Soft
[192.18  -4.77   0.     0.     0.     0.88]  in Section 2, Material: Intermediate
[ 1.8744e+02 -5.9900e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.7000e-01]  in Section 3, Material: Strong

Upon further loading, the work hardening in section 1 causes the flow stress in this region to increase. Finally, the stress in the laminate structure exceeds the yield strength of the material in section 2. The iso-stress condition is fulfilled to some approximation in the structure, although the different sections reveal pronounced differences in the cross-contraction.

In [21]:
'increase load uo to 1% strain'
ubc = 0.01*fem.lenx
fem.bcright(ubc, 'disp') # displacement on rhs nodes
fem.solve()              # solve equations for mechanical equilibrium
fem.calc_global()        # calculate global stress and strain

print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[iel2].eps.round(decimals=5),'in Section 2, Material: ', fem.element[iel2].Mat.name)
print(fem.element[fem.Nel-ny].eps.round(decimals=5),'in Section 3, Material: ', fem.element[fem.Nel-ny].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)' 
      % (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[iel2].epl.round(decimals=6),' in Section 2, Material:',fem.element[iel2].Mat.name)
print(fem.element[iel3].epl.round(decimals=6),' in Section 3, Material:',fem.element[iel3].Mat.name)
fem.plot('peeq', shownodes=False)

2-d Model: iso-stress, uniaxial strain
--------------------------------------
Global total strain: (epl_11, epl_22) = ( 0.010000, -0.004900)
Local total strain (element solution in different sections, Voigt tensor)
[ 0.0117 -0.0055 -0.0026  0.      0.     -0.    ] in Section 1, Material:  Soft
[ 0.01045 -0.00516 -0.00226  0.       0.       0.     ] in Section 2, Material:  Intermediate
[ 0.00833 -0.00451 -0.00164  0.       0.       0.     ] in Section 3, Material:  Strong

Global plastic strain: (epl_11, epl_22) = ( 0.008000, -0.004000)
Local plastic strain (element solution in different sections, Voigt tensor):
[ 0.009291 -0.004615 -0.004676  0.        0.        0.      ]  in Section 1, Material: Soft
[ 0.008445 -0.003909 -0.004536  0.        0.        0.      ]  in Section 2, Material: Intermediate
[ 0.006526 -0.003791 -0.002735  0.        0.        0.      ]  in Section 3, Material: Strong

At the maximum strain, the laminate structure is again fully plastified. The plastic strains in the different sections reflects the different states of work hardening in the associated materials.

In [23]:
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)' 
      % (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[iel2].sig.round(decimals=2),' in Section 2, Material:',fem.element[iel2].Mat.name)
print(fem.element[iel3].sig.round(decimals=2),' in Section 3, Material:',fem.element[iel3].Mat.name)
fem.plot('stress1', mag=100, shownodes=False)

Global stress (MPa): (sig_11, sig_22) = (335.74,  0.00)
Local stress (element solution in different sections, Voigt tensor) (MPa):
[ 3.3943e+02 -4.7900e+00  0.0000e+00  0.0000e+00  0.0000e+00 -1.6000e-01]  in Section 1, Material: Soft
[3.4136e+02 1.5600e+01 0.0000e+00 0.0000e+00 0.0000e+00 3.2000e-01]  in Section 2, Material: Intermediate
[331.29 -14.98   0.     0.     0.     1.3 ]  in Section 3, Material: Strong

The iso-stress assumption is still fulfilled within an error margin of 10%, despite the different cross contractions that are magnified here with a factor of 100 to illustrate the shape under deformation fulfilling uniaxial stress conditions, i.e. with free relaxation of the cross-contractions.

4. Summary

This tutorial demonstrates the homogenization of structures composed of materials with different elastic-plastic properties by non-linear FEA. Models fhe different load cases leading to iso-strain and iso-stress conditions in the structures are set up and the results are analyzed. It is seen that the numerical solution agrees very well with the analyitcal homogenization rules according to the Taylor and Sachs models.

5. Exercises

  • Create materials with different hardening rates, e.g. 100 GPa for mat1, 50 GPa for mat2 and 10 GPa for mat3 and observe how the homogenized stress behaves in that case.
  • Create a model in which the sections are not equally distributed and verify the correctness of the Taylor and Sachs models in this case.
  • Assign a purely elastic material to one of the sections and study the behavior of the structure.