Model Fails Statics Checks

[bjhowie]

Sorry to come again with problems and not solutions, but i'm having some trouble with a proof-of-concept lateral stability model that i have been trying to analyse with PyNite. It is a model of an RC building core and shear wall, with some rigid members connecting the two to represent the rigid floor diaphragm.

The model runs fine, and the stress plots on the quad elements look good, but the model fails statics checks. The global reactions are larger than the sum of the applied forces. The 'check_statics' function built in to the 'analyze' method confirms this.

Any thoughts on what could be causing this?

The code to generate the model is complex so i've pickled the model and attached it in a zip file below:

stability_test.zip

I have subclassed PyNite.FEModel3D so you'll need to use the below code to read the model.

import pickle

infile = ...

class customUnpickler(pickle.Unpickler):
    def find_class(self, module, name):
        if module == "Analysis.AnalysisBase":
            module = "PyNite.FEModel3D"
            name = "FEModel3D"

        return super().find_class(module, name)
    
with open(infile,mode="rb") as f: model=customUnpickler(f).load()

Also, i've noticed that Visualization._PrepContour fails when there are nodes in the model not connected to Quad or Plate elements.
node.contour = sum(node.contour)/len(node.contour) results in a divide by zero error.

[JWock82]

Can you confirm you're using the latest version (0.0.49)? I recently fixed an issue in the statics check related to quadrilateral elements. I want to make sure this is not the same issue.

[bjhowie]

Yes, using 0.0.49

[JWock82]

OK. I'll have a closer look at this.

[JWock82]

Another question. I haven't looked at your model yet. But I did take another look at the code calculating the reactions. I had not updated it to account for spring supports yet. That's a newer feature I recently added. Does your model have spring supports?

[bjhowie]

No, no spring supports. The lowest level of the quad elements have their translations fully fixed, otherwise no other supports.

…

________________________________ From: Craig ***@***.***> Sent: Thursday, 30 September 2021 9:47 AM To: JWock82/PyNite ***@***.***> Cc: bjhowie ***@***.***>; Author ***@***.***> Subject: Re: [JWock82/PyNite] Model Fails Statics Checks (#102) Another question. I haven't looked at your model yet. But I did take another look at the code calculating the reactions. I had not updated it to account for spring supports yet. That's a newer feature I recently added. Does your model have spring supports? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#102 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AK5ZQKMBRCVNENPZITCBRDLUEO6S3ANCNFSM5E2AW5AQ>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[JWock82]

So I've fixed the issue with the divide by zero error. I've also published a new release as v0.0.50 with a few changes to the calculate reactions method. I'm not sure if this fixed it, but I can't check it since the model is pickled. I changed the Node class in version 0.0.50, so I get errors when I run your pickled version. Can you either verify that version 0.0.50 fixes your problem, or send me a new file that's compatible with v0.0.50?

[bjhowie]

Contours seem to be working now. Still having issues with the statics though.

I've re-run and attached the pickled model below. You'll still need to use the custom unpickler above to deal with the subclassed model class.
Let me know if it still doesnt work.

StabilityModel.zip

[JWock82]

I'm trying to narrow the list of possible problems down. Will you confirm a few things for me?

1. Does your model have any quadrilateral surface pressures applied? I didn't see any when I rendered it. It looked like nodal loads only.
2. Are your elements rectangular? If so, can you switch them to rectangular plate elements instead of quads and see if you have the same issue. That would help me determine if it's a problem with the quad element formulation or something else.

[bjhowie]

1. No, node loads only
2. In the model i sent you, they are not perfectly rectangular. I created another model with rectangular plate elements and the issue still seemed to exist, so its likely not related to the quad element formulation

[JWock82]

Update: I’m zeroing in on this bug. I’ve figured out how to replicate it in other models. It seems to occur when the plates aren’t parallel to a global plane. That suggests to me that the transformation matrix for these elements has a bug in some cases. I’ll keep looking at it.

[JWock82]

Alright. I think I've got it figured out. There are two things going on here. First, I had been using the transpose of the transformation matrix instead of the inverse of it to save a little computational overhead. That should work for orthogonal matrices, but when I switch out the transpose with the inverse, the difference between loads and reactions closes a little. I'm not a matrix expert, and maybe this matrix is not orthogonal as I had thought, so I'll switch PyNite to use the inverse since that is mathematically how the algebra works out.

Now for what I believe is the main issue. Plates/quads do not support a drilling degree of freedom. I've built in an artificial weak spring for these elements that prevents model instabilities for drilling loads. The stiffness of this weak spring is somewhat arbitrary, and has no "static" basis. The spring stiffness is calculated as 1/1000 of the smallest bending stiffness in the stiffness matrix. It's a placeholder to prevent a model instability, but it doesn't prevent the user from applying drilling loads to plates.

Your model does not have any drilling loads directly applied, but you do have quads connected to other quads that are not in the same plane and not orthogonal. They appear to be 45 degrees from each other. That means a component of the bending moment in one plate becomes a drilling load on its neighbor. That component becomes "lost" because the plates/quads are not formulated to be loaded that way.

[bjhowie]

That makes some sense. I dont think it's the non-orthogonal quads that are doing it though, the problem seems to be with the connection of the beam elements to the quads. If i remove all of the beam elements the issue is resolved, even with the quads at ~45 degrees. Even if there are no loads applied to or via the beam elements it throws off the base reactions

The takeaway then is that, with the current element formulations, beam elements cannot connect to quad elements at a single point?

Dissapointing outcome. Is there a way that PyNite can at least detect these drilling conditions and throw an exception instead of happily reporting incorrect results?

[JWock82]

You may be interested in the following links:

https://risa.com/risahelp/risa2d/Content/Common_Topics/Modeling%20Tips.htm

https://www.eng-tips.com/viewthread.cfm?qid=472062

It looks like there is a plate formulation out there that includes drilling:

https://onlinelibrary.wiley.com/doi/epdf/10.1002/cnm.1630070102

Maybe at some point I’ll pony up the $50 to purchase the document, and then spend the time to figure out how to implement it.

Based on the engineer-tips discussion above, it looks like some of these formulations are questionable, and diverge when the mesh size is reduced.

In most cases I’m willing to guess these errors are inconsequential. In your case, the floor diaphragm is primarily transferring shear to the shear walls. Drilling applied to the plates would be a secondary load transfer mechanism.

[JWock82]

Also you may try using additional beam elements around the perimeter of the floor attached to the beams and plates. That would resolve the drilling by converting it into a force couple. That’s how my FEA text book suggests handling drilling.

You may also try setting the "J" property of your rigid links to something very small. That would force the load to transfer into the plates using a mechanism other than torsion.

[bjhowie]

Thanks. I found the same eng-tips article and thought it was quite helpful.

I'm not sure that this is the full story though. We shouldnt need to consider any drilling degrees of freedom to make this model solve (unlike the flag pole example in the article). There are many other loads paths that can transmit the load via translational push/pull action, but seemingly any model that is capable of attracting any miniscule amount of drilling load to the quads is thrown off entirely. Maybe i'm missing something.

Reducing torsional stiffness of rigid links makes no difference.

I will keep looking into it.

[JWock82]

I agree there are plenty of stiffer alternate load paths available, which doesn’t make sense why the drilling stiffness has any impact. But when I adjust the drilling stiffness the unresolved load changes. I’ll keep digging into it.

[KASR]

Hello,

Concerning the drilling degree of freedom stiffness, you could try to set the stiffness in the k_m routine by using the minimum eigenvalue of the stress-strain matrix. I commented everything elated to the drilling degrees in the k_b routine and added the following... note that I work with Matlab so I do not know the syntax to compute the eigenvalues using python, but you should get the idea. :)

# Add the term from the unexpanded matrix into the expanded matrix
k_exp[m, n] = k[i, j]

# compute drilling stiffness penalty parameter
#  --> eigenvalues of C_Membrane
k_rz = min( eig( C ) ) 

# Add the drilling degree of freedom's weak spring
k_exp[5, 5] = k_rz
k_exp[11, 11] = k_rz
k_exp[17, 17] = k_rz
k_exp[23, 23] = k_rz

[JWock82]

I'm still looking into this issue. One thing I'm noticing is that the unbalanced forces on all these models seem relatively small. The table below is the statics check from the model you sent me. The unbalanced SumRY and SumRZ forces in the output below are about 3 orders of magnitude smaller than the SumRX forces. The unbalanced SumRMX moments below are one order of magnitude smaller than the SumFX forces, suggesting a moment arm of less than 1/10 of one length unit. I'm not sure what unit system these forces represent so I can't say if the differeces are significant. All the other forces seem pretty well balanced. This may just be a precision error after the math is crunched.

On other models I'm running that have the same issue I'm noticing that the unbalanced forces are minuscule, with zeros out to the 8th decimal place or more. The stress results from the models all compare reasonably to published solutions.

[JWock82]

I looked at this again tonight. The largest of the unbalanced forces is reported by PyNite as 4440 in the table above. The force you've applied to the nodes in your model is 100,000 per node. That means the unbalanced force is only 4.44% of one of the nodal forces you've applied to the model.

I've been trying, but I can't seem to create a model that has a significant force/reaction imbalance. I think more and more that this is a computational precision error related to plates/quads. There is a lot of numerical integration associated with their force and stiffness matrices. I'm not surprised to see a slight imperfection in the results. I'm going to close this issue, but I'm going to keep an eye out for bad statics checks on all my future plate/quad models.

[JWock82]

I've spent the last couple of evenings reading through various formulations for plate drilling. I've also been playing with the drilling stiffness to see what happens. Here's what I've discovered:

1. The unbalanced forces are influenced strongly by the value I use for the drilling stiffness.
2. There are very complex formulations out there for the drilling stiffness terms. These formulations lead to mathematical "correctness". Whether they lead to engineering correctness is debatable.
3. The approach PyNite takes of adding a weak spring is proposed by Bathe in "Finite Element Procedures, 2nd Ed." Example 4.19 (starting on p. 207).
4. This procedure can be used for the DKMQ element which has the same drilling problem as the MITC4 element.
5. Using a drilling stiffness value of near zero leads to irrational displaced shapes and severely imbalanced forces.
6. Using a drilling stiffness value of nearly rigid also leads to irrational displaced shapes and severely imbalanced forces.
7. Bathe recommends using a weak drilling stiffness of 1/1000 of the least diagonal term in the plate's stiffness matrix. PyNite takes this approach, but only considers the diagonal terms related to rotation.
8. PyNite's displaced shapes seem rational with the 1/1000 stiffness term. The force and moment imbalances appear to be insignificant.

I don't have any plans to pursue this issue any further, but I'll reopen it and add a help wanted tag to this issue. I'll also tag this as a "program limitation". If someone wants to take this any further they are welcome to submit a pull-request. In the meantime, I recommend checking the statics for anything that looks significant. I may add a parameter that allows the user to specify the drilling stiffness in the future. That way you can adjust it if the statics check gets out of control. I may also add a feature that looks at the statics check for significant imbalances and reports errors if they get too large.

[JWock82]

@bjhowie try running this model again with the latest version. I’ve corrected an issue I found with quadrilaterals related to their in-plane stiffness.

[bjhowie]

Hey, sorry i never got around to replying to your previous post. Unfortunate outcome, but sounds like it's an inherent limitation of the MITC4 plate formulation.

I ran it again with 0.0.52 but got similar results, in fact the phantom RY reaction increased slightly

The Y and Z reactions are small, but maybe not neglegible depending on the application.

[JWock82]

Thanks bjhowie. I just added more unit testing for plates/quads to compare them against published solutions. Overall I'm finding that the plates are giving very accurate results for out-of-plane, in-plane, and combined forces.

[JWock82]

@bjhowie,
Today I fixed an issue related to statics checks that may resolve this issue. I don't have access to your original file that showcased the issue. Do you still have it handy by chance?

[bjhowie]

I should have it somewhere. I'll have a dig around tonight and see if I can find it.

[bjhowie]

I couldn't find the code that generated the original model, but I used the pickled data attached in the post above to regenerate the model using version 0.0.93.

Your updates definitely seem to have helped the issue. RY and RZ are basically zero. I'm unsure if the RMX reaction is an issue? How is it calculated?

[bjhowie]

OK so i manage to find my old code that generated that model. Doesn't look like your latest changes have helped.

I've used it to generate two more models, one using quad elements and one using plates.
StabilityModels.zip

This is the statics check for the quads model:

And for the plates model:

Fairly similar but the phantom reactions change slightly, which is probably not unexpected.

I've got some members in there that are intended to represent a rigid diaphragm. In the examples above ive given them the material properties of steel, but if I increase the E value to something like 1E15 (effectively infinitely stiff) then the phantom reactions get very large and the RX sum is no longer correct.