Difference between revisions of "Contrib:KeesWouters/plate/variable pressure"

Revision as of 14:28, 13 January 2013

Varying pressure of shells - use of Python

If you are already familiar to Code Aster, enjoy this. If you are new to Code Aster have a look here first please.
A similar description for variable thickness of a plate is given in Contrib:KeesWouters/plate/thickness
november 2010 - SalomeMeca 2010 (Code Aster 10.2) - Ubuntu Meerkat 10.10 64bits.
January 2013 - update to Code Aster version 11.3 - Salome 6.6.0 - Linux Mint 14 Nadia

Static displacement of a shell under pressure

This is just a simple example of the use of shells (coque_3d) for static calculations. The main focus is on varying the pressure on a shell element by a Python list. The pressure on the shell depends on the average Ycog coordinate of the element. for mesh geometry, mesh and material property details see, once again, Contrib:KeesWouters/plate/thickness.

General way of determination of the pressure

The geometry and mesh are generated in Salome in the standard way. In Code Aster the mesh is read, including nodes number, node coordinates, element types and element connectivity. We use quadratic quad8 elements that are converted to quad9 (coque_3d) elements suitable for C-Aster.
Of each element

• determine the nodes are defined by the connectivity matrix,
• determine the coordinates of all the nodes and finally
• determine the centre of gravity and
• generate the thickness of the shell

by Python commands. This includes the y coordinate Ycog. The pressure depends on the centre y-value of each shell element. It is conveinient to use the added, last node of each shell element of the connectivity list. This node has been added in the Code-Aster command file by the MODI_MAILLAGE commnand to modify the QUAD8 to QUAD9 element. The geometry and mesh is shown in the picture below. The main dimensions are: [x1, x2; y1, y2; z1, z2] = [0.00,1.00; 0.00,3.00; 1.25,2.35], inclined plate in y direction.

Extraction of the coordinates

Among others, in Code Aster forum the procedure is discussed.
In Code Aster, after reading the initial mesh, the quad8 elements are converted to quad9 (coque_3d compatible) elements. This mesh is called meshmod (modified mesh) and is used as the basis for the further operations. Note that the last node in the connectivity list of the element is the added, centre node that will be used for applying the pressure value.

First we adapt the DEBUT() command since we need real Python commands:

DEBUT(PAR_LOT='NON');

We import the Python partition module in the command file to read mesh quantities. The Python file for determining the pressure on the shells is a Python script/module/file is ..... .

from Utilitai.partition import *
import sys
sys.path.append('/cae/caexample/shell/shellplate_thlist')   ## define your own path here
from construct_shell_thickness import *

After reading the initial mesh meshinit with quad8 elements

meshinit=LIRE_MAILLAGE(FORMAT='MED',INFO=1,);

and converting it to quad9, coque_3d elements

meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
                     MODI_MAILLE=(_F(TOUT='OUI',
                                     OPTION='QUAD8_9',),),);

info = 0 PressureOnShell = ApplyPressureOnGroup(modiMesh,'shell',info)

ApplyPr = AFFE_CHAR_MECA(MODELE=cocModel,FORCE_COQUE=PressureOnShell);

we define the pressre on each element by calling the following command sequence:

info = 1   ## 0:1:2 - print additional information 0:none 1:general 2:detailed information
PressureOnShell = ApplyPressureOnGroup(modiMesh,'shell',info)

The Python function ApplyPressureOnGroup(mesh,groupstr,info) does the work.
Three arguments for the function need to be given now: the mesh, the group where the pressure needs to be applied shell and info for controlling printed information. (Note that the group shell has been added compared to the function for determining the thickness. This need to be adapted there also.)
We have a look at this in detail. Standard, applying the load, in this case the pressure on a shell is performed by the command AFFE_CHAR_MECA:

ApplyPr = AFFE_CHAR_MECA(MODELE=cocModel,FORCE_COQUE=_F(GROUP_MA='shell',PRES=pressure,),);

will be replaced by:

ApplyPr = AFFE_CHAR_MECA(MODELE=cocModel,FORCE_COQUE=PressureOnShell)

where the dictionary PressureOnShell has been filled by the Python function:

PressureOnShell = ApplyPressureOnGroup(modiMesh,'shell',info)

Now remember that the parameters in the _F(...) operator can be replaced by a Python list, in this case: ellist = {'GROUP_MA': 'shell','PRESS':1e5}. Or more general:

 ellist = [{'MAILLE': 'Mxxx','PRESS':Pxxx}, 
           {'MAILLE': 'Myyy','PRESS':Pyyy}, 
             .......                      
             .......
           {'MAILLE': 'Mzzz','PRESS':Pzzz}]

where 'Mxxx' is the element number and Pxxx is the corresponding pressure on the element. Eg in this example the following list thicknessGroupShell is generated:

 [{'EPAIS': 0.16625, 'MAILLE': 'M15'}, 
  {'EPAIS': 0.22875, 'MAILLE': 'M16'}, 
   ....
   ....
  {'EPAIS': 0.22875, 'MAILLE': 'M25'}, 
  {'EPAIS': 0.10375, 'MAILLE': 'M26'}]

General - extraction of coordinates

In this example, the actual thickness of the shell depends on the y coordinate at the centre of the element. This is extracted as follows. First read the mesh and extract the mesh by the commands:

meshCA   = MAIL_PY()
meshCA.FromAster(mesh)

In meshCA all the quantities can be extracted:

 nonu     = meshCA.dime_maillage[0]    ## number of nodes
 elnu     = meshCA.dime_maillage[2]    ## number of elements
 ncoord   = meshCA.cn                  ## coordinates of the nodes
 connect  = meshCA.co                  ## connectivity of the element
 NodeList = list(meshCA.correspondance_noeuds)   ## says it all
 ElemList = list(meshCA.correspondance_mailles)
 ElemType = list(meshCA.tm)

The procedure is then as follows:

• initialise the dictionary
• step through element list of the group GROUP_MA='shell'
  • read connectivity of the element, ie all attached node numbers: connectivity = meshCA.co
  • step through all node of the element, in case of QUAD9 -correct- 9 nodes
    • determine nodes and its coordinates
    • compute average (y) coordinate or select last, centre node element
    • define pressure on the shell depending on centre y coordinate
    • add the pressure and element to dictionary {...}