Possible Skewed Quad Error

[JWock82]

Quad results for skewed plates might be off. I've had trouble transforming local forces to global forces for these elements, and have come to the conclusion the issue might be in my implementation of the MITC4 element. I'm currently working on a DKMQ element to replace the MITC4 element, which has a simpler implementation, provides better results, and will provide a second check against whether the plate formulation itself is responsible for this force transformation issue. More to come on this.

[bjhowie]

Can ask what sort of issues you've been encountering?

I've been running tests on skewed quad meshes and comparing the results to rectangular meshes and the results seem to align quite well.

[bjhowie]

Can I ask what sort of issues you've been encountering?

I've been running tests on skewed quad meshes and comparing the results to rectangular meshes and the results seem to align quite well.

[JWock82]

@bjhowie This is the issue you and I have been discussing. I'm almost certain I've got the local to global plate stress transformation equations correct, and yet when I do it to skewed quads, the final global plate stress results appear to be off by a fair amount. If my transformation equations are right, that means the quad could be the problem. The unit tests I have set up only check rectangular quads, so an error in skewed quads could be lurking. The derivation I've used for the MITC4 is mathematically intensive, and I've found a more straightforward derivation for a DKMQ element, which is an improvement over the MITC4. I've been meaning to implement it for some time. Now seems like a natural time to do it. I'm going to use the DKMQ as a benchmark for my force transformation. If the DKMQ has problems too, I'll know it wasn't the MITC4 element itself that was the problem.

[bjhowie]

I figured as much. Curious that my local axis modification method gives good results for skewed quads while the global transformation method doesn't seem too. They should be equivalent, but the difference may suggest the issue isn't with the element itself? In either case, cross checking with the DKMQ element is a great idea. Good luck!

[bjhowie]

Craig, well done on finally merging the DKMQ branch. We all really appreciate your efforts! I havent had a chance to run it through all my tests yet, but will over the next week or two.

I've been playing around with my old spring-raft footing model with the new elements and am still observing similarly strange behavior as we were seeing with the MITC4 elements. The below two models are exactly the same size and with exactly the same loads, but the model on the right uses a quad meshing algorithm that generates quads that arent necessarily aligned with the global axes. The global MY quad moments are extracted in both cases (though for the model on the left the My and MY results are basically identical, as you would expect). As you can see, the model with un-aligned quads is not generating results as expected.

I had the quad result extraction working very well with the MITC4 elements using user-defined element local axes and would be keen to implement something similar for the DKMQ elements if you would consider it? I like this approach better as it seems to be the method used by most commercial applications, and it offers the ability to extract results about any arbitrary axis instead of the just the plate local or global axes.

As a side note, I notice the DKMQs are showing a greater concentration of moment at the point load locations than the MITC4 elements.

[angelmusonda]

There is a problem with the transformation from local coordinates to the global coordinate system in the current implementation. Since the bending moment is a second-order tensor, the transformation used here should be as below:

       dir_cos = self.T()[:3, :3]

        Mx = float(Mx)
        My = float(My)
        Mxy = float(Mxy)
        M_local = np.array([
            [Mx, Mxy, 0.0],
            [Mxy, My, 0.0],
            [0.0, 0.0, 0.0]
        ])

        M_global_tensor = dir_cos @ M_local @ dir_cos.T

        Mx_global = M_global_tensor[0, 0]
        My_global = M_global_tensor[1, 1]
        Mxy_global = M_global_tensor[0, 1]

This way, both approaches should give the same results. However, for flexibility, I would prefer @bjhowie's approach.

[JWock82]

I will make this my next priority. I think this may be the solution I’ve been looking for. I must admit I have struggled understanding the concept of tensors. Perhaps that’s where I need to study up.

Not to knock @bjhowie, who has contributed some great ideas to this project and has been instrumental in its progress, but I can’t explain why he was able to get his method to work. I think he may have stumbled upon a situation where it just happened to work out, perhaps because he’s working in a global plane that allows it to work, or some other reason I don’t understand? He redefined the local x-axis for the plate before deriving it. In all the papers for the plate formulation, the local x-axis runs between the first two nodes. All the formulas are based on that local axis assumption. For a local axis transformation to be valid, it would have to come after the plate is derived, rather than before. In other words, we should be transforming the forces, rather than the element itself.

[bjhowie]

I will admit that I came up with my approach on the assumption that it didn't actually matter how you defined the local X axis, provided you were consistent with your transformations. I suppose thats not necessarily guaranteed.

I have only rigorously tested it on plates aligned with global XY plane. It's possible there are issues when considering other plate alignments, though I project the local axis vector onto the plane of the plate so there theoretically should be no difference unless the plate normal and the local axis Vector happen to match.

Still my preferred approach due to the flexibility it provides.

[JWock82]

@bjhowie maybe it does work, but if so I don’t understand why, so I’ve been hesitant to “canonize” that approach for all the other users. Sometimes mathematics surprises me.

[JWock82]

@angelmusonda this seems to have corrected the issue for bending moments. I'm still experiencing the issue for shears. From what I've been able to gather, those are not tensors, and would not need to use a tensor transformation.

[angelmusonda]

@JWock82 Indeed, shear force should be transformed as a first-order tensor (i.e., a vector). Therefore, the original approach was correct. The resulting maps from a skewed mesh may differ from those of regular quadrilaterals because shear results are highly sensitive to mesh quality. I have verified this behaviour even in Robot Structural Analysis. Regarding the other issue, I reviewed Katili’s 1993 paper and observed that he did not explicitly state that the local x-axis must be defined from node 1 to node 2 for his formulation to be valid. It appears that the choice of local x-axis is arbitrary, provided that the local z-axis is normal to the element plane. I compared the bending moment and shear results using both approaches for regular and skewed meshes and found them to be identical as shown below.

The solid slab is at an angle to the Z-AXIS and is loaded with constant pressure load.

APPROACH 1: USING A LOCAL BASIS WHOSE X-AXIS LIES ALONG VECTOR FROM NODE 1 TO NODE 2

Global Bending Moments - Regular Quads

Global Bending Moments - Skewed Quads

Global Shear Force - Regular Quads

Global Shear Force - Skewed Quads

APPROACH 2: USING A USER-SUPPLIED LOCAL AXIS Global/Local Bending Moments - Regular Quads

Global/Local Bending Moments - Skewed Quads

Global/Local Shear Force - Regular Quads

Global/Local Shear Force - Skewed Quads

My Conclusion The results are identical for both approaches; therefore, either can be used. However, additional testing by more users is encouraged to further verify this finding. Shear results are sensitive to mesh distortion, and further research is needed to understand the reasons for this behaviour.

[JWock82]

@angelmusonda thanks for this further verification of what @bjhowie was finding. Having two people come to the same conclusion with separate models is encouraging. Here's a couple of models I've been testing against each other, one skewed and one not skewed:

from Pynite import FEModel3D
from Pynite.Rendering import Renderer
import math
import random
from numpy import isclose
   

# Create the model
quad_model1 = FEModel3D()
quad_model2 = FEModel3D()

# Define geometry
t = 1  # ft
mesh_size = 1  # ft
a = 10  # ft
b = 15  # ft

# Define a material
E = 57000*math.sqrt(4500)*12**2  # psf
G = 0.4*E  # psf
nu = 1/6
rho = 150  # psf
quad_model1.add_material('Concrete', E, G, nu, rho)
quad_model2.add_material('Concrete', E, G, nu, rho)

# Generate the mesh of quads
quad_model1.add_rectangle_mesh('MSH1', mesh_size, a, b, t, 'Concrete', kx_mod=1, ky_mod=1,
                                element_type='Quad')
quad_model2.add_rectangle_mesh('MSH1', mesh_size, a, b, t, 'Concrete', kx_mod=1, ky_mod=1,
                                element_type='Quad')
quad_model1.meshes['MSH1'].generate()
quad_model2.meshes['MSH1'].generate()

# Add supports to the sides and base of the wall
for node in quad_model1.nodes.values():
    if node.X == 0 or node.X == a or node.Y == 0:
        quad_model1.def_support(node.name, True, True, True, True, True, True)

# Add supports to the sides and base of the wall
for node in quad_model2.nodes.values():
    if node.X == 0 or node.X == a or node.Y == 0:
        quad_model2.def_support(node.name, True, True, True, True, True, True)

# Skew the quads in model 1
for node in quad_model1.nodes.values():

    # Adjust all but the top and bottom rows of nodes
    # if not isclose(node.Y, 0) and not isclose(node.Y, 15):

        # Generate a random number between -0.25 and +0.25
        # rand = random.uniform(-0.25, 0.25)

        # Adjust the node's Y-coordinate
        # node.Y += rand*mesh_size

    # Adjust all but the left, center and right columns of nodes
    if not isclose(node.X, 0) and not isclose(node.X, a/2) and not isclose(node.X, a):

        # Generate a random number between -0.25 and +0.25
        rand = random.uniform(-0.25, 0.25)

        # Adjust the node's X-coordinate
        node.X += rand*mesh_size

# Add uniform surface pressure to the elements
for element in quad_model1.quads.values():
    p = 100
    quad_model1.add_quad_surface_pressure(element.name, p, 'L')

# Add uniform surface pressure to the elements
for element in quad_model2.quads.values():
    p = 100
    quad_model2.add_quad_surface_pressure(element.name, p, 'L')

# Add a load combination to the model
quad_model1.add_load_combo('Live', {'L': 1.0})
quad_model2.add_load_combo('Live', {'L': 1.0})

# Analyze the model
quad_model1.analyze()
quad_model2.analyze()

# rndr1 = Renderer(quad_model1)
# rndr1.annotation_size = 0.25
# rndr1.combo_name = 'Live'
# rndr1.color_map = 'QY'
# rndr1.render_loads = False
# rndr1.render_model()

# rndr2 = Renderer(quad_model2)
# rndr2.annotation_size = 0.25
# rndr2.combo_name = 'Live'
# rndr2.color_map = 'My'
# rndr2.render_loads = False
# rndr2.render_model()

# Get the min/max global forces in the quads
M_min1 = -quad_model1.meshes['MSH1'].max_moment('MY', 'Live')
M_max1 = -quad_model1.meshes['MSH1'].min_moment('MY', 'Live')

M_min2 = quad_model2.meshes['MSH1'].min_moment('My', 'Live')
M_max2 = quad_model2.meshes['MSH1'].max_moment('My', 'Live')

print('MY Moments:')
print('MYmin - Skewed: ', M_min1, ' Unskewed: ', M_min2)
print('MYmax - Skewed: ', M_max1, ' Unskewed: ', M_max2)

M_min1 = quad_model1.meshes['MSH1'].min_moment('MX', 'Live')
M_max1 = quad_model1.meshes['MSH1'].max_moment('MX', 'Live')

M_min2 = quad_model2.meshes['MSH1'].min_moment('Mx', 'Live')
M_max2 = quad_model2.meshes['MSH1'].max_moment('Mx', 'Live')

print('')
print('MX Moments:')
print('MXmin - Skewed: ', M_min1, ' Unskewed: ', M_min2)
print('MXmax - Skewed: ', M_max1, ' Unskewed: ', M_max2)

Q_min1 = quad_model1.meshes['MSH1'].min_shear('QY', 'Live')
Q_max1 = quad_model1.meshes['MSH1'].max_shear('QY', 'Live')

Q_min2 = quad_model2.meshes['MSH1'].min_shear('Qy', 'Live')
Q_max2 = quad_model2.meshes['MSH1'].max_shear('Qy', 'Live')

print('')
print('QY Shears:')
print('QYmin - Skewed: ', Q_min1, ' Unskewed: ', Q_min2)
print('QYmax - Skewed: ', Q_max1, ' Unskewed: ', Q_max2)

Here's the output it creates:

MY Moments: MYmin - Skewed: -556.8890800594435 Unskewed: -556.686825434148 MYmax - Skewed: 165.9885156433813 Unskewed: 164.36916880174425

MX Moments: MXmin - Skewed: -862.2819217178005 Unskewed: -853.4854240702 MXmax - Skewed: 449.0580837212343 Unskewed: 443.96122140472005

QY Shears: QYmin - Skewed: -108.25528885805105 Unskewed: -50.371637743448986 QYmax - Skewed: 357.8403315810722 Unskewed: 360.45804154695253

Everything appears to be in order, except the minimum shears are about double what I was expecting.

[JWock82]

I'm going to close this issue. It seems the plates are working properly now. I do think there is room for improvement with allowing for an arbitrary/custom local x-axis, but that is a different issue altogether than the OP for this issue. I'd be willing to entertain a PR for that if someone wants to add that functionality.