Difference between revisions of "Contrib:KeesWouters/shell/plotcoq3d"

From CAELinuxWiki
Jump to: navigation, search
m (Download the files)
m (Input files for the FE Analysis and references)
Line 181: Line 181:
  
 
Reference:<br/>
 
Reference:<br/>
[[http://www.code-aster.org/V2/doc/v9/man_u/u4/u4.81.04.pdf CALC_CHAMP]]
+
[[http://www.code-aster.org/V2/doc/v9/man_u/u4/u4.81.04.pdf CALC_CHAMP]] <br/>
[[http://www.code-aster.org/V2/doc/v9/man_u/u4/u4.81.05.pdf POST_CHAMP]]
+
[[http://www.code-aster.org/V2/doc/v9/man_u/u4/u4.81.05.pdf POST_CHAMP]] <br/>
 +
[[http://www.code-aster.org/forum2/viewtopic.php?id=18551]] <br/>
  
 
April 2013
 
April 2013
 
Enjoy - kees
 
Enjoy - kees

Revision as of 14:45, 28 April 2013

Post processing of COQUE-3D results

[under construction]

This contribution is only usefull for Code-Aster version 11.3 (maybe earlier version as 11.1 and 11.2 but certainly no previous releases).

In this contribution the main (and only) focus is on the postprocessing of the displacements and stresses of shell COQUE-3D elements. The centre nodes of tria7 and quad9 elements used for the discription of the displacements field of COQUE_3D poses some nasty problems in the postprocessing with Salome. To overcome this problem the results are projected on the original quadratic mesh without the centre node.

Definition of the geometry

This time the geometry has been drawn in FreeCAD. It is good enough at this moment (april 2013) to draw some more or less complex geometry, although, in this case, it wouldnot have been much more effort to construct it in Salome. Anyway, I import the brep geometry from FreeCAD into Salome and mesh it there. It is a planar shell with four holes. Main dimensions are Length x Width = 100 x 50. The holes with diameter of 6 are centred along the short centre line of have offset along the long centre line, see picture below. The thickness is defined in the Code Aster *.comm file and can be changed any time: thickness th = 0.02.
As usual no explicit dimensions are used here. But units have to be consistent and will be if you imagine:

  • length and displacments [m]
  • forces [N] and moments [Nm] and
  • Youngs' modulus and stresses [Pa]

Boundary conditions and loads:

  • all dof fixed on short side Lxm:
    • (DX=0.0, DY=0.0, DZ=0.0, DRX=0.0, DRY=0.0, DRZ=0.0)
  • 2 'loads': prescribed in-plane x displacement or out-out-of-plane vertical displacement fixed on short side Lxp:
    • (DZ=Zdisp) or
    • (DX=Xdisp)

Material:
Material steel

  • E = 210e9 Pa
  • Nu = 0.28
  • rho = 7850 kg/m3 (not needed)

Kw freecad sketch plane1.png * Kw coque3d mesh.png

A rough discription of the Code-Aster commands

The general flow is as follows:

  1. initialise
  2. import the initial mesh (initMesh)
  3. create a modified mesh (modiMesh)
    • create COQUE_3D elements with centre nodes
    • orient all element in the same direction
  4. creaste models with both meshes: iniModel and modModel
  5. define Coque_3d characteristics (thickness on group shell)
  6. define material (steel)
  7. aplly boundary conditions
  8. main calculation: RESU=MECA_STATIQUE(...)
  9. post process the result RESU:
    • create displacement and stress fields for modified model: CALC_CHAMP
    • create displacement and stress fields for modified model at bottom, centre and top layer: POST_CHAMP
    • create displacement and stress equivalent fields for modified model at bottom, centre and top layer: CALC_CHAMP
    • create displacement and stress fields for initial model at bottom, centre and top layer: PROJ_CHAMP. This last step is required to remove the data at the centre nodes that Salome is unable to coop with.
    • print results of modified and initial model (results of the modified model are only printed to show that they are not fully what you expect)

Steps 1 to 8 are fairly standard and are not descript in more detail here. We will concentrate on step 9). The stresses can be provided at the top, the centre and the bottom of the shell plane. These layers are denoted by 'SUP', 'MOY', and 'INF' in Code-Aster (superieur, moyenne et inferieur je pense) and depends on the normal vector of the shell elements (see step 3 in the general flow). These can be changed by the ORIENT command. In the next part all the commands could be given for the three layers, but it has been described here only for 'SUP' to keep it tidier. In the *.comm file at the end the complete command file can be downloaded with all layers excecuted.

Short description of CALC_CHAMP, POST_CHAMP and PROJ_CHAMP

We will use the commands CALC_CHAMP, POST_CHAMP and PROJ_CHAMP several times. Now a short description from the C-A-documentation is given:

  • Short description of CALC_CHAMP, from U4.81.04.pdf:
    • Create or supplement a result by computing fields by element or with the nodes (forced, strains, ...). The produced result concept either is created, or modified, i.e. the call to CALC_CHAMP is done in the following way:
      • resu = CALC_CHAMP(RESULTAT = resu..., reuse = resu,...), or
      • resu1 = CALC_CHAMP(RESULTAT = resu,...)
    • this last means that POST_CHAMP returns and extend its original field (resu) or create a new field (resu1). POST_CHAMP is not re-entrance.
  • Short description of POST_CHAMP, from U4.81.05.pdf:
    • specific Post-processing for the structural elements (shells, beams,...):
    • extraction of a field for a subpoint
    • computation of the minimum/maximum on all of the subpoints of a point taken into account of the eccentring of the plates for computation of the force
  • Short description of PROJ_CHAMP, from U4.72.05.pdf:
    • the goal of the operator is to project the fields of a data structure result on another mesh. This command can be used for example to transfer on a “mechanical” mesh, the result of a thermal computation carried out on a “thermal” mesh different. [In this case we want to get rid of the centre node of the coque-3d mesh that carries no displacements and just carry out a mechanical to mechanical projection.]

Calculate the mechanical behaviour of the structure

#main calculation
RESU=MECA_STATIQUE(MODELE=cocModel,
                    CHAM_MATER=material,
                    CARA_ELEM=shellch,
                    OPTION=('SIEF_ELGA'),   ##'SIEF_ELGA','DEPL','SICO_ELNO','SIGM_ELNO','SIEQ_ELNO'
                    EXCIT=(_F(CHARGE=clamped),
                           _F(CHARGE=ApplyPr),),);

Extracting the stresses of the coque_3d model

RESU=CALC_CHAMP(reuse =RESU,
               RESULTAT=RESU,
               CONTRAINTE=('SIEF_ELNO','EFGE_NOEU','EFGE_ELNO',),
               CRITERES='SIEQ_ELNO',
               EXCIT=(_F(CHARGE=clamped),
                      _F(CHARGE=ApplyPr),),);

Extracting the stresses of the coque_3d model on layers

Use the commands POST_CHAMP and again CALC_CHAMP to extract the local stresses Sxx, Syy, Szz, Sxy, Sxz and Syz and the equivalent stresses (VMIS, S1, S2, S3,TRESCA) from the result RESU issued in the previous command: SIEQ_SUP (equivalent stresses) and SIEF_SUP. POST_CHAMP only retrieves the local stresses (to be checked) and the subsequent CALC_CHAMP retrieves the equivalent stresses (to be checked). Note the all commands can be extended for the centre (MOY) and bottom (INF) layer in stead of the top layer (SUP). Finally the retrieved stresses and displacements are written to file (IMPR_RESU, all layers are choosen here, extend the POST_CHAMP and CALC_CHAMP for all layers if you want to use this).


SIEQ_SUP=POST_CHAMP(RESULTAT=RESU,
                   EXTR_COQUE=_F(NOM_CHAM='SIEQ_ELNO',
                                 NUME_COUCHE=1,
                                 NIVE_COUCHE='SUP',),);
SIEF_SUP=POST_CHAMP(RESULTAT=RESU,
                   EXTR_COQUE=_F(NOM_CHAM='SIEF_ELNO',
                                 NUME_COUCHE=1,
                                 NIVE_COUCHE='SUP',),);
SIEQ_SUP=CALC_CHAMP(reuse =SIEQ_SUP,
                   RESULTAT=SIEQ_SUP,
                   CRITERES='SIEQ_NOEU',);
SIEF_SUP=CALC_CHAMP(reuse =SIEF_SUP,
                   RESULTAT=SIEF_SUP,
                   CONTRAINTE='SIEF_NOEU',);
......
......
IMPR_RESU(FORMAT='MED',
         UNITE=80,
         RESU=(_F(RESULTAT=RESU,),
               _F(RESULTAT=SIEQ_SUP,),
               _F(RESULTAT=SIEQ_MOY,),
               _F(RESULTAT=SIEQ_INF,),
               _F(RESULTAT=SIEF_SUP,),
               _F(RESULTAT=SIEF_MOY,),
               _F(RESULTAT=SIEF_INF,),),);

Extracting the stresses of the initial model on layers

Alright, when using Salome the above results are disappointing (you may check the file written to unit 80 above, but be warned). To get really nice picture (and that is what FEM is all about, isnot it?) we have to get ride of the centre node of the mesh node. That is now fairly easy: project the results from the coque_3d model to the initial model by using PROJ_CHAMP. Then check the file written to unit 81, et voila, nice pictures.

QSUP_INI=PROJ_CHAMP(RESULTAT=SIEQ_SUP,
                   MODELE_1=cocModel,    # project fields of this model to
                   MODELE_2=iniModel,); # field of model_2
FSUP_INI=PROJ_CHAMP(RESULTAT=SIEF_SUP,
                   MODELE_1=cocModel,    # project fields of this model to
                   MODELE_2=iniModel,); # field of model_2
initRes=PROJ_CHAMP(RESULTAT=RESU,
                   MODELE_1=cocModel,    # project fields of this model to
                   MODELE_2=iniModel,); # field of model_2
IMPR_RESU(FORMAT='MED',
         UNITE=81,
         RESU=(_F(RESULTAT=initRes,NOM_CHAM=('DEPL'),),
               #_F(MAILLAGE=initMesh,RESULTAT=RESU,NOM_CHAM=('DEPL'),),
               _F(RESULTAT=QSUP_INI,),
               _F(RESULTAT=QMOY_INI,),
               _F(RESULTAT=QINF_INI,),
               _F(RESULTAT=FSUP_INI,),
               _F(RESULTAT=FMOY_INI,),
               _F(RESULTAT=FINF_INI,),),);
IMPR_RESU(FORMAT='MED',
         UNITE=82,
         RESU=(_F(RESULTAT=initRes),),);
FIN();

That's about it. Except for saying grace to TdS, jeanpierreaubry and Theobolt for asking and supplying help in the C-A-forum: [[1]]

Results of Ux displacement

The load is a displacement of 0.1 in x direction. The calculated Ux displacement (overlayed with the mesh in undeformed state) and the von Mises stresses in the centre layer are shown in the following pictures.
Kw dispux ux.png * Kw vnmises ux.png

Results of pressure load

In the picture below the vector displacement due to a pressure of 1 bar (1e5 Pa) is shown. Note that this is quite an unrealistic case, as the displacements are far too big. Kw plane uvect press.png

The next pictures show the von Mises and Sxx stresses in the top, centre and bottom layer of the shell. In the picture the centre stress is translated 55 in y direction and the bottom stress is translated 2*55 units. The von Mises stresses in the top and bottom layer are equal. In the centre layer the stress is practically zero. This is what we expect for a shell under pressure.
The stress component Sxx has different sign in the top layer compared to the bottom layer. Again, the centre layer stress is zero.

Kw plane vmises press.png

Kw plane sxx press.png


Download the files

Input files for the FE Analysis and references

Earlier version Code Aster 9&10
Input files:

  • Python script for defining the geometry and mesh (xxxxxx.py), load by File --> load script (cntrl T in the object browser), refresh (F5) after running
    • the mesh file should be save by filename Mshell by the previous script
  • ASTK file (xxx.astk, you need to edit the path to your requirements ...)
  • command file (xxx.comm)

Download the files here:

Reference:
[CALC_CHAMP]
[POST_CHAMP]
[[2]]

April 2013 Enjoy - kees