Blade stiffened shell improvements

[A-CGray]

This PR contains various improvements to the TACSBladeStiffenedShellConstitutive constitutive class:

1. Adds stiffener column buckling and crippling failure modes
2. Material failure criterion is now computed at bottom as well as the top of the stiffener
3. Adds capability to choose which of the failure modes are used in the overall failure calculation (material failure, global panel buckling, inter-stiffener/local panel buckling, stiffener column buckling, stiffener crippling). At the python level, these can be set through the setFailureModes method
4. Related to 3, the unit tests for this constitutive model now contain separate tests for each failure mode's derivatives as well as tests for the full combined failure envelope
5. Rather than having them hardcoded in the evalFailure and failure sensitivity methods, I have moved the global and local panel buckling calculations to their own methods. This should reduce the amount of code duplication required in @sean-engelstad 's Gaussian Process Buckling Constraints in Blade Stiffened Shell Constitutive Subclass #311 PR

[timryanb]

@sean-engelstad, could you review this PR when you get a chance and see if these new methods can be utilized in your GP panel work (PR #311)?

[sean-engelstad]

Hi @timryanb I should be able to review it later today

[sean-engelstad]

Overall I really like the increased modularity as it will be helpful to reduce code duplication in my subclass. Would like to see some results on applying only the stiffener crippling failure on a flat plate pure axial load case to see what the final stiffener aspect ratios are. Just to double check we implemented it correctly.

[sean-engelstad]

Sorry didn't mean to click approve on the PR earlier. Want to see printout of the failure values for where stiffener column buckling is active. Also want to see on a flat plate (one TACS component case) when only stiffener crippling mode is active, what is the critical stiffener aspect ratio SAR = stiffenerHeight / stiffenerThick.

[sean-engelstad]

@A-CGray if you add a mode and setter for writing out different values to the f5 file or in the DVs (like what I did in my PR), then you can see the failure values. Also, I think you should change the way I did it and have dv1 be argmax(fails) basically so you can see the index of which failure criterion is maximum.

[sean-engelstad]

Here's what I get with the closed-form solutions in TacsGPBladeStiffenedShellConstitutive subclass when I run a copy of your optimization script with stiffener crippling on and no stiffener pitch DV. This should be equivalent to:

python optimize_stiffened_plate.py --includeStiffenerBuckling

Design Vars
{'dvs.dv_struct': array([0.01196416, 0.05817722, 0.00388759])}

Nonlinear constraints
{'CompAndShear.struct_post.eval_funcs.KSFailure': array([1.]),
 'Pressure.struct_post.eval_funcs.KSFailure': array([0.46583251]),
 'Pressure.struct_post.eval_funcs.MaxDispZ': array([0.002])}

Linear constraints
None

Objectives
{'CompAndShear.struct_post.mass_funcs.Mass': array([28.20693225])}

Optimization terminated successfully    (Exit mode 0)
            Current function value: 28.206932245353478
            Iterations: 31
            Function evaluations: 51
            Gradient evaluations: 31
Optimization Complete
-----------------------------------
Optimal sizing:
================================
Stiffener pitch: 125.0 mm
Panel thickness: 11.964156918285498 mm
Skin ply fractions:
[0.13953 0.36047 0.36047 0.13953]
Stiffener height: 58.17722485065222 mm
Stiffener thickness: 3.887590239767674 mm
Stiffener ply fractions:
[0.59402 0.16239 0.16239 0.0812 ]

[A-CGray]

Thanks @sean-engelstad , can you also try running without stiffener crippling? I have a suspicion that's what's causing most of the difference between the two models here

[A-CGray]

For the plate with the default design variables values:

Stiffener pitch: 125.0 mm
Panel thickness: 2.1717 mm
Skin ply fractions:
[0.13953 0.36047 0.36047 0.13953]
Stiffener height: 57.0 mm
Stiffener thickness: 6.4772727272727275 mm
Stiffener ply fractions:
[0.59402 0.16239 0.16239 0.0812 ]

I get these values for each failure mode in the compression+shear load case:

• Skin failure: 0.202345
• Stiffener failure: 0.079251
• Global panel buckling: 1.332269
• Inter-stiffener buckling: 0.428721
• Stiffener column buckling: 0.292454
• Stiffener crippling: 0.113508

These are the values not including the safety factor

[A-CGray]

Update for @timryanb

• Me and @sean-engelstad compared our failure mode values for the default DVs in this plate case
• The biggest difference is actually in the predicted critical loads for global and local panel buckling, my critical axial and shear loads for global buckling are almost an order of magnitude higher than Sean's.
• Looking at the equations we use to compute the critical axial load, mine computes the bending stiffness of the skin+stiffener cross section about it's centroid while it looks like the formula Sean is using may not account for the offset of the skin and stiffener from this centroid, which might explain the much lower critical load.
• Another difference is that Sean's equation is for a plate pinned on all edges, while my model uses the Euler column buckling equation, so there is no stiffening from the sides of the panel.
• We tried to figure out whose model was closer to reality by comparing our critical load predictions to TACS eigenvalue buckling analyses of plates with explicitly modelled stiffeners.
• For all the plates we could successfully run eigenvalue analyses on, our models were actually relatively close (e.g within 10%), however these were all with extremely low-aspect ratio (AR=1) stiffeners
• However, we couldn't successfully run eigenvalue analyses on plates with higher stiffener aspect ratios or plates where the stiffeners were a larger fraction of the overall stiffness. For example, we couldn't successfully run an eigenvalue analysis on the plate above where our models differ by an order of magnitude

[timryanb]

Thanks for looking into this @A-CGray and @sean-engelstad. An order of magnitude difference in strength is a bit concerning.

How are you guys modeling the stiffeners in the eigenvalue analysis? If the goal is to model high aspect ratio stiffeners then the better approach would be to model them as 1D beams, rather than shell elements. I don't think that TACS beams currently have buckling analysis support, so we may have to resort to doing the analysis in NASTRAN?

Since this PR will make a breaking change to how this class calculates failure going forward, I really want to make sure we get things right.

[sean-engelstad]

@timryanb I think I know how to fix this, there is an inconsistency in stiffener and panel D11 centroid in my code. I will also try and put together a higher quality finite element model of the stiffened panel (validated against the literature this week) that we can compare against. @A-CGray as far as I can tell your euler buckling prediction looks correct, and has a similar form as what is used in NASA SP-8007 for cylinder smeared stiffener models

[timryanb]

Sounds good. Thanks @sean-engelstad

[sean-engelstad]

My FEA models of stiffened panels in my ml_buckling repo are working better for stiffener aspect ratios of around 20 now. There was an inconsistency in not using the overall centroid before. And that threw off my shift-and-invert eigenvalue guesses leading me to not be able to run higher stiffener ARs. My closed-form solution is matching better against the FEA models now for larger stiffeners (with 2% error at high ARs). Need to update my TACS code to compute panel Dij about overall centroid and differentiate that step. Then I can post an updated benchmark design of the stiffened panel Ali made

This plot shows the comparison of my closed-form (CF) to the FEA models for gamma = 11.25 (fairly large stiffener with stiffener AR = 20). At high ARs, the stiffener is very long and the buckling loads match well.

At lower ARs, I found that the presence of the stiffener distorts the mode shape away from the CF sine-sine shape and smooths out the mode-switching behavior (causing a difference with the CF although CF is still conservative).

[sean-engelstad]

Comparison of the closed-form solution (TACSGP const. class) and an FEA model

Updated stiffened panel verification case that I put together in my repo, . Here's a panel design and its inputs

Stiffened panel geometry object:
        a = 6.25773160690871
        b = 1.0
        h = 0.01
        num_stiff = 1
        stiffener pitch = 0.5
        h_w = 0.07663302688934862 (height of stiffener)
        t_w = 0.003831651344467431 (thick of stiffener)
        AR = 6.25773160690871
        SR = 100.0
        stiff_AR = 20.0
Plate material:
Composite material object
        E11 = 1.38e11
        nu12 = 0.3
        _E22 = 1.38e11
        _G12 = E11 / 2 / (1+nu12)
        _G23 = None
        _G13 = None
        _ply_angles = [0]
        _ply_fractions = [1.0]
        material_name = None
        symmetric = True
        ref_axis = [1 0 0]
Stiffener material:  same as panel
Advanced parameters in the stiffened panel analysis obj
        rho_0 = 5.9999999999306555
        xi plate = 0.9193240803560768
        xi stiff = 0.9999999999999999
        gamma = 10.000000000000062
        delta = 0.058726208102236284
        zeta plate = 0.00028571428571428574
        zeta stiff = 0.007142857142857143
        stiffAR = 20.0
Mesh + Case Settings
        exx = 9.28592e-05
        exy = 0.00000e+00
        eyy = 0.00000e+00
        num nodes = 2803
In plane loads
        Nxx=128145.6574428344
        Nxy=0.0

Axial

Predicted by TACS GP closed-form: N11,cr = 1.141e6 N/m
Predicted by FEA: N11,cr = 1.0058e6 N/m
Rel Error: 3.58%

Shear

Predicted by TACS GP closed-form: N12,cr = 1.009e6 N/m
Predicted by FEA: N12,cr = 1.018e6 N/m
Rel Error: 0.9%

Updated verification plot from a parametric study

---

Here's also a plot I made to compare CF and FEA in my code for a bunch of models with a single stiffener in the middle. The closed-form solutions are most accurate at unstiffened case or high aspect ratios (although shear is only accurate at high and low ARs). For lower ARs, the mode distortion from the stiffener end constraints results in the global mode distorting to a semi global-local mode and modal assurance criterion doesn't pick it up. At intermediate ARs, the mode distortion of the stiffener results in somewhat higher buckling loads (extra elastic support to the panel than the CF predicts).

Overall, since the buckling loads of the closed-form solution match well at high ARs and are conservative for other ARs - I think these are good results. The Timoshenko : Theory of Elasticity book makes a similar comment on a panel CF of this sort => that they should match well at high ARs and be conservative otherwise.

[timryanb]

Thanks for doing all of this verification work @sean-engelstad. Your results seem to be in good agreement with theory and FEA. Did any of these changes bring your and @A-CGray's results closer together in agreement or are we still waiting on that result?

[sean-engelstad]

@timryanb, I made a new comparison case here. So @A-CGray could you change the plate to this design and post your comparison results? I could also do the optimization with the old stiffened panel => but need to update my derivatives in the const class first.

[A-CGray]

@sean-engelstad , sorry to be a pain but is it possible to run a buckling analysis on a panel with the same cross section but that's wider than it is long? Or at least one where the edges of the plate parallel to the stiffener aren't constrained? I'm pretty certain my global buckling prediction will not be at all meaningful for this case because it assumes that the panel is infinitely wide, whereas in this case the buckling mode is clearly dictated by the panel width.

[sean-engelstad]

@A-CGray, @timryanb here's an updated benchmark using a lower aspect ratio panel using the TACSGPBladeStiffenedShellConstitutive class compared to my TACS FEA analysis conducted in the sean-engelstad/ml_buckling repo I maintain. I also included a plot from my paper on what happens at low aspect ratios for these panels - namely global-local mode mixing.

1. Global-Local Mode Mixing Study

This figure that I've added on stiffened panel verification in my paper shows that once rho_0 the affine aspect ratio is low enough and gamma the EI ratio of stiffener to panel is large enough, global to local mode mixing begins. The greater the number of $N_s$ stiffeners there are, the more you can delay the onset of local modes and preserve a mostly global mode shape. The closed-form is significantly conservative in this regime for metals, but not as conservative for highly orthotropic laminates (if the ply fraction is mostly a 0 deg ply).

2. Low Aspect Ratio Verification Case

I've selected a low aspect ratio panel with enough stiffeners to ensure the mode shapes remain fairly global although there is still some modal distortion. Also I used a stiffener aspect ratio of 5 here, I could go higher, but I found that I sometimes need lower stiffener aspect ratios with more stiffeners otherwise crippling dominates a lot of the solved eigenmodes. This is a simply supported case BC on all sides.

Stiffened panel geometry object:
        a = 0.20277833751554175
        b = 1.0
        h = 0.01
        num_stiff = 8
        stiffener pitch = 0.1111111111111111
        w_b = 0.0
        t_b = 0.0
        h_w = 0.01797546789982607
        t_w = 0.003595093579965214
        AR = 0.20277833751554175
        SR = 100.0
        stiff_AR = 5.0
Plate material:
Composite material object
        E11 = 138000000000.0
        nu12 = 0.3
        _E22 = 138000000000.0
        _G12 = 53076923076.92307
        _G23 = None
        _G13 = None
        _ply_angles = [0.0]
        _ply_fractions = [1.0]
        material_name = None
        symmetric = True
        ref_axis = [1 0 0]
Stiffener material:
Composite material object
        E11 = 138000000000.0
        nu12 = 0.3
        _E22 = 138000000000.0
        _G12 = 53076923076.92307
        _G23 = None
        _G13 = None
        _ply_angles = [0.0]
        _ply_fractions = [1.0]
        material_name = None
        symmetric = True
        ref_axis = [1 0 0]

Advanced parameters in the stiffened panel analysis obj
        rho_0 = 0.2
        xi plate = 0.9727850217307755
        xi stiff = 1.0
        gamma = 1.0000000000000016
        delta = 0.058161140319181945
        zeta plate = 0.00028571428571428574
        zeta stiff = 0.11428571428571428
        stiffAR = 5.0

2.1. Axial Loading:

Here the FEA eigenvalue is quite a bit higher than the closed-form but they are still the same order of magnitude. I think it makes sense the FEA eigenvalue is higher since the stiffeners seem to have influenced the mode shape and appear to be adding some extra restraint to the panel.

Mesh + Case Settings
        exx = 8.78028e-05
        exy = 0.00000e+00
        eyy = 0.00000e+00
        num nodes = 3057
In plane loads
        Nxx=121167.88013431858
        Nxy=0.0
MAC identified global mode 0 with shape (m,n)=(1,1) and max_similarity=0.9981
Eigenvalues in ml_buckling code:
        my CF min eigenvalue= 51.985570043461585
        FEA min eigenvalue= 88.47617896320556

Namely the predicted closed-form vs FEA critical loads by TACS are:

TACS GP closed-form ,  N11,cr = 6,057,629 N/m
TACS FEA prediction,     N11,cr = 10,719,644 N/m
relative error = 43.49%

2.2. Shear Loading

Again the closed-form is somewhat conservative to the FEA here, but still in the same ballpark.

Mesh + Case Settings
        exx = 0.00000e+00
        exy = 2.38255e-04
        eyy = 0.00000e+00
        num nodes = 3057
In plane loads
        Nxx=0.0
        Nxy=126458.42133227034
No shear mode MAC rn just heuristic - first assumed global mode 2 taken

Mode type predicted as global
        my CF min eigenvalue= 150.36762781854705
        FEA min eigenvalue= 194.89827769163952

TACS GP closed-form ,  N12,cr = 18,523,089  N/m
TACS FEA prediction,     N11,cr = 24,532,852  N/m
relative error = 24.4%

The stiffened shear mode for this case is then:

Now if you're wondering whether this is truly a global mode, the unstiffened shear mode for a panel of the same dimensions is the following. I believe that since higher gamma affects the aspect ratio by a 4thRoot(1+gamma) factor that an extra shear half-wave appeared in the above stiffened panel.

[sean-engelstad]

Oops I accidentally closed this

[sean-engelstad]

@timryanb @A-CGray
Update: I discovered the reason why the FEA buckling loads are higher than my closed-form solution by up to 2x or 3x in some cases. As the panel is stiffened, the mode shape goes from SSSS (all simply supported) to approximately SCSC with the x_1 edges becoming clamped. I believe that for low aspect ratio panels as well, the panel width is less important and so the fact that it is nearly clamped along the length direction => it now more closely matches the Euler clamped beam which has 4x the eigenvalue of the simply supported beam.

[A-CGray]

@timryanb , just ran @sean-engelstad 's case, here's what I get:

N1Crit = 6.188025e+06, N12Crit = 2.037140e+07

Comparing all 3 methods:

    TACS GP closed-form: N11,cr =  6,057,629 N/m  N12,cr = 18,523,089 N/m
TACSBladeStiffenedShell: N11,cr =  6,188,025 N/m  N12,cr = 20,371,400 N/m
    TACS FEA prediction: N11,cr = 10,719,644 N/m  N12,cr = 24,532,852 N/m