Enter ABD matrices directly

Is there a way to directly enter ABD matrix values into a shell element? A colleague of mine recreated this tutorial in ABAQUS to create a .inp file https://www.youtube.com/watch?v=oWWGcRNtOms
This contains the section:

** Section: Section-1
*Shell General Section, elset=Set-1
711541., 23317., 711541., 0., 0., 43060., -1.58519e+06, 0.
0., 5.92951e+06, 0., 1.58519e+06, 0., 194306., 5.92951e+06, 0.
0., 0., 0., 0., 358833., 
*Transverse Shear Stiffness
29442., 29442., 0.
*End Instance

However this doesn’t seem to run with Calculix. So is there a way to enter the values for a general shell section?

Here (chapter 7. Input deck format) you can find a list of keywords supported by CalculiX: https://www.dhondt.de/ccx_2.21.pdf

Those are not available and shells can’t be defined this way in CalculiX. In fact, all shell elements (apart from US3) aren’t true shells - they are expanded to solids.


Hmmm, this seems like it’s going to be harder than I realised! I’m trying model a sparse network of unidirectional tows to form a lattice to reinforce a plate (the arrangement of the lattice to be generated by an optimisation routine). It’s a similar situation to this diagram (from the ABAQUS manual):

So I’ll have lots of varying thicknesses elements, in fact when I have a full lattice almost every element will be a different thickness and number of layers to some of it’s neighbours. At this stage I’m only interested in linear elastic behaviour, so I don’t need to worry about the stresses within layers.
So my plan had been to:

  • work out the number of tows and orientations of the plies in each element (each element will have to be in it’s own set)
  • Perform the CLT calculation to work out the ABD matrix for each element
  • Apply these to shell elements and solve

However it looks like I can’t do this approach in Calculix. I can apply a composite lay up directly to each element in Calculix. I can’t see this working though, what happens for example in the above diagram where ply 3 goes over ply 2? There’s going to be hundreds of areas such as that in my model.

One way would be to model them as solid elements right away. From the schematic you shared, you could project the edges of each ply, create surface elements, and then extrude them in the direction needed. You would not be able to enter the ABD matrix directly. If you are looking to do something like this, it may be better to write your own Matlab or Python script.

I’m pressing on with trying to use the composite solid elements. But I’ve got a problem with orientations, how do I set the orientation of each ply separately? I’ve defined my orientations thus:

1.,0.,0., 0.,1.,0.
3, 0.

1.,0.,0., 0.,1.,0.
3, 90.

1.,0.,0., 0.,1.,0.
3, 45.

If I then do:

** Sections
1,, RED_TRH50
1,, RED_TRH50

Then I get a deformed shape as I’d expect (I can see bend-twist coupling)

However if I do:

** Sections
*SHELL SECTION, ELSET=MaterialSolidElementGeometry2D, COMPOSITE

Then it runs without errors but just looks like it’s not applying any orientation, the deformation is the same whether I set OR0, OR45 or OR90. What am I doing wrong in this case?

You are using shell elements here, not solid elements directly.

Model each of the plies as solid elements in cgx. Then you can use the neigh command in cgx to create TIE constraints between the plies.

Whenever a ply is crossing over a previous ply, the elements have to move up and over, and then down. This would be tedious to do manually, but using a scripting language you can let the computer do the work for you. :slight_smile:

Since you are laying down unidirectional tows, it should not be difficult to automatically generate orientations if you are using second order hex elements (C3D20). Just make sure that the thickness of the elements in the tow is the smallest dimension. Then you can easily determine the local z-axis for the element from the node distances like this Python code:

# A hex element has three basic vectors.
bases = [
    sub(nodes[4], nodes[0]),
    sub(nodes[1], nodes[0]),
    sub(nodes[3], nodes[0]),
bases.sort(key=lambda v: -sum(j * j for j in v))
# The normal is perpendicular to the two longest bases.
normal = normalize(cross(bases[0], bases[1]))

(nodes is a dictionary. The node number is the key and each value is a 3-tuple of float values which represent the coordinates of the node in the global coordinate system as a vector. sub is a function that calculates the difference between two vectors. Likewise, cross calculates the cross-product and normalize normalizes the length of a vector.)

Using the cross-product between the two longest base vectors, we find the z-axis for the element.

If the element is also slightly longer then wide, the local x-axis can also be determined; it is the longest base vector. The local y-axis follows from the cross-product between x and z. I use a similar technique in my calculix-orient script.

1 Like