Abaqus to Code_Aster : plane stress

Here is a really simple comparison between Abaqus and Code_Aster syntax. This calculation can be found on this site.

The case here is a shell element with a point load in its plane :

femesh.cantilever

Expected results from Calculix (a free FEM software compatible with the Abaqus syntax) are :

Screenshot from 2017-03-19 21-09-47

After calculation, here are Code_Aster results :

cantilevern

As you can see, results are very similar but not identical. A possible reason is that shell elements are not as standard as beam elements in their behaviour. (the actual reason seems to be that Calculix translates the CPS4 elements into C3D8 hexa, causing locking as indicated here, thanks Pierre J. for explaining this!)

We can compare results at relevant nodes :

Node Abaqus/Calculix Disp Code_Aster Disp Diff %
N2 DX -1.88E-08 -1.68499301421E-08 10%
N2 DY -2.5188E-08 -2.31688649492E-08 8%
N3 DY -4.9401E-08 -4.7396662606E-08 4%
N5 DX -2.0527E-09 -1.66253225287E-09 19%
N5 DY -2.1897E-08 -2.08142601378E-08 5%
N6 DX -2.7304E-09 -1.54379242024E-09 43%
N6 DY -5.5466E-08 -5.24009052452E-08 6%
N8 DX 2.1645E-08 1.95674878507E-08 10%
N8 DY -2.0987E-08 -2.11438422919E-08 1%
N9 DX 3.4694E-08 3.04818811083E-08 12%
N9 DY -7.7596E-08 -7.15339782224E-08 8%

The translation between Abaqus and Code_Aster can be made by using the following correspondances :

Abaqus keyword Code_Aster keyword
*NODE See “cantilevern.mail” file: COOR_3D
*ELEMENT, TYPE=CPS4, ELSET See “cantilevern.mail” file: QUAD4 and GROUP_MA + AFFE_MODELE(MODELISATION=’DKT‘) + AFFE_CHAR_CINE to stay on XY plane
*SOLID SECTION, ELSET, MATERIAL AFFE_CARA_ELEM(COQUE) + AFFE_MATERIAU
*MATERIAL *DENSITY *ELASTIC, TYPE=ISO DEFI_MATERIAU(ELAS)
*BOUNDARY AFFE_CHAR_CINE
*STEP *STATIC MECA_STATIQUE or STAT_NON_LINE
*CLOAD AFFE_CHAR_MECA(FORCE_NODALE)
*NODE PRINT IMPR_RESU(FORMAT=’RESULTAT’)
*NODE FILE IMPR_RESU(FORMAT=’MED’)

The Abaqus input file looks like this :


*HEADING
 SQUARE PLATE SUBJECTED TO POINT LOAD
**
*NODE
       1,      0.0,        0.0
       2,      0.5,        0.0
       3,      1.0,        0.0
       4,      0.0,        0.5
       5,      0.5,        0.5
       6,      1.0,        0.5
       7,      0.0,        1.0
       8,      0.5,        1.0
       9,      1.0,        1.0
**
**
*ELEMENT, TYPE=CPS4, ELSET=SQUARE_P
       1,       1,       2,       5,       4
       2,       2,       3,       6,       5
       3,       4,       5,       8,       7
       4,       5,       6,       9,       8
**
** square_plate
**
*SOLID SECTION, ELSET=SQUARE_P, MATERIAL=STEEL
        0.05,
**
** steel
** Date: 26-Jul-94           Time: 14:14:09
**
*MATERIAL, NAME=STEEL
**
*DENSITY
        7680,
**
*ELASTIC, TYPE=ISO
     2.1E+11,         0.3
**
** fixed_edge
**
*BOUNDARY, OP=NEW
       1, 1,,          0.
       1, 2,,          0.
       4, 1,,          0.
       4, 2,,          0.
       7, 1,,          0.
       7, 2,,          0.
**
*RESTART,WRITE
**
** step  1,Default
**
*STEP, AMPLITUDE=RAMP
Linear Static Analysis
**
** This load case is the default load case that always appears
**
*STATIC
**
** force_1
**
*CLOAD, OP=NEW
       9, 2,       -100.
*NODE PRINT, FREQ=1
U
RF
*NODE FILE, FREQ=1
U
RF
*END STEP

As usual, we need separate files for Code_Aster, in particular a mesh (file “cantilevern.mail”) :



TITRE    
	cantilevern
	SQUARE PLATE SUBJECTED TO POINT LOAD
FINSF

COOR_3D

	N1  0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
	N2  5.0000000000e-01 0.0000000000e+00 0.0000000000e+00
	N3  1.0000000000e+00 0.0000000000e+00 0.0000000000e+00
	N4  0.0000000000e+00 5.0000000000e-01 0.0000000000e+00
	N5  5.0000000000e-01 5.0000000000e-01 0.0000000000e+00
	N6  1.0000000000e+00 5.0000000000e-01 0.0000000000e+00
	N7  0.0000000000e+00 1.0000000000e+00 0.0000000000e+00
	N8  5.0000000000e-01 1.0000000000e+00 0.0000000000e+00
	N9  1.0000000000e+00 1.0000000000e+00 0.0000000000e+00
FINSF


QUAD4 

 E1  N1 N2 N5 N4
 E2  N2 N3 N6 N5
 E3  N4 N5 N8 N7
 E4  N5 N6 N9 N8 
FINSF

GROUP_MA NOM = SQUARE_P
    % Original part : SQUARE_P
  E4 E1 E3 E2
FINSF

FIN

We also need the Code_Aster command file (file “cantilevern.comm”) :


# cantilevern
# SQUARE PLATE SUBJECTED TO POINT LOAD

DEBUT(PAR_LOT='NON',);

MAIL=LIRE_MAILLAGE(FORMAT='ASTER',
                   VERI_MAIL=_F(VERIF='OUI',),);

MODMECA=AFFE_MODELE(MAILLAGE=MAIL,
                    AFFE=_F(GROUP_MA='SQUARE_P',
                            PHENOMENE='MECANIQUE',
                            MODELISATION='DKT',),);

MAT00001=DEFI_MATERIAU(ELAS=_F(E=2.1e+11,
                               NU=0.3,
                               RHO=7680.0,),);

CHMAT=AFFE_MATERIAU(MAILLAGE=MAIL,
                    AFFE=_F(GROUP_MA='SQUARE_P',
                            MATER=MAT00001,),);

CINE0=AFFE_CHAR_CINE(MODELE=MODMECA,
                     MECA_IMPO=(_F(NOEUD=('N1', 'N4', 'N7'),
                                   DY=0.0,),
                                _F(NOEUD=('N1', 'N4', 'N7'),
                                   DX=0.0,),
                                _F(GROUP_MA='SQUARE_P',
                                   DZ=0.0,
                                   DRX=0.0,
                                   DRY=0.0,),),);

CHARGE1=AFFE_CHAR_MECA(MODELE=MODMECA,
                       FORCE_NODALE=_F(NOEUD='N9',
                                       FY=-100.0,),);

CAEL=AFFE_CARA_ELEM(MODELE=MODMECA,
                    COQUE=_F(GROUP_MA='SQUARE_P',
                             EPAIS=0.05,
                             VECTEUR=(0.9,0.1,0.2,),),);

RESU1=MECA_STATIQUE(MODELE=MODMECA,
                    CHAM_MATER=CHMAT,
                    CARA_ELEM=CAEL,
                    EXCIT=(_F(CHARGE=CHARGE1,),
                           _F(CHARGE=CINE0,),),
                    OPTION='SANS',
                    SOLVEUR=_F(NPREC=12,),);

TEST_RESU(RESU=(_F(RESULTAT=RESU1,
                   NUME_ORDRE=1,
                   NOM_CHAM='DEPL',
                   NOM_CMP='DY',
                   NOEUD='N9',
                   CRITERE='RELATIF',
                   VALE_CALC=-7.15339782224E-08,
                   PRECISION=0.1,
                   REFERENCE='SOURCE_EXTERNE',
                   VALE_REFE=-7.7596e-08,),),);

IMPR_RESU(FORMAT='RESULTAT',
          RESU=_F(RESULTAT=RESU1,
                  NOM_CHAM='DEPL',),);

IMPR_RESU(FORMAT='MED',
          UNITE=80,
          RESU=_F(MAILLAGE=MAIL,
                  RESULTAT=RESU1,
                  NOM_CHAM='DEPL',),);

FIN(FORMAT_HDF='OUI',);

You will also need a file to keep everything together ( file “cantilevern.export”). This file is the one needed to start the calculation on command line (“as_run cantilevern.export”) :


P actions make_etude
P memjob 262144
P memory_limit 256.0
P mode interactif
P mpi_nbcpu 1
P mpi_nbnoeud 1
P ncpus 1
P testlist verification sequential
P time_limit 60.0
P tpsjob 2
P version testing
A memjeveux 32.0
A tpmax 60.0
F comm cantilevern.comm D 1
F mail cantilevern.mail D 20
F mess cantilevern.mess R 6
F resu cantilevern.resu R 8
F rmed cantilevern.rmed R 81

8 thoughts on “Abaqus to Code_Aster : plane stress

  1. Hello Luca,

    You indicate:
    “As you can see, results are very similar but not identical. A possible reason is that shell elements are not as standard as beam elements in their behaviour.”

    You may already know it, but as indicated in CalculiX website:
    “The CPS4 and CPS4R elements are expanded into C3D8 and C3D8R elements, respectively.”
    http://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node41.html

    And C3D8 elements (8 node hexa element – full integration) have a well known tendency of locking.
    Would you translate them in Code_ASter with 8 node hexa with MODELISATION=3D, you would obtain closer results.

    You can see the locking effect: CalculiX results show less displacement vs Code_Aster, ie a stiffer structure. To reduce it, you can either try:
    – a much finer mesh.
    – or turning formulation to CPS4R

    Real Abaqus CPS4 elements would probably behave more like those of Code_Aster.

    Bests,

    Pierre

    PS: you don’t mention in the website Vega tools?

    Like

    1. Thank you very much Pierre, you are certainly right, I did not know about this difference between Calculix and Abaqus 🙂 I will try to do as you suggest by using MASSIF in Code_Aster, even if the objective of this post was to show how one can use Code_Aster towards Abaqus users, as the title suggest 😉 About Vega, a finite element converter that can be used to do these things automatically, it would have been a little surprise later : someone here would like to create and share a “Vega for Windows binary version” and asked me to write these two articles… Also, the actual version of Vega online cannot (yet) handle the Abaqus syntax, while one can add it rather easily, I think.

      Best Regards

      Like

  2. Hello Luca,

    I dare an additional small comment :).
    In your articles, when you compare results, the last column is systematically called “Err”, meaning “Error” I guess.
    In my opinion, you shouldn’t name it that way. One or the other (Abaqus or Code_Aster, or Nastran or Code_Aster) doesn’t have necessarily the right results.
    Naming it “Diff” for “Difference” would be more appropriate, don’t you think?

    Congratulation for this work.

    Best regards,
    Pierre

    Like

    1. Hello Pierre,
      I was thinking about it this morning! It is exactly as you say and I will change the column names. Also, I think that it would be nicer to have cases with a known mathematical solution, in this case we could show two “Errors” : Nastran and Code_Aster, against the “true” solution, I think that this is the right way to compare two FEM codes… I thought about two sources for these cases : the NASA inp files in the Nastran95 package and the Code_Aster validation cases, translated into Nastran or Abaqus, what do you think about that 🙂 ?

      Like

Leave a comment