Quad Issues

[JWock82]

Quads still seem to have some issues when they are not aligned with the global coordinate system. Currently working on a fix to the element formulation.

[KASR]

Hello, I also noticed this when I tried using the pynite code in Matlab. I used a different procedure to obtain the transformation matrix. I tried to translate the matlab code to python. You can find the routine below, I'm not certain about the used node convention so i tried to include some ASCII art to show the node numbering and the element axis system.

On a side note, I do K_global = T' * K_local * T to transform the local matrices to the global ones (referring to an earlier note by you where you changed it to the inverse)

 def T(self):
 '''
 Returns the coordinate transformation matrix for the quad element.

            n *---------------* k
              |               |
              |               |
              |      ye       |
              |      |_ xe    |
              |               |
              |               |
            i *---------------* j
                 shell geometry
 '''

 # find 2 vectors in the plane of the shell
 v1_x = 0.5*( self.j_node.X + self.k_node.X - self.n_node.X - self.i_node.X)
 v1_y = 0.5*( self.j_node.Y + self.k_node.Y - self.n_node.Y - self.i_node.Y)
 v1_z = 0.5*( self.j_node.Z + self.k_node.Z - self.n_node.Z - self.i_node.Z)

 v2_x = 0.5*( self.n_node.X + self.k_node.X - self.j_node.X - self.i_node.X)
 v2_y = 0.5*( self.n_node.Y + self.k_node.Y - self.j_node.Y - self.i_node.Y)
 v2_z = 0.5*( self.n_node.Z + self.k_node.Z - self.j_node.Z - self.i_node.Z)

 # Normalize v1
 v1 = array([v1_x, v1_y, v1_z])
 v1 = v1/norm(v1)

 # Use Gram–Schmidt to find v2
 v2 = array([v2_x, v2_y, v2_z])
 alpha = dot(v2, v1)
 v2 = v2 - (alpha*v1)
 v2 = v2/norm(v2)

 # Use the cross product to find v3
 v3 = cross(v1, v2)

 # Create the direction cosines matrix.
 '''           
               -                 -     
              | v1[0] v1[1] v1[2] |
    dirCos =  | v2[0] v2[1] v2[2] |
              | v3[0] v3[1] v3[2] |
               -                 -  
 '''
 dirCos = array([v1,v2,v3])

 # Build the transformation matrix.
 T = zeros((24, 24))
 T[0:3, 0:3] = dirCos
 T[3:6, 3:6] = dirCos
 T[6:9, 6:9] = dirCos
 T[9:12, 9:12] = dirCos
 T[12:15, 12:15] = dirCos
 T[15:18, 15:18] = dirCos
 T[18:21, 18:21] = dirCos
 T[21:24, 21:24] = dirCos

 # Return the transformation matrix.
 return T

[KASR]

For sake of completeness, i use the same v1 and v2 to find the local nodal coordinates. You can find the code below. Again sorry for my lack of python knowledge. But hopefully you get the idea.

def _local_coords(self):
'''
Calculates or recalculates and stores the local (x, y) coordinates for each node of the
quadrilateral.

            n *---------------* k
              |               |
              |               |
              |      ye       |
              |      |_ xe    |
              |               |
              |               |
            i *---------------* j

                 shell geometry
 '''

 # find 2 vectors in the plane of the shell
 v1_x = 0.5*( self.j_node.X + self.k_node.X - self.n_node.X - self.i_node.X)
 v1_y = 0.5*( self.j_node.Y + self.k_node.Y - self.n_node.Y - self.i_node.Y)
 v1_z = 0.5*( self.j_node.Z + self.k_node.Z - self.n_node.Z - self.i_node.Z)

 v2_x = 0.5*( self.n_node.X + self.k_node.X - self.j_node.X - self.i_node.X)
 v2_y = 0.5*( self.n_node.Y + self.k_node.Y - self.j_node.Y - self.i_node.Y)
 v2_z = 0.5*( self.n_node.Z + self.k_node.Z - self.j_node.Z - self.i_node.Z)

 # normalize v1
 v1 = array([v1_x, v1_y, v1_z])
 v1 = v1/norm(v1)

 # use Gram–Schmidt process to define v2
 v2 = array([v2_x, v2_y, v2_z])
 alpha = dot(v2, v1)
 v2 = v2 - (alpha*v1)
 v2 = v2/norm(v2)

 # Get the global coordinates for each node
 xglob_i = array([self.i_node.X, self.i_node.Y, self.i_node.Z])
 xglob_j = array([self.j_node.X, self.j_node.Y, self.j_node.Z])
 xglob_k = array([self.k_node.X, self.k_node.Y, self.k_node.Z])
 xglob_n = array([self.n_node.X, self.n_node.Y, self.n_node.Z])

# Calculate the local (x, y) coordinates for each node
self.x1 = dot(xglob_i, v1)
self.x2 = dot(xglob_j, v1)
self.x3 = dot(xglob_k, v1)
self.x4 = dot(xglob_n, v1)

self.y1 = dot(xglob_i, v2)
self.y2 = dot(xglob_j, v2)
self.y3 = dot(xglob_k, v2)
self.y4 = dot(xglob_n, v2)

[JWock82]

Does this work any better for you?

I'm having a hard time following how you came up with v1 and v2. I'm having a hard time conceptualizing what those formulas represent.

Awesome ASCII art by the way.

[JWock82]

FYI: I'm in the process of rebuilding the Quad3D element to use a simpler MITC4 formulation I found. I'm mostly done, but am still getting hung up on something I can't seem to isolate when the elements form a cylinder. That would lead me to believe the transformation matrix is at fault, but the formulation as I have it makes sense to me.

[KASR]

Does this work any better for you?

I'm having a hard time following how you came up with v1 and v2. I'm having a hard time conceptualizing what those formulas represent.

Awesome ASCII art by the way.

to give you a bit of background information. I was looking for a good shell element for topology optimization purposes. I came across your implementation after some extensive google searching and decided to translate it towards Matlab. Your detailed comments made it very understandable.

However, when I did some tests using half a sphere with some loads on, and compared the deformations with NX Nastran (v2021) I noticed that the results were not ok. To give an indication of the changes I made:

• adapted the order of the nodes and the shape functions so that order is as shown below
• adapted the computation of the drilling degrees of freedom (i use the method proposed by Hughes in 1989, if you want i can send you the paper)
• adapted the evaluation of the local coordinates and the transformation matrix

If you want I can e-mail you the resulting Matlab code. In the weekend I will reevaluated the half-sphere model so that I can visualize the results and the comparison with Nastran.

However, it is important to note that I still assume that each shell is flat, even when they are on a cylinder or sphere or curved structure. This means that if the mesh is too coarse an error is introduced. Fixing this is a bit more complicated since the transformation matrix then requires different axis systems (or direction vectors as bathe calls them) for each node. Essentially the transformation matrix now computes the element axis system, and it is then shared with the nodes. However if the shell is curved, the out-of-plane axis at the nodes will be pointing in another direction. Hence an error is introduced.

I have considered to first compute the axis for the elements, and then average these for the nodes in order to get a closer match. (I hope this makes a bit of sens without making some plots). In the weekend i will also make some plots to visualize what i mean.

        4 *---------------* 3
          |               |
          |               |
          |      ye       |
          |      |_ xe    |
          |               |
          |               |
        1 *---------------* 2

[JWock82]

I found the issue. The membrane stiffness was too low, which was requiring the bending stiffness to make up the difference. I had an extra 1/4 factor in the membrane stiffness B matrix calculations. This explains why the results looked good for a flat plate with no in-plane loads, but not good for a cylinder that had hoop tension. The fix is very simple and I'll roll it out soon.

The transformation matrix was not the problem. The Quad3D element's transformation matrix is an identical formulation to the Plate3D element's, which has been performing fine. I went so far as to bring in the Plate3D element's T matrix code and the problem was still there.

[JWock82]

This error is now fixed as of v0.0.52