Topic: Adding new elements

Dear all,
is there a documentation about how to add new elements. I am working in shipbuilding and we are usign more complex elements (stiffened panels) for calculation. In addition does Quad Shaped Shell element exist?
Thanks in advance, Svemir.

2

Re: Adding new elements

Hi,

there is a simple, commented example of implementation of 2d truss element in oofem programmers manual (http://www.oofem.org/resources/doc/prog … ammer.html). It should give you enough information on how to implement any element. The description is slightly outdated, for example, the element defines the shape functions itself, while now we have interpolation classes for that, but in general it will give the oveerall picture.

There is no  Quad Shaped Shell element. What can help (as an example of implementation): quads for plane stress and strain and triangular shell element.

Re: Adding new elements

Hi,
we are trying to add new element but it seams that it will take more time than we expect (and have now). Is it possible to get more help? We are ready to pay for that kind of support. Reason is that we are looking for free FEM solver that we could include in our software (http://www.bvbcafe.com). Now it is running with NASTRAN, ANSYS and LS DYNA but more people are asking for free solver. Thanks in advance,Svemir.

Re: Adding new elements

Dear svemir,
Do you have any reference (a PDF preferably) of the exact type of shell element you wish to include?
Shell elements are always a bit special, due to locking behaviors and such.

Re: Adding new elements

Hi Mikael,
For test purpose it is enough to implement same 4 node shell elements as triangles which are already implemented in OOFEM (cct, cct3d,tr_shell01).

Re: Adding new elements

In that case, why do you need new elements? Exporting quads as triangles could be done directly in your software (and it could possibly be useful when using Abaqus or Ansys as well).
Although it's very simple, I'm reluctant to add a bundled element consisting of 2 triangles when it's so simple to just add two triangles in the input file.

Re: Adding new elements

Problem is what in Shipbuilding it is not allowed to use triangular elements. For example average coarse full ship model consists of 100000 elements (quads, beams, rods and triangles). Number of triangles is limited to few hundreds.

8

Re: Adding new elements

Svemir,

we could provide support, in terms of consultancy and help with the implementation. However, due to the limited resources, I could not commit myself to the implementation. On the other hand, I supervise one MSc student who will implement several shell elements in oofem in coming months. We will start with triangular element (as this is required by the project supporting that work), but later, some quad element(s) can be implemented as well.

Re: Adding new elements

Borek,
thank you for your information. Since I have to have answer on question which open source FEM we support until end of January waiting for a few months is not an option. I am trying to understand the way how OOFEM is organized. There is some progress in that direction. Complete code is very well written what helps me a lot. Once when I will understand it adding new elements and testing them will be easy. If we decide to support OOFEM then we will discuss about consultancy. Kind regards, Svemir.

P.S. Maybe Mikael will have some resources to help me in decission making.

Re: Adding new elements

Svemir,
If you could point me to a paper describing the element you which to be included, I could probably add it in a day.
If you don't you specifically know what theory you want for the element, then I would at least like to know;

  • Large or small displacements?

  • Is it enough for planar geometry, or do the shells need to be curved?

  • Do you need plastic material behavior, or just elastic structural elements?

If the answer is "small", "planar", and "elastic", then I'm quite confident it could be added quickly.

Re: Adding new elements

Mikael,
answer is "small" "planar" and "elastic" for this phase. I will be able to compare it with other packages having the same elements.
During the weekend I will work on it and let's see how far I will come. Kind regards, Svemir.

Re: Adding new elements

I read up on some shell elements, and it seems to me that a Mindlin plate + membrane element is probably what you are looking for .
With these simpler shell elements you end up with 5 degrees of freedom per node (just like rershell).
Problem is how to deal with the "missing" torsional DOF (which would likely cause problems in a naive attempt to solve in OOFEM).

I wonder, for the commercials tools, do you have to give special treatment to corners (6 dofs) and smooth/flat surfaces (5 dofs) manually?

Re: Adding new elements

6th dof is not a problem since our models are such that global stiffness matrix is always OK. In case of trivial examples user have to constraint 6th dof manually. Our tools cannot be used by dummies.

Re: Adding new elements

Svemir,
I've been looking at adding a simple Mindlin plate in 2D to start with. The descriptions I can find in litterature doesn't seem to match what CCT does (and CCT looks more reasonable). Hopefully, there will be some colleagues at the department this week so I can ask them for advice.

Then its a simple matter or adding a simple membrane stiffness, and use the local coordinate system to get the 3D material. After that it's just a question of implementing the needed boundary conditions (deadweight, normal pressure, etc.)
The type of element I had in mind is often called "MITC4".

Re: Adding new elements

Mikael,
MITC4 is good for starting point. That element is different in formulation then CCT that is why it is not so easy for me to add new one. I found more appropriate plane stress quad element as starting point.

Re: Adding new elements

Hello,
I am in process of adding CCQPlate Element. Prior to start testing it I was testing 1 m x 1 m steel plate loaded with force of -1000 N i vertical direction. Results that I am getting are strange. Translation in Z direction on position of Maximal load is smaller then along the sides of plate? Any idea about what can be wrong?

Input file is given below:

cct.out
Development Test of CCQ elements -> force in center of plate (plate 1x1 meter)
LinearStatic nsteps 1
domain 2dMindlinPlate
OutputManager tstep_all dofman_all element_all
ndofman 9 nelem 8 ncrosssect  1 nmat 1 nbc 2 nic 0 nltf 1
node 1 coords 3  -0.5  -0.5   0.0  bc 3 1 0 0
node 2 coords 3   0.0  -0.5   0.0  bc 3 0 0 0
node 3 coords 3   0.5  -0.5   0.0  bc 3 1 0 0
node 4 coords 3  -0.5   0.0   0.0  bc 3 0 0 0
node 5 coords 3   0.0   0.0   0.0  bc 3 0 0 0  load 1 2   
node 6 coords 3   0.5   0.0   0.0  bc 3 0 0 0
node 7 coords 3  -0.5   0.5   0.0  bc 3 1 0 0
node 8 coords 3   0.0   0.5   0.0  bc 3 0 0 0
node 9 coords 3   0.5   0.5   0.0  bc 3 1 0 0
CCTPlate 1 nodes 3 1 2 5 mat 1 NIP 1 crossSect 1
CCTPlate 2 nodes 3 2 3 6 mat 1 NIP 1 crossSect 1
CCTPlate 3 nodes 3 1 5 4 mat 1 NIP 1 crossSect 1
CCTPlate 4 nodes 3 2 6 5 mat 1 NIP 1 crossSect 1
CCTPlate 5 nodes 3 4 5 8 mat 1 NIP 1 crossSect 1
CCTPlate 6 nodes 3 5 6 9 mat 1 NIP 1 crossSect 1
CCTPlate 7 nodes 3 4 8 7 mat 1 NIP 1 crossSect 1
CCTPlate 8 nodes 3 5 9 8 mat 1 NIP 1 crossSect 1
SimpleCS 1 thick 0.01
IsoLE 1 d 0. E 2.1E11 n 0.3 tAlpha 0.000012
BoundaryCondition  1 loadTimeFunction 1 prescribedvalue 0
NodalLoad 2 loadTimeFunction 1 Components 3 0.0 0.0 -1000.0
ConstantFunction 1 f(t) 1.0

Outpu is given below:
############################################################## www.oofem.org ###
       ####    ####   ######  ######  #    #
     #    #  #    #  #       #       ##  ##
    #    #  #    #  #####   #####   # ## #
   #    #  #    #  #       #       #    #
  #    #  #    #  #       #       #    #   OOFEM ver. @VERSION@
  ####    ####   #       ######  #    #    Copyright (C) 1994-2008 Borek Patzak
################################################################################


Starting analysis on: Fri Jan 11 23:00:53 2013

development
Domain type: 2dmindlinplate, default ndofs per node is 3, per side is 0



==============================================================
Output for time  1.00000000e+000
==============================================================
Output for domain   1


DofManager output:
------------------
Node           1 (       1):
  dof 1   d  0.00000000e+000
  dof 2   d -7.43934900e-003
  dof 3   d  3.55983097e-003
Node           2 (       2):
  dof 1   d -6.36409779e-004
  dof 2   d  3.12716251e-003
  dof 3   d -2.71155296e-003
Node           3 (       3):
  dof 1   d  0.00000000e+000
  dof 2   d  7.90019768e-003
  dof 3   d -1.53367475e-003
Node           4 (       4):
  dof 1   d -1.92001511e-003
  dof 2   d  1.47350849e-003
  dof 3   d -1.20930012e-003
Node           5 (       5):
  dof 1   d -2.60967463e-016
  dof 2   d  1.12005251e-003
  dof 3   d -8.18373704e-003
Node           6 (       6):
  dof 1   d  1.92001511e-003
  dof 2   d  1.47350849e-003
  dof 3   d -1.20930012e-003
Node           7 (       7):
  dof 1   d  0.00000000e+000
  dof 2   d  7.90019768e-003
  dof 3   d -1.53367475e-003
Node           8 (       8):
  dof 1   d  6.36409779e-004
  dof 2   d  3.12716251e-003
  dof 3   d -2.71155296e-003
Node           9 (       9):
  dof 1   d  0.00000000e+000
  dof 2   d -7.43934900e-003
  dof 3   d  3.55983097e-003




Element output:
---------------
element 1 (       1) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000 -1.7334e-006  3.7402e-007  0.0000e+000 -1.2543e-002  4.0142e-003  0.0000e+000  0.0000e+000  0.0000e+000 -3.2077e-002
              stresses  0.0000e+000  0.0000e+000  0.0000e+000 -1.1667e+003  2.5174e+002  0.0000e+000 -2.1805e+002  4.8344e+000  0.0000e+000  0.0000e+000  0.0000e+000 -2.1591e+002
element 2 (       2) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000  2.7454e-006 -2.2600e-007  0.0000e+000  2.3558e-003  1.2853e-002  0.0000e+000  0.0000e+000  0.0000e+000 -8.8973e-003
              stresses  0.0000e+000  0.0000e+000  0.0000e+000  1.8479e+003 -1.5212e+002  0.0000e+000  1.1946e+002  2.6077e+002  0.0000e+000  0.0000e+000  0.0000e+000 -5.9886e+001
element 3 (       3) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000 -3.3454e-006 -2.7239e-006  0.0000e+000 -1.3949e-002 -1.7826e-002  0.0000e+000  0.0000e+000  0.0000e+000 -8.8314e-003
              stresses  0.0000e+000  0.0000e+000  0.0000e+000 -2.2517e+003 -1.8334e+003  0.0000e+000 -3.7109e+002 -4.2328e+002  0.0000e+000  0.0000e+000  0.0000e+000 -5.9442e+001
element 4 (       4) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000  2.3334e-006 -3.3670e-006  0.0000e+000  1.3949e-002  4.0142e-003  0.0000e+000  0.0000e+000  0.0000e+000 -1.1651e-002
              stresses  0.0000e+000  0.0000e+000  0.0000e+000  1.5706e+003 -2.2663e+003  0.0000e+000  2.9141e+002  1.5767e+002  0.0000e+000  0.0000e+000  0.0000e+000 -7.8422e+001
element 5 (       5) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000  2.3334e-006 -3.3670e-006  0.0000e+000 -1.3949e-002 -4.0142e-003  0.0000e+000  0.0000e+000  0.0000e+000  1.1651e-002
              stresses  0.0000e+000  0.0000e+000  0.0000e+000  1.5706e+003 -2.2663e+003  0.0000e+000 -2.9141e+002 -1.5767e+002  0.0000e+000  0.0000e+000  0.0000e+000  7.8422e+001
element 6 (       6) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000 -3.3454e-006 -2.7239e-006  0.0000e+000  1.3949e-002  1.7826e-002  0.0000e+000  0.0000e+000  0.0000e+000  8.8314e-003
              stresses  0.0000e+000  0.0000e+000  0.0000e+000 -2.2517e+003 -1.8334e+003  0.0000e+000  3.7109e+002  4.2328e+002  0.0000e+000  0.0000e+000  0.0000e+000  5.9442e+001
element 7 (       7) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000  2.7454e-006 -2.2600e-007  0.0000e+000 -2.3558e-003 -1.2853e-002  0.0000e+000  0.0000e+000  0.0000e+000  8.8973e-003
              stresses  0.0000e+000  0.0000e+000  0.0000e+000  1.8479e+003 -1.5212e+002  0.0000e+000 -1.1946e+002 -2.6077e+002  0.0000e+000  0.0000e+000  0.0000e+000  5.9886e+001
element 8 (       8) :
  GP  1.1  :  strains   0.0000e+000  0.0000e+000  0.0000e+000 -1.7334e-006  3.7402e-007  0.0000e+000  1.2543e-002 -4.0142e-003  0.0000e+000  0.0000e+000  0.0000e+000  3.2077e-002
              stresses  0.0000e+000  0.0000e+000  0.0000e+000 -1.1667e+003  2.5174e+002  0.0000e+000  2.1805e+002 -4.8344e+000  0.0000e+000  0.0000e+000  0.0000e+000  2.1591e+002




    R E A C T I O N S  O U T P U T:
    _______________________________


    Node        1 iDof  1 reaction  5.0000e+002    [bc-id: 1]
    Node        3 iDof  1 reaction -5.0000e+002    [bc-id: 1]
    Node        7 iDof  1 reaction  5.0000e+002    [bc-id: 1]
    Node        9 iDof  1 reaction -5.0000e+002    [bc-id: 1]

User time consumed by solution step 1: 0.000 [s]


Finishing analysis on: Fri Jan 11 23:00:53 2013

Real time consumed: 000h:00m:00s
User time consumed: 000h:00m:00s

Re: Adding new elements

Dear svemir,
When you apply loads directly on nodes, they apply to the dofs directly, so "Components 3 0.0 0.0 -1000.0" is interpreted as a -1000 moment load around y-rotation.

I suggest using the vtkxml export module for visualization (which in this case is a great help to find bugs)

Here is a corrected input file;

cct.out
Development Test of CCQ elements -> force in center of plate (plate 1x1 meter)
LinearStatic nsteps 1 nmodules 1
vtkxml tstep_all domain_all primvars 1 1 stype 1
domain 2dMindlinPlate
OutputManager tstep_all dofman_all element_all
ndofman 9 nelem 8 ncrosssect  1 nmat 1 nbc 2 nic 0 nltf 1
node 1 coords 3  -0.5  -0.5   0.0  bc 3 1 0 0
node 2 coords 3   0.0  -0.5   0.0  bc 3 0 0 0
node 3 coords 3   0.5  -0.5   0.0  bc 3 1 0 0
node 4 coords 3  -0.5   0.0   0.0  bc 3 0 0 0
node 5 coords 3   0.0   0.0   0.0  bc 3 0 0 0  load 1 2   
node 6 coords 3   0.5   0.0   0.0  bc 3 0 0 0
node 7 coords 3  -0.5   0.5   0.0  bc 3 1 0 0
node 8 coords 3   0.0   0.5   0.0  bc 3 0 0 0
node 9 coords 3   0.5   0.5   0.0  bc 3 1 0 0
CCTPlate 1 nodes 3 1 2 5 mat 1 crossSect 1
CCTPlate 2 nodes 3 2 3 6 mat 1 crossSect 1
CCTPlate 3 nodes 3 1 5 4 mat 1 crossSect 1
CCTPlate 4 nodes 3 2 6 5 mat 1 crossSect 1
CCTPlate 5 nodes 3 4 5 8 mat 1 crossSect 1
CCTPlate 6 nodes 3 5 6 9 mat 1 crossSect 1
CCTPlate 7 nodes 3 4 8 7 mat 1 crossSect 1
CCTPlate 8 nodes 3 5 9 8 mat 1 crossSect 1
SimpleCS 1 thick 0.01
IsoLE 1 d 0. E 2.1E11 n 0.3 tAlpha 0.000012
BoundaryCondition  1 loadTimeFunction 1 prescribedvalue 0
NodalLoad 2 loadTimeFunction 1 Components 3 -1000.0 0.0 0.0
ConstantFunction 1 f(t) 1.0

Re: Adding new elements

Mikael,
thank you for your help about loading. Now it is working fine. My ccq element is also working fine for this test (only linear static is implemented). Now next step is ccq3d which I implement based using same logic as cct3d. Unfortunatelly I do not know how to test. Problem is domain keyword. Should I use 3DShell domain and then for bc put 6 Ux Uy Uz Rx Ry Rz?
Thanks in advance, Svemir.

Re: Adding new elements

The domain keyword only determines the default unknowns in each node. With arbitrary rotation in 3D then you can't just include the "plate" dofs so that is why there only exists 3dShell (i.e. 3dPlate wouldn't make any sense).

Note that this will include all 3D dofs, which are unlikely to be controlled by plates only, unless your plates are completely interconnected.
You most likely want to go for a shell element, but even in those cases, they still only control 5 dofs, so you need to either;
1. Mess about with local coordinate systems in each node, and prescribe the last, uncontrolled dof (often called the "drilling dof", as it's the rotation normal to the surface)
2. Add a ficticous so "drilling-stiffness".

I'm currently finishing a 2d quadratic mindlin plate, to which I plan to add a drilling stiffness. Then its a simple matter to add a membrane and add arbitrary rotations to allow it to work in 3D.

Re: Adding new elements

Mikael,
Please find attached mine quad plate element and give some comment.
Thanks Svemir.

Post's attachments

ccq3d.zip 16.55 kb, 4 downloads since 2013-01-13 

You don't have the permssions to download the attachments of this post.

Re: Adding new elements

svemir,
First, I'd like to point out that "CCT" stands for constant curvature triangle. For a bilinear interpolation, you just won't have constant curvature, so the very name of the element is misleading. Otherwise, this seems to be exactly the element I have implemented as well.

There are some functions in your code which I do not believe that you have fully checked. Such that the N matrix is 2*8, instead of 3*12.
I suspect you copied these from the plane stress 2d quadrilateral element, which won't work at all.
Most of these functions you have implemented you won't ever need. For example, the NLB matrix is just for large deformations.

I just finished the element I'm working on, which is the bilinear mindlin element in 2D, essentially what you have in CCQ
As far as I can tell, I got the quadratic mindlin element to work now. However, it shows extreme shear locking in the equivalent example above.
I'm trying out reduced integration for the transverse shear strains, but then you get spurious zero stiffness modes. There seems to be methods to deal with this, but it looks like more effort than it is worth.
It's OK for thicker plates; Length/thickness ~ 10, but it is really really bad for Length/thickness ~100 (like in your example).
The only option I see is to go higher order (at least as far as Mindlin theory is concerned).

Re: Adding new elements

Here is an example for the Quad1Mindlin element;

quad1mindlin.out
Testing
NonLinearStatic nsteps 1 nmodules 1 controlmode 1 rtolv 1.e-8 rtold 1e8
vtkxml tstep_all domain_all primvars 1 1 stype 1
domain 2dMindlinPlate
OutputManager tstep_all dofman_all element_all
ndofman 9 nelem 4 ncrosssect  1 nmat 1 nbc 3 nic 0 nltf 1
node 1 coords 3  -0.5  -0.5   0.0  bc 3 1 1 1
node 2 coords 3   0.0  -0.5   0.0  bc 3 0 0 0
node 3 coords 3   0.5  -0.5   0.0  bc 3 0 0 0
node 4 coords 3  -0.5   0.0   0.0  bc 3 1 1 1
node 5 coords 3   0.0   0.0   0.0  bc 3 0 0 0
node 6 coords 3   0.5   0.0   0.0  bc 3 0 0 0
node 7 coords 3  -0.5   0.5   0.0  bc 3 1 1 1
node 8 coords 3   0.0   0.5   0.0  bc 3 0 0 0
node 9 coords 3   0.5   0.5   0.0  bc 3 0 0 0
Quad1Mindlin 1 nodes 4 1 2 5 4 mat 1 crossSect 1
Quad1Mindlin 2 nodes 4 2 3 6 5 mat 1 crossSect 1 boundaryLoads 2 3 2
Quad1Mindlin 3 nodes 4 4 5 8 7 mat 1 crossSect 1
Quad1Mindlin 4 nodes 4 5 6 9 8 mat 1 crossSect 1 boundaryLoads 2 3 2
SimpleCS 1 thick 0.01
IsoLE 1 d 0. E 2.1E11 n 0.3 tAlpha 0.000012
BoundaryCondition  1 loadTimeFunction 1 prescribedvalue 0
NodalLoad 2 loadTimeFunction 1 Components 3 -1.0e6 0.0 0.0
ConstantEdgeLoad 3 loadTimeFunction 1 ndofs 3 loadtype 3 Components 3 -1.0e6 0.0 0.0
ConstantFunction 1 f(t) 1.0

Re: Adding new elements

Mikael,
are you using shear coefficient in material matrix Ds=E/(2*k*(1+ni))*[[1 0][0 1]] (k=1.2) for our problems.
Where I can download your elements and add it to my project. I am using Visual Studio for development. Thanks, Svemir.

Re: Adding new elements

You can check out the development code using GIT from http://www.oofem.org/git/oofem.git  (msysgit, tortoisegit, not sure which is the best tool for windows)
The latest version uses CMake for build scripts, which directly supports Visual Studio.

Regarding the tangent matrix. The elements (both Quad1Mindlin and CCT) lets the material routine decide the tangent.
See LinearElasticMaterial::give2dPlateStiffMtrx for the actual implementation.

Re: Adding new elements

An update, I'm working part time on a flat quadrilateral shell for 3d analysis. I have to work out integration of the ficticous drilling dof stiffness, the sixth unconstrained dof, and then test the element.

From what I've read, the best element would use a serendipity interpolation for the rotations, and full quadratic for displacements. OOFEM currently only have implementations for serendipity interpolation, so I'll add full quadratic (9 node) quadrilateral interpolation to my TODO-list.