Contrib:KeesWouters/plasticity

From CAELinuxWiki
Revision as of 19:41, 7 February 2011 by Keeswouters (Talk | contribs) ('''V shaped construction under plastic deformation''')

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

V shaped construction under plastic deformation

To start with, this contribution mainly focuses on the use of Salome and Code Aster, not on the results and the mechanical justifications of the code that has been used. So no garantees that the results will be correct upto the fifth decimal. In fact, they wonot. I do hope though that this information is useful. For me it has been, because I had to think about some commands and look through the documentation and learn from that. In case of mistakes, errors and the like, please notify me, or better, you are invited to correct them yourself. Enjoy.

For version 10.X of Code Aster have a look here too

The construction is defined by an extrusion of a V shaped face.

Kw vshape1.jpg * Kw vshape2.jpg

Three surfaces are defined on the geometry: Pleft and Pright for applying boundary conditions (dx=dy=dz=0) and Pforce for defining a prescribed displacement in y direction (dy-->dy).

The definition of the geometry as well as the meshing is given in the Python input file.
In the geometry module of Salome select File --> Load Script (or ctrl T) and thereafter select the file: Media : kw_gm_vshape.zip.

General procedure of the calculation

The general procedure of the calculation is as follows:

  • read the 3D tetrahedral mesh, define models, meshes, groups of surfaces
  • apply the non linear material properties of the construction
  • define the boundary conditions
  • define the maximum load, in this case a maximum displacement at the plane Pforce
  • define 'time' and multiplication factor for 'load' increments
  • perform the non linear analysis
  • define results in a med-file for Salome, in this case, displacements, vonMises and plastic stresses
  • print results in a file, notably the displacements of the plane Pforce and its corresponding forces.


Kw res vshape.png

In the figure the results of four load cycles are given: the maximum displacement of each calculation is 0.080, 0.085, 0.095 and 0.100 mm. The plastic deformation at a load free construction after the load cycle is roughly 0.005, 0.009, 0.018 and 0.022 mm.
Each dot in the curves indicate a output point of the calculation. In most case 21 points are given, these have been defined in the 'time' increment function.

Non linear material properties

The linear material behaviour of steel is defined by the Young's modulus E = 210 GPa.

Kw sigma eps steel.png

The plastic material behaviour is defined by TRACTION and a table of (eps_i - sigma_i), as indicated in the figure. In the keyword DEFI_MATERIAU apart from Young's modulus and poisson ratio, now TRACTION is given for the non linear (eps-sigma) relation. The command AFFE_MATERIAU assigns the material to the complete mesh.

#Material properties
Traction=DEFI_FONCTION(NOM_PARA='EPSI',
                       VALE=(1.6e-3,300.0,
                             2.0e-2,400.0,
                             1.0e-1,450.0,
                             1.5e-1,473.0,
                             2.0e-1,500.0,
                             3.0e-1,550.0,),
                       INTERPOL='LIN',
                       PROL_DROITE='LINEAIRE',
                       PROL_GAUCHE='EXCLU',);
#define plastic behaviour of steel by Traction
steel=DEFI_MATERIAU(ELAS=_F(E=2.1e5,NU=0.27,),
                    TRACTION=_F(SIGM=Traction,),);
#assign material steel to the whole construction
matprops=AFFE_MATERIAU(MAILLAGE=pmesh,AFFE=_F(TOUT='OUI',MATER=steel,),);

Remarks from CA forum:

  • [1] Normally, the first point of a traction (tensile) curve is (eps=0.0,sigma=0.0). Note that in Code_Aster, when defining a traction curve with the command DEFI_FONCTION, the first point (0,0) must not be defined.


  • [2] stress-strain curve:
    • when inserting the material with elasto-plastic behaviour what kind of curve must be inserted? the true stress-true strain curve or engineering stress-strain curve? When I read the strain value this are logarithmic (true strain) or engineering?
  • It depends on your model and hypothesis on strains : small or large
    • if your model leads to small displacements and strains, then Cauchy stresses are equivalent to Poila Kirchoff stresses; the engineering curve is able to model the stress-strain curve up 5%.
    • if you have to take into account large displacements, but small strains, strain tensors are Green-Lagrange, and you have to define DEFORMATION=GREEN. But engineering stress-strain curve is still available
    • if you have to model large strains, then (with DEFORMATION='SIMO_MIEHE') stress-strain curve is the true one: Cauchy stress and logarithmic strains. The when you compute the stresses (with STAT_NON_LINE) they are Cauchy stresses, and the strains (with CALC_ELEM) are computes with the hypothesis of large strains.
  • About the definition of plastic zone with command DEFI_FONCTION, if you insert the epsilon and sigma value, the first point of the curve must be the strain value (for example about 0.002 for the steel) and second value the yield stress or must you insert the plastic residual value (ABAQUS SW use this formulation so the first point of strain must be 0 and the sigma value the sigma yeld value)?
    • The plastic curve is defined as follow: stress = function( total strain). That means the first abscissa corresponds to the elastic strain: epsilon = yield_stress/Emodulus

Boundary conditions and loads

In Code Aster boundary conditions and loads are treated in the shame way. So no distinction need to be made between them. In the mesh we defined three groups for applying boundary conditions and loads: for fixed boundary conditions Pleft and Pright. In the command file they are united to one group Pfix by the command:

pmesh=DEFI_GROUP(reuse =pmesh,
                   MAILLAGE=pmesh,
                   CREA_GROUP_MA=(_F(NOM='Pfix',UNION= ('Pleft','Pright',),),),);

Of course, this could have been done in Salome just as well.
The boundary condition itself are given by:

LoadFix=AFFE_CHAR_MECA(MODELE=pmode,
                   FACE_IMPO=(_F(GROUP_MA='Pfix',
                                 DX=0.0,DY=0.0,DZ=0.0,),),);

The load will be applied in steps defined by a 'time' function and multiplification factor on the load.

ydisp = 0.1000
LoadPres=AFFE_CHAR_MECA(MODELE=pmode,
                   FACE_IMPO=(_F(GROUP_MA='Pforce',DY=ydisp,),),);
tsteps = 21
time=DEFI_LIST_REEL(DEBUT=0.0,
                   INTERVALLE=_F(JUSQU_A=2.1,NOMBRE=tsteps,),
                   INFO=2,TITRE='time',);
ramp=DEFI_FONCTION(NOM_PARA='INST',
                  VALE=(0.00,0.00,
                        1.00,1.00,
                        1.10,1.00,
                        2.10,0.00,),
                  INFO=2,TITRE='ramp',);


  • Time:
the time is defined from 0.0 to 2.1 s, in 21 points or stepsize of 0.1 s
  • Load: the load multiplication factor is defined by
ramping from 0.0 to 1.0 during the first 1 s
a const level at 1.0 during 0.1 s and
decreasing from 1.0 to 0.0 in the last second.
Kw time ramp.jpg

At some step, e.g. step 14, the 'time' is 1.4 s and the corresponding multiplication factor for the load is 0.7, i.e. for the y displacement DY = 0.7*0.1000 = 0.070 mm.

Perform the calculation

The command for the static non linear case is STAT_NON_LINE

Presul=STAT_NON_LINE(MODELE=pmode,
                     CHAM_MATER=matprops,
                     EXCIT=(_F(CHARGE=LoadFix,),
                            _F(CHARGE=LoadPres,FONC_MULT=ramp,),),
                     COMP_INCR=_F(RELATION='VMIS_ISOT_TRAC',
                                  DEFORMATION='SIMO_MIEHE',
                                  TOUT='OUI',),
                     INCREMENT=_F(LIST_INST=time,
                                  SUBD_PAS=10,
                                  SUBD_PAS_MINI=0.005,),
                     NEWTON=_F(REAC_INCR=1,
                               MATRICE='TANGENTE',
                               REAC_ITER=1,),
                     CONVERGENCE=_F(ITER_GLOB_MAXI=20,),
                     ARCHIVAGE=_F(PAS_ARCH=1,),);

The keywords MODELE and CHAM_MATER define the model and material properties assigned in previous commands (not discussed here, see command file below).
The loadcase is defined by the keyword EXCIT and takes as parameters the boundary conditions and prescribed displacements of the plane Pforce. The multiplication factor is defined by the function FONC_MULT=ramp, where ramp is the previously discussed function.
The COMP_INCR defines the type of deformation (SIMO_MIEHO) and constitutive law are to be used.
In INCREMENT the previously defined 'time' function with the stepsizes must be given. When the solution does not converge using the original stepsize, possible substeps (SUBD_PAS) are defined. The SUBD_PAS_MINI (minimum stepsize) should be smaller than the original stepsize divided by SUBD_PAS. Eg. original stepsize is 0.1 s, SUBD_PAS=10, so SUBD_PAS_MINI=0.005 should be less then 0.1/10 = 0.01.
.....

Calculation with more substeps

In the example below the same calculation has been carried out with the courser mesh (maximum size is 5.0). Then the calculation obviously needs more substeps near the plastic region compared to the finer mesh.

Kw sigma eps steel 10steps.jpg

From step 7 to 11, i.e. from displacement 0.07 to 0.1 and back from 0.1 to 0.09 mm the time steps are divided into 10 substeps. During time step 10, 1.0 to 1.1 s, where the displacement remains 0.1 mm, the code needs only one step. See also result file at the end.

Results for post processing by Salome

If the calculation converges the results can be extracted. Here the displacements, von Mises stresses and plastic strains are stored for visualisation by Salome.

# Compute (VonMises) Stress / Strain
Presul=CALC_ELEM(reuse =Presul,
                 RESULTAT=Presul,
                 OPTION=('SIEF_ELNO_ELGA','EPSI_ELNO_DEPL','EPSP_ELNO','EQUI_ELNO_SIGM','EQUI_ELNO_EPSI',),);


# Write Results to MED file
# displacements DEPL (x, y and z component)
# vonMisses stresses EQUI_ELNO_SIGM (vonmisses component only)
# plastic strains EPSP_ELNO (six components epsxx, epsyy, ...)
IMPR_RESU(FORMAT='MED',
         UNITE=21,   ## change to 80 for default MED file extension ....
         RESU=(_F(MAILLAGE=pmesh,
                  RESULTAT=Presul,
                  NOM_CHAM='DEPL',
                  NOM_CMP=('DX','DY','DZ',),),
               _F(MAILLAGE=pmesh,
                  RESULTAT=Presul,
                  NOM_CHAM='EQUI_ELNO_SIGM',
                  NOM_CMP='VMIS',),
                _F(MAILLAGE=pmesh,
                  RESULTAT=Presul,
                  NOM_CHAM='EPSP_ELNO',
                  NOM_CMP=('EPXX','EPYY','EPZZ','EPXY','EPXZ','EPYZ',),),),);


Kw displ500.jpg
displacement at step 10 (maximum displacement)

Kw epspl500.jpg.jpg
plastic deformation (strain) at step 19 - roughly at the load/force free condition

Write displacements and nodal forces to table and file

Since we applied a displacement on the plane Pforce, which is easier to control the convergence, we are interested in the reaction forces on this plane.

First we determine the nodal forces and reaction on the nodes ('FORC_NODA','REAC_NODA')

# calculate nodal equivalents and reaction forces
Presul=CALC_NO(reuse =Presul,
            RESULTAT=Presul,
            OPTION= ('SIGM_NOEU_DEPL','EQUI_NOEU_SIGM','FORC_NODA','REAC_NODA',),);


Then define a group of nodes from the mesh plane Pforce

pmesh=DEFI_GROUP(reuse =pmesh,
                MAILLAGE=pmesh,
                CREA_GROUP_NO=_F(NOM='Nreac',GROUP_MA='Pforce',),);

Create a table for displacements ('DEPL'), reaction forces ('REAC_NODA') and nodal forces ('FORC_NODA'). In all three cases the three directions x, y and z (RESULTANTE=('DX','DY','DZ',)) will be extracted. Since no applied forces are present on this selected group of nodes, basically both fields are the same. ([REAC_NODA = FORC_NODA - applied load] or [REAC_NODA = FORC_NODA - chargements appliqués]).

TB_nodf=POST_RELEVE_T(ACTION=(_F(OPERATION='EXTRACTION',
                               INTITULE='displacements',
                               RESULTAT=Presul,
                               NOM_CHAM='DEPL',
                               TOUT_ORDRE='OUI',
                               GROUP_NO='Nforce',
                               RESULTANTE=('DX','DY','DZ',),
                               MOYE_NOEUD='NON',),
                            _F(OPERATION='EXTRACTION',
                               INTITULE='ReactionForces',
                               RESULTAT=Presul,
                               NOM_CHAM='REAC_NODA',
                               TOUT_ORDRE='OUI',
                               GROUP_NO='Nreac',
                               RESULTANTE=('DX','DY','DZ',),
                               MOYE_NOEUD='NON',),
                            _F(OPERATION='EXTRACTION',
                               INTITULE='NodalForces',
                               RESULTAT=Presul,
                               NOM_CHAM='FORC_NODA',
                               PRECISION=0.000001,
                               GROUP_NO='Nreac',
                               RESULTANTE=('DX','DY','DZ',),),),
                            TITRE='[Reaction] Nodal Forces',);

Finally, write the data to a file (define a file in ASTK with file unit 26, in this case).

IMPR_TABLE(TABLE=TB_nodf,
        FORMAT='TABLEAU',
        UNITE=26,
        SEPARATEUR=' * ',
        TITRE='displacements at nodes on group Nforce',);

The picture with the force - displacements at the start has been created with data of this table.

For more information regarding available fields, quantities and components that can be extracted from the results, see:

Download input and output files

Only suitable for CAster 9 and below: Media:kw_inout_vshape.zip, as detailed out above.

  • The file gm_beama.py defines geometry and mesh of the v-shape. After loading the script (cntrl T) into Salome, click refresh in the tree and save the vbeam mesh. For mesh size (MaxVolsize, line 93) I use 5.0 (coarser mesh, about 700 elements) and 1.0 (finer mesh, about 7k5 elements).
  • The file vbeam3.comm defines all the commands for Code Aster. You might play a bit with the settings of SUBD_PAS --> 20 and SUBD_PAS_MINI --> 0.0001 if the calculation does not converge. But in that case the mess file gives reasonable hints what to do with these parameters. Mentioning the mess file: always have a look at that, it gives useful information about progress and status of the calculation.
  • vbeam7k.astk is the astk input file.
  • The file vbdispl.txt is an example of the table defined by the last commands and can be used to construct the displacement - force graph.

Versions for Code Aster 10.X: Media:kw_vbeam_ca10.zip, slightly different syntax. See also


So far.
Any comments welcome.
augustus 2009