A-CGray commented 5 months ago

This PR contains various improvements to the TACSBladeStiffenedShellConstitutive constitutive class:

Adds stiffener column buckling and crippling failure modes
Material failure criterion is now computed at bottom as well as the top of the stiffener
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
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
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 #311 PR

timryanb commented 4 months ago

@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 commented 4 months ago

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

sean-engelstad commented 4 months ago

@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.

timryanb commented 4 months ago

@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.

@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, there should be a cleaner way of doing this. Maybe allowing for multiple failure values to be returned in the f5 file (similar to how dvs are handled)? But writing failure information out to the dv field seems confusing

sean-engelstad commented 4 months ago

@timryanb Yeah if there is a way to add extra outputs to the f5 file that say fail0 or panelStress, etc. and then failIndex for the argmax(fails) part that would be better.

timryanb commented 4 months ago

In that case, let's make this part a separate PR/feature. I'll add it to the issues page

timryanb commented 4 months ago

See #320 @sean-engelstad @A-CGray

A-CGray commented 4 months ago

Just to confirm @timryanb @sean-engelstad , the plan is to merge this PR, then make the necessary updates to #311 and merge that?

timryanb commented 4 months ago

That's fine by me @A-CGray @sean-engelstad. I'd still like to see a benchmark comparison between the two classes for the same model.

A-CGray commented 4 months ago

Sure, for reference @sean-engelstad , here are the results with my constitutive model:

Fixed stiffener pitch no stiffener crippling constraint

>>>python optimize_stiffened_plate.py 
.
.
.
Objectives
{'CompAndShear.Mass': array([10.79339379])}

Optimization terminated successfully    (Exit mode 0)
            Current function value: 10.793393789903986
            Iterations: 8
            Function evaluations: 10
            Gradient evaluations: 8
Optimization Complete
-----------------------------------
Optimal sizing:
================================
Stiffener pitch: 125.0 mm
Panel thickness: 2.9308217613181076 mm
Skin ply fractions:
[0.13953 0.36047 0.36047 0.13953]
Stiffener height: 115.00000000000001 mm
Stiffener thickness: 1.6478083873616798 mm
Stiffener ply fractions:
[0.59402 0.16239 0.16239 0.0812 ]

Fixed stiffener pitch with stiffener crippling constraint

>>>python optimize_stiffened_plate.py --includeStiffenerBuckling
.
.
.
Objectives
{'CompAndShear.Mass': array([10.89286276])}

Optimization terminated successfully    (Exit mode 0)
            Current function value: 10.89286276157974
            Iterations: 13
            Function evaluations: 17
            Gradient evaluations: 13
Optimization Complete
-----------------------------------
Optimal sizing:
================================
Stiffener pitch: 125.0 mm
Panel thickness: 2.948897670317985 mm
Skin ply fractions:
[0.13953 0.36047 0.36047 0.13953]
Stiffener height: 114.20492349292782 mm
Stiffener thickness: 1.6794607381336215 mm
Stiffener ply fractions:
[0.59402 0.16239 0.16239 0.0812 ]

With stiffener pitch DV, no stiffener crippling

>>>python optimize_stiffened_plate.py --useStiffPitchDV
.
.
.
Objectives
{'CompAndShear.Mass': array([9.72264797])}

Optimization terminated successfully    (Exit mode 0)
            Current function value: 9.722647972607012
            Iterations: 9
            Function evaluations: 11
            Gradient evaluations: 9
Optimization Complete
-----------------------------------
Optimal sizing:
================================
Stiffener pitch: 159.99999999997792 mm
Panel thickness: 3.701976915981655 mm
Skin ply fractions:
[0.13953 0.36047 0.36047 0.13953]
Stiffener height: 150.00000000000003 mm
Stiffener thickness: 0.8902831226396368 mm
Stiffener ply fractions:
[0.59402 0.16239 0.16239 0.0812 ]

With stiffener pitch DV and stiffener crippling

>>>python optimize_stiffened_plate.py --useStiffPitchDV --includeStiffenerBuckling
.
.
.
Objectives
{'CompAndShear.Mass': array([10.88611761])}

Optimization terminated successfully    (Exit mode 0)
            Current function value: 10.886117606613414
            Iterations: 12
            Function evaluations: 17
            Gradient evaluations: 12
Optimization Complete
-----------------------------------
Optimal sizing:
================================
Stiffener pitch: 124.0118771982418 mm
Panel thickness: 2.9291042508631016 mm
Skin ply fractions:
[0.13953 0.36047 0.36047 0.13953]
Stiffener height: 114.01187897381736 mm
Stiffener thickness: 1.6777439670951486 mm
Stiffener ply fractions:
[0.59402 0.16239 0.16239 0.0812 ]

sean-engelstad commented 2 months ago

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 commented 2 months ago

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 commented 2 months ago

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 commented 2 months ago

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 commented 2 months ago

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 commented 2 months ago

@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 commented 2 months ago

Sounds good. Thanks @sean-engelstad

sean-engelstad commented 2 months ago

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. gamma_11_25 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). distorted_mode

sean-engelstad commented 1 month ago

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, ml_buckling . 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%

Screenshot from 2024-09-16 20-01-55

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 commented 1 month ago

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 commented 1 month ago

@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 commented 1 month ago

@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 commented 1 month ago

@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 commented 1 month ago

Oops I accidentally closed this

sean-engelstad commented 1 month ago

@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. bc-study drawio

A-CGray commented 1 month ago

@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

smdogroup / tacs

Blade stiffened shell improvements #319

Fixed stiffener pitch no stiffener crippling constraint

Fixed stiffener pitch with stiffener crippling constraint

With stiffener pitch DV, no stiffener crippling

With stiffener pitch DV and stiffener crippling

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

Axial

Shear

Updated verification plot from a parametric study

1. Global-Local Mode Mixing Study

2. Low Aspect Ratio Verification Case

2.1. Axial Loading:

2.2. Shear Loading