Difference between revisions of "Contrib:KeesWouters/staticshell"

Revision as of 14:51, 8 May 2009

Static behaviour of rectangular plate

The dynamics of a rectangular plate can be predicted quite well with shell (coque_3d) elements. I also liked to know the accuracy of the displacements for this kind of structures. A comparison between solid elements and shell elements has been made.

Dimensions of the plate The length and width of the plate are 1*1 m. The thickness of the shell elements has been increased to 0.02 m to enable meshing with solid elements (tetrahedrons and hexahedrons).

Boundary conditions and load Two types of boundary conditions have been given:

1. clamped at one side (all displacements and rotations are suppressed)
2. clamped at the four sides (all displacements and rotations are suppressed)

The load of the structure is a pressure of 0.1 bar on the surface of the plate

Material properties The The Young's module and poisson ratio for steel are set to 210 GPa and 0.28.

Displacement at the tip using beam theory

According to the shell theory the displacement of a beam the tip under uniform pressure is:

• w = q*L^4/(8*D)

• uniform pressure q = 10 kPa
• length of plate L = 1 m
• Youngs' module E = 210 GPa
• poisson ratio nu = 0.28
• bending stiffness D = Eh^3/(12*(1-nu*nu)) = 151910 Nm and
• displacement at the tip w = 8.2286 mm

Salome model of the rectangular plate

Shell elements

The geometry and meshing follows the same procedure as for the dynamics: geometry of the plate and generating the mesh.

Boundary conditions

For these calculations all the degree of freedoms (dofs) of one edge of the plate have been suppressed: three displacements and three rotations. This has been done in the corresponding comm file with the fillowing rule:

clamped=AFFE_CHAR_MECA(MODELE=modelc,
                       DDL_IMPO=(_F(GROUP_MA=('Linex0',),
                                   DX=0.0,
                                   DY=0.0,
                                   DZ=0.0,
                                   DRX=0.0,
                                   DRY=0.0,
                                   DRZ=0.0,),),
                        FORCE_COQUE=_F(GROUP_MA='shell',
                                   PRES=10.0e3,),);

Also the applied load has been given here in the last two lines. The group Linex0, of course, defines one edge of the plate.

For this load case the calculated displacement at the tip of the plate is

w_ca = 8.460 mm for a mesh with 40 nodes (3*3 quad8 elements) and
w_ca = 8.473 mm for a mesh with 2300 nodes.

In both cases this is slightly more than the theoretical value based on the shell theory: +2.8 % and +3.0 % for 40 and 2300 nodes respectively.

Comparison to solid elements

So it may be noted that the shell elements converge fairly rapidly to the correct value (here they even overestimate the displacement, can that be explained?). We expect that for solid elements convergence will be much slower.

In the graph below the results of some calculations are shown. We have used tetrahedron and hexahedron elements, both linear and quadratic.

On the x-axis the number of nodes (thousands of nodes) are given. They range from 40 for the shell mesh to nearly 150k for the solid elements. On the vertical axis the calculated displacement is given where the convergence value is about 0.0082 m.

Note that the 150k nodes is about the maximum my hardware is able to handle. It needs about 3.5 GB ram for this calculation and my physical limit is 4 GB now.

So it is fair to say that linear elements will be much too stiff for normal cases. Therefor the displacements of the structure will be predicted much too low. There has been a discussion about this many times on the CAElinux forums, but I didnot expect that linear elements give such a large stiffness compared to quadratic elements.

Mesh with solid elements

For comparison solid elements have been used to determine the behaviour of the plate. It is very easy in Salome to change from linear to quadratic elements so I in most cases the calculation has been carried out with linear and quadratic elements using the same mesh grid. Saving the mesh into a med file before and after conversion to quadratic elements is all that needs to be done. Starting the Code Aster calculations by ASTK is just hitting the run button if all file and directories are set.

Boundary conditions for the solid elements

Because the solid elements have no rotational degree of freedom the clamped edges must be given by fixed displacement of the side planes. The boundary conditions and loads are now defined in the comm file by:

clamp1=AFFE_CHAR_MECA(MODELE=ModQuad,
                      DDL_IMPO=(_F(GROUP_MA='f1',
                                   DX=0.0,
                                   DY=0.0,
                                   DZ=0.0,),),
                    PRES_REP=_F(GROUP_MA='surfbot',
                                   PRES=10000.0,),);

The group f1 is one of the side planes of the plate and all displacements are suppressed by the DDL_IMPO statement.

The complete command file is given here: Media:commandfile.comm.zip

Hexahedron elements

The following python file defines the geometry and the mesh of the plate. The mesh consists of hexahedron elements. The number of nodes and elements is largely controlled by the parameter Mfineness at the start of the file. Currently set to 0.15 it yields exactly 30k nodes. The number of layer along the thickness of the plate is 4. Setting different values for Mfineness (in the range from 0.0 to 0.6) yields between 2 and 12 layers. For large number of nodes only linear elements can be solved due to lack of memory on my hardware.

Media:geom_mesh_solidplate.zip

Tetrahedron elements

Now change the meshing part in the python file to:

elemVol = 0.12
TMesh = smesh.Mesh(Extrusion1) 
Tetrahedron_Netgen = TMesh.Tetrahedron(algo=smesh.NETGEN)
Max_Element_Volume = Tetrahedron_Netgen.MaxElementVolume(elemVol)
Max_Element_Volume.SetMaxElementVolume(elemVol)
smeshObj = smesh.smesh.CreateHypothesis('NETGEN_Parameters_2D','NETGENEngine')
Netgen_1D_2D = smesh.smesh.CreateHypothesis('NETGEN_2D','NETGENEngine')
status = TMesh.AddHypothesis(Netgen_1D_2D)
isDone = TMesh.Compute()

which basicly set the Netgen 3D and 2D parameters with the parameter elemVol. Setting elemVol=0.12 yields a linear mesh with 1.8k nodes with only one layer along the height of the plate, see picture.

Comparison of results

Finally the displacement at the tip for various numbers of nodes and the corresponding number of layers are given in the table below. This table has been used to compose the graph before.

mesh type	w_ca	        Rel %      KNodes   layers

shell quad	0.008460	102.81	   0.04	
shell quad	0.008473	102.96	   2.3	

tetra LIN	0.001049	12.75	   6	   1
tetra LIN	0.002515	30.57	   16	   2
tetra LIN	0.003300	40.10	   26	   2

tetra Quad	0.008423	102.37	   37	   1
tetra Quad	0.008459	102.81	   100	   2
tetra Quad	0.008461	102.83	   163	   2

hexa LIN	0.001984	24.11	   1	   2
hexa LIN	0.004738	57.58	   5 	   3
hexa LIN	0.005193	63.11	   8	   4
hexa LIN	0.006756	82.10	   24	   5
hexa LIN	0.007128	86.62	   37	   6
hexa LIN	0.007698	93.55	   90	   8
hexa LIN	0.007911	96.14	   156	  10

hexa Quad	0.008415	102.26	3.5	2
hexa Quad	0.008450	102.69	19	3
hexa Quad	0.008454	102.74	30	4
hexa Quad	0.008461	102.82	91	5
hexa Quad	0.008462	102.84	143	6

--kwarky 14:44, 7 May 2009 (CEST)