Quad Stress Clarification & Input Checks

[SoundsSerious]

I have a need to predict failure on mesh elements (which I am fairly rusty on) and have taken a stab at improving the documentation of Quad.member() function. I thought it would be more useful to put this in a PR vs an Issue so as to reduce your workload.

I've studied your code as well as some background on MITC4, and believe your member() function is designed to return the normal & shear stresses which I need for my use case.

I've made two changes here which shouldn't affect results

1. Add Bound Protection to shear,moment,member functions to ensure (r,s) are within -1,1

2. Stress Description to Quad3D.member() function as to what I think is appropriate (but I may be mistaken)

[SoundsSerious]

I'm thinking it might be a good idea to add a von-mises stress calculation to plate and quad elements, perhaps the max stress of the element (i think thats 9 points to consider per quad) to simplify use.

This reference has a calculation for von-mises stress assuming Sigma Z=0: https://docs.bentley.com/LiveContent/web/STAAD.Pro%20Help-v15/en/GUID-777CEFB1-4F20-457A-8594-C3F8288A1AF5.html

VM = 0.707[(SMAX - SMIN)^2 + SMAX^2 + SMIN^2]^(1/2) where SMAX / SMIN are the principle stresses.

To get the principle stresses, I suppose my question here would be is calculating Txz and Tyz is as simple as taking Qx and Qy and dividing by the plate thickness, and then solving as usual except with SigmaZ=0 for principle stresses?

[SoundsSerious]

Took an initial stab at the von mises analysis. Before I commit it I'd like to get your input on if you indeed want to add this, and how you'd like to structure this as part of the Quad / Plate, or perhaps make a class mixin to provide some common interface in light of the reorganization efforts.

#Mesh Utils
def quad_stress_tensor(quad,r=0,v=0,combo='gravity'):
    """determines stresses from a plate or quad"""
    q = quad
    S_xx,S_yy,Txy = q.membrane(r,s,combo)
    Qx,Qy = q.shear(r,s,combo)

    S_xx,S_yy,Txy = float(S_xx),float(S_yy),float(Txy)

    T_zx = float(Qx) / q.t
    T_yz = float(Qy) / q.t

    return np.array([   [S_xx, Txy, T_zx],
                        [Txy, S_yy, T_yz],
                        [T_zx, T_yz, 0.] ])

def node_pt(node):
    return np.array([node.X,node.Y,node.Z])

def node_info(quad,material):
    #make point arrays
    ind = node_pt(quad.i_node)
    jnd = node_pt(quad.j_node)
    nnd = node_pt(quad.n_node)
    mnd = node_pt(quad.m_node)

    #area of triangle = 0.5 (V1 x V2)
    C1 = np.cross(jnd-ind,nnd-ind)
    C2 = np.cross(jnd-mnd,nnd-mnd)

    #info!
    A = float((norm(C1) + norm(C2))*0.5)
    V = A*quad.t
    #mass and cost aren't currently supported features in PyNite 
    #mass= V * material.density

    #output
    out = dict(
            A=A,
            cog = (ind+jnd+mnd+mnd)/4,
            V = V,
            #mass=mass,
            #cost = mass* material.cost_per_kg
            )
    return out

def calculate_von_mises_stress(stress_tensor):
    e_val, e_vec = np.linalg.eigh(stress_tensor) #eigh has improved performance for symmetric matricies
    p3, p2, p1 = np.sort(e_val)
    return 0.707 * ((p1-p2)**2 + (p2-p3)**2 + (p1-p3)**2)**0.5

def calculate_quad_vonmises(quad,combo) -> dict:
    """looks at 9 points over a quad to determine its von mises stresses
    :returns: a dictinary of """
    qvm = {}
    for r,s in itert.product([-1,0,1],[-1,0,1]): 
        S = quad_stress_tensor(q,r,s,combo)
        qvm[(r,s)] = calculate_von_mises_stress(S)
    return qvm

[SoundsSerious]

I have made one more improvement to the mesh code, I found that AnnulusMesh had its generate() commented out. I have uncommented it and it seems to work as intended.

[JWock82]

I think by Quad.member() you mean Quad.membrane(). That method returns the in-plane stresses, often referred to as membrane stresses. Those stresses are the axial stresses for a plate, and the shear stresses (see image below).

From these stresses one can calculate the principal stresses and the von mises stresses.

There should be no need to turn to a commercial program's documentation for von mises stress calculation. Just about every text on mechanics of materials has those formulae published. I'd prefer not to use other programs as a reference point for Pynite. That could make Pynite a target for copyright infringement. I'm trying to develop everything I can from textbooks.

I've thought about adding von mises stresses a few times. I think a simple method built into the quad element object called Quad.von_mises() could do it with only a few lines of code that reference Quad.membrane(). I'm not sure what you're doing with tensors, but you should be able to calculate the von-mises stresses at any given point in the quad using the (r, s) coordinate system to extract membrane stress values and convert them to von-mises values. I think it's just $\sqrt{\sigma_a^2-\sigma_a\sigma_b+\sigma_b^2}$ where $\sigma_a$ and $\sigmab$ are the principal stresses calculated from $\sigma{a,b} = \dfrac{\sigma_x+\sigma_y}{2} \pm \sqrt{\left(\dfrac{\sigma_x-\sigmay}{2}\right)^2+\tau{xy}^2}$. Maybe I'm oversimplifying it?

I've toyed with the idea of adding this feature. Von mises stresses can be useful for an isotropic material like steel where they can be compared directly to an allowable stress. Even then, I think mechanical engineers use this much more than structural engineers do. The American Institute of Steel Construction has a lot to say about plate failure, so we rarely find ourselves using von-mises as a failure criterion directly. Usually our equations for capacity are based on test-results. Von-mises can provide a false sense of security as it does not consider the possibility of buckling modes. Von-mises also ignores post-yield and post-buckling behavior such as tension-field action. When it comes to concrete, the rebar typically runs two orthogonal dimensions that are not aligned with the element's principal axes (not to be confused with the element's local axes), making von-mises stresses irrelevant.

There are a lot of methods in your code to get to one stress result. I'd prefer if you could condense it into a single method that calculates von-mises stresses. That will help keep everything related to von-mises stresses under one roof, making the code easier to follow.

[SoundsSerious]

Sorry for the confusion. You're correct I did mean Quad.membrane() , I'm not a full time structures guru so my terminology could use some improvement.

So first off it sounds like this is really more of a "design philosophy" thing rather than a PR, so maybe its a better idea for me to close this and move to a discussion since no one likes open PR's or worse issues :)

I think you're description of the problem is spot on, I think outside my simple world of isotropic materials it seems like a bigger integration issue to maintain failure models & calcs. Its a completely separate problem.

That being said I think if there was such a thing that would be really excellent! I was thinking a good integration might be in the materials space & failure models. I know sectionproperties has almost verbatim the same materials definition that I combined [here] (https://github.com/Ottermatics/ottermaticslib/blob/master/ottermatics/eng/solid_materials.py) into SolidMaterial . I'm attempting to combine the material properties and failure interface into the same object so the interface stress = f(strain) and stiffness(space) = f(temp,P,...) becomes possible, as well as some kind of failure(stress,cycles...) boolean. That being said thats a completely separate problem for me too, ha!

[JWock82]

I think von-mises makes sense to include. It's a good metric for identifying potential trouble spots in a model.

As far as philosophies for design go, there's really not much personal philosophy involved. In structural engineering (I'm talking civil structures, rather than aerospace or mechanical structures), our designs are highly regulated by building codes, rather than directly by principles of mechanics. We don't get to philosophize much outside of what the codes tell us to do. The codes are based on mechanics, but supplemented greatly by data from experimental tests and research and experience that identify all kinds of failure modes and serviceability criteria.

Steel, aluminum, wood, concrete, and masonry all have different codes governing design with similar but different equations. All of these codes in the US have switched to an LRFD approach, which looks at the ultimate post-yield failure mechanism under factored loads, rather than taking an allowable stress approach under unfactored loads. For example, instead of using the elastic section modulus (S) to calculate bending strength (Myield=Fy*S), we use the plastic section modulus (Mplastic=FyZ). And when we're talking about earthquakes, we are going well beyond yield, to the point of paying close attention to member ductility which indicates how far past yield we can go.

With a different code for each material, and countries using different codes, implementing design in Pynite would be a daunting task. I've chosen to focus primarily on analysis for the time being, at least until I feel I've got all the basics covered.

Along these same lines, I am currently playing around with non-linear materials in Pynite, trying to see how difficult it would be to implement elastoplastic behavior in Pynite models. So far it looks promising - at least for framing elements (other than plates). I would love to get it doing pushover analysis for earthquakes.

[SoundsSerious]

Thats a great synopsis! I have alot to learn about the world of structures. I am having a blast using the library though, and see it as a great tool for optimizing for light weight & strong designs which is really important in the vehicle space. Having von-mises stress would be great for that purpose.

When I was saying more of a philosophy question I was thinking more along the lines of what belongs in PyNite vs what is out of scope. Its easy enough right now to mod a fork for my custom use cases so no worries there, I suppose I am wondering where contributions are needed?

I have also been thinking about non-linear materials alot and I have some background in numerical / solver design so maybe this is something I can help with. Recalling my own matlab beam & plate solver from my college days its all about modifying the stiffness matrix, I see the current convergence strategy is to toggle the active boolean on an member. I wonder if it might be faster to define something that is piecewise tensor based on strain rather than to modify out of the inner loop. Maybe also breaking down the stiffness matrix into a non-linear subproblem.