From 577ffb8970788325ab9c5eda8b13d5c7d76b4bb6 Mon Sep 17 00:00:00 2001 From: Alasdair Gray Date: Mon, 17 Jun 2024 20:32:33 -0400 Subject: [PATCH] Remove changes related to `BladeStiffenedShellConstitutive` --- .../TACSBladeStiffenedShellConstitutive.cpp | 1556 +++++------------ .../TACSBladeStiffenedShellConstitutive.h | 322 +--- tacs/constitutive.pyx | 33 +- tacs/cpp_headers/constitutive.pxd | 49 +- ...test_blade_sitffened_shell_constitutive.py | 111 +- .../test_shell_blade_stiffened_plate_quad.py | 6 +- 6 files changed, 496 insertions(+), 1581 deletions(-) diff --git a/src/constitutive/TACSBladeStiffenedShellConstitutive.cpp b/src/constitutive/TACSBladeStiffenedShellConstitutive.cpp index ae6be3260..bc6c02630 100644 --- a/src/constitutive/TACSBladeStiffenedShellConstitutive.cpp +++ b/src/constitutive/TACSBladeStiffenedShellConstitutive.cpp @@ -21,7 +21,7 @@ bladeFSDT model from previous versions of TACS developed by Graeme Kennedy. #include "TACSMaterialProperties.h" #include "TACSShellConstitutive.h" -const char* const TACSBladeStiffenedShellConstitutive::constName = +const char* TACSBladeStiffenedShellConstitutive::constName = "TACSBladeStiffenedShellConstitutive"; // ============================================================================== @@ -37,9 +37,7 @@ TACSBladeStiffenedShellConstitutive::TACSBladeStiffenedShellConstitutive( TacsScalar _stiffenerHeight, int _stiffenerHeightNum, TacsScalar _stiffenerThick, int _stiffenerThickNum, int _numStiffenerPlies, TacsScalar _stiffenerPlyAngles[], TacsScalar _stiffenerPlyFracs[], - int _stiffenerPlyFracNums[], TacsScalar _flangeFraction, - bool _includeGlobalBuckling, bool _includeLocalBuckling, - bool _includeStiffenerColumnBuckling, bool _includeStiffenerCrippling) { + int _stiffenerPlyFracNums[], TacsScalar _flangeFraction) { this->panelPly = _panelPly; this->panelPly->incref(); @@ -48,12 +46,6 @@ TACSBladeStiffenedShellConstitutive::TACSBladeStiffenedShellConstitutive( this->kcorr = _kcorr; - // --- Failure mode flags --- - this->includeGlobalBuckling = _includeGlobalBuckling; - this->includeLocalBuckling = _includeLocalBuckling; - this->includeStiffenerColumnBuckling = _includeStiffenerColumnBuckling; - this->includeStiffenerCrippling = _includeStiffenerCrippling; - this->numDesignVars = this->numGeneralDV = this->numPanelDV = this->numStiffenerDV = 0; @@ -202,22 +194,22 @@ TACSBladeStiffenedShellConstitutive::TACSBladeStiffenedShellConstitutive( // Arrays for storing failure values, need values for each ply angle at the // top and bottom of the panel and at the tip of the stiffener this->panelPlyFailValues = new TacsScalar[2 * _numPanelPlies]; - this->stiffenerPlyFailValues = new TacsScalar[2 * _numStiffenerPlies]; + this->stiffenerPlyFailValues = new TacsScalar[_numStiffenerPlies]; // Arrays for storing failure strain sensitivities this->panelPlyFailStrainSens = new TacsScalar*[2 * _numPanelPlies]; - this->stiffenerPlyFailStrainSens = new TacsScalar*[2 * _numStiffenerPlies]; + this->stiffenerPlyFailStrainSens = new TacsScalar*[_numStiffenerPlies]; for (int ii = 0; ii < 2 * _numPanelPlies; ii++) { this->panelPlyFailStrainSens[ii] = new TacsScalar[this->NUM_STRESSES]; } - for (int ii = 0; ii < 2 * _numStiffenerPlies; ii++) { + for (int ii = 0; ii < _numStiffenerPlies; ii++) { this->stiffenerPlyFailStrainSens[ii] = new TacsScalar[TACSBeamConstitutive::NUM_STRESSES]; } // Arrays for storing ply failure sensitivities this->panelPlyFailSens = new TacsScalar[2 * this->numPanelPlies]; - this->stiffenerPlyFailSens = new TacsScalar[2 * this->numStiffenerPlies]; + this->stiffenerPlyFailSens = new TacsScalar[this->numPanelPlies]; } // ============================================================================== @@ -767,7 +759,7 @@ void TACSBladeStiffenedShellConstitutive::addStressDVSens( // correct) this->addPanelStressDVSens(scale, strain, psi, &dfdx[this->panelDVStartNum]); - // Transform the psi vector the same way we do for strains + // Transform the psi vector the same way we do for stains TacsScalar stiffenerPsi[TACSBeamConstitutive::NUM_STRESSES]; this->transformStrain(psi, stiffenerPsi); @@ -849,45 +841,45 @@ TacsScalar TACSBladeStiffenedShellConstitutive::evalFailure( // Compute the failure values for each failure mode of the stiffened panel TacsScalar TACSBladeStiffenedShellConstitutive::computeFailureValues( const TacsScalar e[], TacsScalar fails[]) { - // Initialize the failure values to some very large and negative that won't - // contribute to the KS aggregate - for (int ii = 0; ii < this->NUM_FAILURES; ii++) { - fails[ii] = DUMMY_FAIL_VALUE; - } + // --- Panel material failure --- + fails[0] = this->computePanelFailure(e); + + // --- Stiffener material failure --- TacsScalar stiffenerStrain[TACSBeamConstitutive::NUM_STRESSES]; this->transformStrain(e, stiffenerStrain); + fails[1] = this->computeStiffenerFailure(stiffenerStrain); - // --- Panel material failure --- - if (this->includeMaterialFailure) { - fails[0] = this->computePanelFailure(e); + // --- Local panel buckling --- + // Compute panel stiffness matrix and loads + TacsScalar panelStiffness[NUM_TANGENT_STIFFNESS_ENTRIES], + stress[NUM_STRESSES]; + this->computePanelStiffness(panelStiffness); + const TacsScalar *A, *D; + this->extractTangentStiffness(panelStiffness, &A, NULL, &D, NULL, NULL); + this->computePanelStress(e, stress); - // --- Stiffener material failure --- - fails[1] = this->computeStiffenerFailure(stiffenerStrain); - } + // Compute the critical local loads + TacsScalar N1Crit, N12Crit; + TacsScalar D11 = D[0], D12 = D[1], D22 = D[3], D66 = D[5], + L = this->stiffenerPitch; + N1Crit = this->computeCriticalLocalAxialLoad(D11, D22, D12, D66, L); + N12Crit = this->computeCriticalShearLoad(D11, D22, D12 + 2.0 * D66, L); - // --- Local panel buckling --- - if (this->includeLocalBuckling) { - fails[2] = this->evalLocalPanelBuckling(e); - } + // Compute the buckling criteria + fails[2] = this->bucklingEnvelope(-stress[0], N1Crit, stress[2], N12Crit); // --- Global buckling --- - if (this->includeGlobalBuckling) { - fails[3] = this->evalGlobalPanelBuckling(e); - } - - // --- Stiffener column buckling --- - if (this->includeStiffenerColumnBuckling) { - fails[4] = this->evalStiffenerColumnBuckling(stiffenerStrain); - } + TacsScalar D1, D2, D3; + this->computeCriticalGlobalBucklingStiffness(&D1, &D2, &D3); + L = this->panelLength; - // --- Stiffener crippling --- - if (this->includeStiffenerCrippling) { - fails[5] = this->evalStiffenerCrippling(stiffenerStrain); - } + N1Crit = computeCriticalGlobalAxialLoad(D1, L); + N12Crit = this->computeCriticalShearLoad(D1, D2, D3, L); - TacsScalar ksFail = ksAggregation(fails, this->NUM_FAILURES, this->ksWeight); + this->evalStress(0, NULL, NULL, e, stress); + fails[3] = this->bucklingEnvelope(-stress[0], N1Crit, stress[2], N12Crit); - return ksFail; + return ksAggregation(fails, this->NUM_FAILURES, this->ksWeight); } // Evaluate the derivative of the failure criteria w.r.t. the strain @@ -897,82 +889,93 @@ TacsScalar TACSBladeStiffenedShellConstitutive::evalFailureStrainSens( memset(sens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); TacsScalar fails[this->NUM_FAILURES], dKSdf[this->NUM_FAILURES]; + // First compute the sensitivity of the panel failure value + TacsScalar panelFailSens[this->NUM_STRESSES]; + fails[0] = this->evalPanelFailureStrainSens(e, panelFailSens); - // Initialize the failure values to some very large and negative that won't - // contribute to the KS aggregate - for (int ii = 0; ii < this->NUM_FAILURES; ii++) { - fails[ii] = DUMMY_FAIL_VALUE; - } - - TacsScalar stiffenerStrainSens[TACSBeamConstitutive::NUM_STRESSES], - stiffenerMatFailSens[this->NUM_STRESSES]; - memset(stiffenerStrainSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - memset(stiffenerMatFailSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - - TacsScalar stiffenerStrain[TACSBeamConstitutive::NUM_STRESSES]; + // And now for the stiffener failure value, first in terms of the beam + // strains, and then transformed back to shell strains + TacsScalar stiffenerStrain[TACSBeamConstitutive::NUM_STRESSES], + stiffenerStrainSens[TACSBeamConstitutive::NUM_STRESSES], + stiffenerFailSens[this->NUM_STRESSES]; this->transformStrain(e, stiffenerStrain); + fails[1] = this->evalStiffenerFailureStrainSens(stiffenerStrain, + stiffenerStrainSens); + this->transformStrainSens(stiffenerStrainSens, stiffenerFailSens); - // --- Material failure --- - TacsScalar panelFailSens[this->NUM_STRESSES]; - memset(panelFailSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - if (this->includeMaterialFailure) { - // First compute the sensitivity of the panel failure value - fails[0] = this->evalPanelFailureStrainSens(e, panelFailSens); + // --- Local panel buckling --- + // Compute panel stiffness matrix and loads + TacsScalar panelStiffness[NUM_TANGENT_STIFFNESS_ENTRIES], + panelStress[NUM_STRESSES]; + this->computePanelStiffness(panelStiffness); + const TacsScalar *APanel, *DPanel; + this->extractTangentStiffness(panelStiffness, &APanel, NULL, &DPanel, NULL, + NULL); + this->computePanelStress(e, panelStress); - // And now for the stiffener failure value, first in terms of the beam - // strains, and then transformed back to shell strains - fails[1] = this->evalStiffenerFailureStrainSens(stiffenerStrain, - stiffenerStrainSens); - this->transformStrainSens(stiffenerStrainSens, stiffenerMatFailSens); - } + // Compute the critical local loads (no need to compute their sensitivities + // because they're not dependent on the strain)) + TacsScalar N1CritLocal, N12CritLocal; + TacsScalar D11 = DPanel[0], D12 = DPanel[1], D22 = DPanel[3], D66 = DPanel[5], + L = this->stiffenerPitch; + N1CritLocal = this->computeCriticalLocalAxialLoad(D11, D22, D12, D66, L); + N12CritLocal = this->computeCriticalShearLoad(D11, D22, D12 + 2.0 * D66, L); - // --- Local panel buckling --- - TacsScalar localBucklingSens[this->NUM_STRESSES]; - memset(localBucklingSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - if (this->includeLocalBuckling) { - fails[2] = this->evalLocalPanelBucklingStrainSens(e, localBucklingSens); - } + // Compute the buckling criteria and it's sensitivities + TacsScalar N1LocalSens, N12LocalSens, N1CritLocalSens, N12CritLocalSens; + fails[2] = this->bucklingEnvelopeSens( + -panelStress[0], N1CritLocal, panelStress[2], N12CritLocal, &N1LocalSens, + &N1CritLocalSens, &N12LocalSens, &N12CritLocalSens); // --- Global buckling --- - TacsScalar globalBucklingSens[this->NUM_STRESSES]; - memset(globalBucklingSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - if (this->includeGlobalBuckling) { - fails[3] = this->evalGlobalPanelBucklingStrainSens(e, globalBucklingSens); - } - - // --- Stiffener column buckling --- - TacsScalar stiffenerBucklingSens[this->NUM_STRESSES]; - memset(stiffenerBucklingSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - if (this->includeStiffenerColumnBuckling) { - fails[4] = evalStiffenerColumnBucklingStrainSens(stiffenerStrain, - stiffenerStrainSens); - this->transformStrainSens(stiffenerStrainSens, stiffenerBucklingSens); - } + TacsScalar stiffness[NUM_TANGENT_STIFFNESS_ENTRIES], stress[NUM_STRESSES]; + this->computeStiffness(stiffness); + const TacsScalar *A, *B, *D, *As; + TacsScalar drill; + this->extractTangentStiffness(stiffness, &A, &B, &D, &As, &drill); + this->computeStress(A, B, D, As, drill, e, stress); + TacsScalar N1GlobalSens, N1CritGlobalSens, N12GlobalSens, N12CritGlobalSens; + TacsScalar D1, D2, D3; + this->computeCriticalGlobalBucklingStiffness(&D1, &D2, &D3); + L = this->panelLength; + TacsScalar N1CritGlobal = computeCriticalGlobalAxialLoad(D1, L); + TacsScalar N12CritGlobal = this->computeCriticalShearLoad(D1, D2, D3, L); - // --- Stiffener crippling --- - TacsScalar stiffenerCripplingSens[this->NUM_STRESSES]; - memset(stiffenerCripplingSens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); - if (this->includeStiffenerCrippling) { - fails[5] = - evalStiffenerCripplingStrainSens(stiffenerStrain, stiffenerStrainSens); - this->transformStrainSens(stiffenerStrainSens, stiffenerCripplingSens); - } + fails[3] = this->bucklingEnvelopeSens( + -stress[0], N1CritGlobal, stress[2], N12CritGlobal, &N1GlobalSens, + &N1CritGlobalSens, &N12GlobalSens, &N12CritGlobalSens); // Compute the sensitivity of the aggregate failure value to the individual // failure mode values TacsScalar fail = ksAggregationSens(fails, this->NUM_FAILURES, this->ksWeight, dKSdf); - // Sum up the strain sens for each failure mode + // Compute the total sensitivity due to the panel and stiffener material + // failure memset(sens, 0, this->NUM_STRESSES * sizeof(TacsScalar)); for (int ii = 0; ii < this->NUM_STRESSES; ii++) { - sens[ii] = - dKSdf[0] * panelFailSens[ii] + dKSdf[1] * stiffenerMatFailSens[ii] + - dKSdf[2] * localBucklingSens[ii] + dKSdf[3] * globalBucklingSens[ii] + - dKSdf[4] * stiffenerBucklingSens[ii] + - dKSdf[5] * stiffenerCripplingSens[ii]; + sens[ii] = dKSdf[0] * panelFailSens[ii] + dKSdf[1] * stiffenerFailSens[ii]; } + // Now add the sensitivity of the buckling criteria + + // local buckling + N1LocalSens *= dKSdf[2]; + N12LocalSens *= dKSdf[2]; + sens[0] += N1LocalSens * -APanel[0] + N12LocalSens * APanel[2]; + sens[1] += N1LocalSens * -APanel[1] + N12LocalSens * APanel[4]; + sens[2] += N1LocalSens * -APanel[2] + N12LocalSens * APanel[5]; + + // Global buckling + N1GlobalSens *= dKSdf[3]; + N12GlobalSens *= dKSdf[3]; + sens[0] += N1GlobalSens * -A[0] + N12GlobalSens * A[2]; + sens[1] += N1GlobalSens * -A[1] + N12GlobalSens * A[4]; + sens[2] += N1GlobalSens * -A[2] + N12GlobalSens * A[5]; + sens[3] += N1GlobalSens * -B[0] + N12GlobalSens * B[2]; + sens[4] += N1GlobalSens * -B[1] + N12GlobalSens * B[4]; + sens[5] += N1GlobalSens * -B[2] + N12GlobalSens * B[5]; + return fail; } @@ -985,72 +988,215 @@ void TACSBladeStiffenedShellConstitutive::addFailureDVSens( // stuff. It should be rewritten to use only forward or only backward // differentiation in future. const int numDV = this->numDesignVars; - const TacsScalar t = this->panelThick; // Compute the failure values and then compute the // sensitivity of the aggregate failure value w.r.t. them TacsScalar fails[this->NUM_FAILURES], dKSdf[this->NUM_FAILURES]; - this->computeFailureValues(strain, fails); + TacsScalar fail = this->computeFailureValues(strain, fails); ksAggregationSens(fails, this->NUM_FAILURES, this->ksWeight, dKSdf); + // Sensitivity of the panel failure value to it's DVs + this->addPanelFailureDVSens(strain, scale * dKSdf[0], + &dfdx[this->panelDVStartNum]); + + // Add the direct sensitivity of the stiffener failure value w.r.t DVs + // Sensitivity of the panel failure value to it's DVs TacsScalar stiffenerStrain[TACSBeamConstitutive::NUM_STRESSES]; this->transformStrain(strain, stiffenerStrain); + this->addStiffenerFailureDVSens(stiffenerStrain, scale * dKSdf[1], + &dfdx[this->stiffenerDVStartNum]); - // Sensitivity of the panel failure value to it's DVs - if (this->includeMaterialFailure) { - this->addPanelFailureDVSens(strain, scale * dKSdf[0], - &dfdx[this->panelDVStartNum]); + // Add the sensitivity of the stiffener failure value w.r.t. the DVs + // due to the dependence of the stiffener strains on the DVs + TacsScalar stiffenerFailStrainSens[TACSBeamConstitutive::NUM_STRESSES]; + this->evalStiffenerFailureStrainSens(stiffenerStrain, + stiffenerFailStrainSens); + this->addStrainTransformProductDVsens(stiffenerFailStrainSens, strain, + scale * dKSdf[1], dfdx); + + // --- Local panel buckling --- + // Compute panel stiffness matrix and loads + TacsScalar panelStiffness[NUM_TANGENT_STIFFNESS_ENTRIES], + panelStress[NUM_STRESSES]; + this->computePanelStiffness(panelStiffness); + const TacsScalar *A, *D; + this->extractTangentStiffness(panelStiffness, &A, NULL, &D, NULL, NULL); + this->computePanelStress(strain, panelStress); - // Add the direct sensitivity of the stiffener failure value w.r.t DVs - // Sensitivity of the panel failure value to it's DVs - this->addStiffenerFailureDVSens(stiffenerStrain, scale * dKSdf[1], - &dfdx[this->stiffenerDVStartNum]); + // Compute the critical local loads and their sensitivities + TacsScalar N1Crit, N12Crit; + TacsScalar D11 = D[0], D12 = D[1], D22 = D[3], D66 = D[5], + L = this->stiffenerPitch; - // Add the sensitivity of the stiffener failure value w.r.t. the DVs - // due to the dependence of the stiffener strains on the DVs - TacsScalar stiffenerFailStrainSens[TACSBeamConstitutive::NUM_STRESSES]; - this->evalStiffenerFailureStrainSens(stiffenerStrain, - stiffenerFailStrainSens); - this->addStrainTransformProductDVsens(stiffenerFailStrainSens, strain, - scale * dKSdf[1], dfdx); - } + // Create arrays for the sensitivities of the critical loads: + // [dN/dD11, dNdD22, dNdD12, dNdD66, dNdL] + TacsScalar N1CritSens[5], N12CritSens[5]; + N1Crit = this->computeCriticalLocalAxialLoadSens( + D11, D22, D12, D66, L, &N1CritSens[0], &N1CritSens[1], &N1CritSens[2], + &N1CritSens[3], &N1CritSens[4]); + N12Crit = this->computeCriticalShearLoadSens( + D11, D22, D12 + 2.0 * D66, L, &N12CritSens[0], &N12CritSens[1], + &N12CritSens[2], &N12CritSens[4]); - // --- Local panel buckling --- - if (this->includeLocalBuckling) { - this->addLocalPanelBucklingDVSens(elemIndex, scale * dKSdf[2], pt, X, - strain, dvLen, dfdx); + // N12CritSens[2] is currently dN12Crit/d(D12 + 2D66) + N12CritSens[3] = 2.0 * N12CritSens[2]; + + // Compute the buckling criteria and it's sensitivities to the applied and + // critical loads + TacsScalar dfdN1Local, dfdN12Local, dfdN1CritLocal, dfdN12CritLocal; + fails[2] = this->bucklingEnvelopeSens(-panelStress[0], N1Crit, panelStress[2], + N12Crit, &dfdN1Local, &dfdN1CritLocal, + &dfdN12Local, &dfdN12CritLocal); + dfdN1Local *= dKSdf[2]; + dfdN12Local *= dKSdf[2]; + dfdN1CritLocal *= dKSdf[2]; + dfdN12CritLocal *= dKSdf[2]; + + // Convert sensitivity w.r.t applied loads into sensitivity w.r.t DVs + TacsScalar dfdPanelStress[NUM_STRESSES]; + memset(dfdPanelStress, 0, NUM_STRESSES * sizeof(TacsScalar)); + dfdPanelStress[0] = -dfdN1Local; + dfdPanelStress[2] = dfdN12Local; + this->addPanelStressDVSens(scale, strain, dfdPanelStress, + &dfdx[this->panelDVStartNum]); + + // Convert the sensitivities of the critical loads w.r.t the D matrix entries + // to sensitivities of the buckling failure criteria w.r.t the D matrix + // entries + TacsScalar localBucklingDMatSens[4]; + for (int ii = 0; ii < 4; ii++) { + localBucklingDMatSens[ii] = + N1CritSens[ii] * dfdN1CritLocal + N12CritSens[ii] * dfdN12CritLocal; } // --- Global buckling sensitivity --- - if (this->includeGlobalBuckling) { - this->addGlobalPanelBucklingDVSens(elemIndex, scale * dKSdf[3], pt, X, - strain, dvLen, dfdx); - } - - // --- Stiffener column buckling --- - if (this->includeStiffenerColumnBuckling) { - TacsScalar stiffenerStress[TACSBeamConstitutive::NUM_STRESSES]; - this->computeStiffenerStress(stiffenerStrain, stiffenerStress); - const TacsScalar columnBucklingLoad = - this->computeStiffenerColumnBucklingLoad(); - const TacsScalar stiffenerAxialLoad = -stiffenerStress[0]; - this->addStiffenerColumnBucklingDVSens(dKSdf[4] * scale, strain, - stiffenerStrain, stiffenerAxialLoad, - columnBucklingLoad, dfdx); - } - - // --- Stiffener crippling --- - if (this->includeStiffenerCrippling) { - // Direct dependence of the stiffener crippling on the DVs - this->addStiffenerCripplingDVSens(scale * dKSdf[5], stiffenerStrain, - &dfdx[this->stiffenerDVStartNum]); - // Sensitivity of the stiffener crippling due to effect of DVs on the - // stiffener strains - TacsScalar stiffenerCripplingStrainSens[TACSBeamConstitutive::NUM_STRESSES]; - this->evalStiffenerCripplingStrainSens(stiffenerStrain, - stiffenerCripplingStrainSens); - this->addStrainTransformProductDVsens(stiffenerCripplingStrainSens, strain, - scale * dKSdf[5], dfdx); + TacsScalar stress[NUM_STRESSES]; + this->evalStress(0, NULL, NULL, strain, stress); + TacsScalar dfdN1Global, dfdN12Global, dfdN1CritGlobal, dfdN12CritGlobal; + TacsScalar D1, D2, D3; + this->computeCriticalGlobalBucklingStiffness(&D1, &D2, &D3); + L = this->panelLength; + N1Crit = computeCriticalGlobalAxialLoad(D1, L); + N12Crit = this->computeCriticalShearLoad(D1, D2, D3, L); + + this->bucklingEnvelopeSens(-stress[0], N1Crit, stress[2], N12Crit, + &dfdN1Global, &dfdN1CritGlobal, &dfdN12Global, + &dfdN12CritGlobal); + dfdN1Global *= dKSdf[3]; + dfdN12Global *= dKSdf[3]; + dfdN1CritGlobal *= dKSdf[3]; + dfdN12CritGlobal *= dKSdf[3]; + + // Add the sensitivity of the buckling failure criteria due to the dependence + // of the applied loads on the DVs + TacsScalar dfdStress[NUM_STRESSES]; + memset(dfdStress, 0, NUM_STRESSES * sizeof(TacsScalar)); + dfdStress[0] = -dfdN1Global; + dfdStress[2] = dfdN12Global; + this->addStressDVSens(elemIndex, scale, pt, X, strain, dfdStress, numDV, + dfdx); + + // Propogate the sensitivity of the buckling failure criteria w.r.t the + // critical loads back to the DVs + TacsScalar dfdD1, dfdD2, dfdD3, dfdPanelLength; + this->computeCriticalShearLoadSens(D1, D2, D3, L, &dfdD1, &dfdD2, &dfdD3, + &dfdPanelLength); + dfdD1 *= dfdN12CritGlobal; + dfdD2 *= dfdN12CritGlobal; + dfdD3 *= dfdN12CritGlobal; + dfdPanelLength *= dfdN12CritGlobal; + // N1CritGlobal = M_PI * M_PI * D1 / (L * L); + dfdD1 += dfdN1CritGlobal * (M_PI * M_PI / (L * L)); + dfdPanelLength += dfdN1CritGlobal * (-2.0 * M_PI * M_PI * D1 / (L * L * L)); + + if (this->panelLengthNum >= 0) { + dfdx[this->panelLengthLocalNum] += scale * dfdPanelLength; + } + + TacsScalar globalBucklingspSens, globalBucklingtpSens, globalBucklinghsSens, + globalBucklingtsSens, globalBucklingQstiffSens[NUM_Q_ENTRIES], + globalBucklingQpanelSens[NUM_Q_ENTRIES]; + this->computeCriticalGlobalBucklingStiffnessSens( + dfdD1, dfdD2, dfdD3, &globalBucklingspSens, &globalBucklingtpSens, + &globalBucklinghsSens, &globalBucklingtsSens, globalBucklingQpanelSens, + globalBucklingQstiffSens); + + // --- Panel thickness sensitivity --- + TacsScalar t = this->panelThick; + if (this->panelThickNum >= 0) { + int dvNum = this->panelThickLocalNum; + // --- Local buckling contribution --- + TacsScalar dMatdt = 0.25 * t * t; // d/dt(t^3/12) = t^2/4 + for (int ii = 0; ii < this->numPanelPlies; ii++) { + TacsScalar* Q = &(this->panelQMats[ii * NUM_Q_ENTRIES]); + dfdx[dvNum] += + scale * dMatdt * this->panelPlyFracs[ii] * + (localBucklingDMatSens[0] * Q[0] + localBucklingDMatSens[1] * Q[3] + + localBucklingDMatSens[2] * Q[1] + localBucklingDMatSens[3] * Q[5]); + } + // --- Global buckling contribution --- + dfdx[dvNum] += scale * globalBucklingtpSens; + } + + // --- Panel Ply fraction sensitivities --- + for (int ii = 0; ii < this->numPanelPlies; ii++) { + if (this->panelPlyFracNums[ii] >= 0) { + int dvNum = this->panelPlyFracLocalNums[ii]; + // --- Local buckling contribution --- + TacsScalar DFactor = t * t * t / 12.0; + TacsScalar* Q = &(this->panelQMats[ii * NUM_Q_ENTRIES]); + // Do df/dx += df/dD11 * dD11/dx + df/dD22 * dD22/dx + df/dD12 * dD12/dx + + // df/dD66 * dD66/dx + dfdx[dvNum] += + scale * DFactor * + (localBucklingDMatSens[0] * Q[0] + localBucklingDMatSens[1] * Q[3] + + localBucklingDMatSens[2] * Q[1] + localBucklingDMatSens[3] * Q[5]); + + // --- Global buckling contribution --- + dfdx[dvNum] += scale * (globalBucklingQpanelSens[0] * Q[0] + + globalBucklingQpanelSens[1] * Q[1] + + globalBucklingQpanelSens[2] * Q[2] + + globalBucklingQpanelSens[3] * Q[3] + + globalBucklingQpanelSens[4] * Q[4] + + globalBucklingQpanelSens[5] * Q[5]); + } + } + + // --- Stiffener height contribution --- + if (this->stiffenerHeightNum >= 0) { + int dvNum = this->stiffenerHeightLocalNum; + // --- Global buckling contribution --- + dfdx[dvNum] += scale * globalBucklinghsSens; + } + + // --- Stiffener thickness contribution --- + if (this->stiffenerThickNum >= 0) { + int dvNum = this->stiffenerThickLocalNum; + // --- Global buckling contribution --- + dfdx[dvNum] += scale * globalBucklingtsSens; + } + + // --- Stiffener ply fraction contributions --- + for (int ii = 0; ii < this->numStiffenerPlies; ii++) { + if (this->stiffenerPlyFracNums[ii] >= 0) { + int dvNum = this->stiffenerPlyFracLocalNums[ii]; + TacsScalar* Q = &(this->stiffenerQMats[ii * NUM_Q_ENTRIES]); + + // --- Global buckling contribution --- + dfdx[dvNum] += scale * (globalBucklingQstiffSens[0] * Q[0] + + globalBucklingQstiffSens[1] * Q[1] + + globalBucklingQstiffSens[2] * Q[2] + + globalBucklingQstiffSens[3] * Q[3] + + globalBucklingQstiffSens[4] * Q[4] + + globalBucklingQstiffSens[5] * Q[5]); + } + } + + // --- Stiffener pitch sensitivity --- + if (this->stiffenerPitchNum >= 0) { + dfdx[this->stiffenerPitchLocalNum] += + scale * (N12CritSens[4] * dfdN12CritLocal + + N1CritSens[4] * dfdN1CritLocal + globalBucklingspSens); } } @@ -1178,6 +1324,8 @@ void TACSBladeStiffenedShellConstitutive::transformStrain( void TACSBladeStiffenedShellConstitutive::transformStrainSens( const TacsScalar stiffenerStrainSens[], TacsScalar panelStrainSens[]) { + memset(panelStrainSens, 0, + TACSBeamConstitutive::NUM_STRESSES * sizeof(TacsScalar)); TacsScalar z = this->computeStiffenerCentroidHeight() - 0.5 * this->panelThick; @@ -1706,11 +1854,11 @@ void TACSBladeStiffenedShellConstitutive::computeStiffenerStiffness( void TACSBladeStiffenedShellConstitutive::computeEffectiveModulii( const int numPlies, const TacsScalar QMats[], const TacsScalar plyFracs[], - TacsScalar* const E, TacsScalar* const G) { + TacsScalar* E, TacsScalar* G) { *E = 0.; *G = 0.; for (int plyNum = 0; plyNum < numPlies; plyNum++) { - const TacsScalar* const Q = &(QMats[plyNum * NUM_Q_ENTRIES]); + const TacsScalar* Q = &(QMats[plyNum * NUM_Q_ENTRIES]); *E += plyFracs[plyNum] * (Q[0] - Q[1] * Q[1] / Q[3]); *G += plyFracs[plyNum] * Q[5]; } @@ -1719,24 +1867,21 @@ void TACSBladeStiffenedShellConstitutive::computeEffectiveModulii( // code above, but for some reason (probably related to floating point // arithmetic), it produces results that don't quite match complex step // derivatives w.r.t the ply fractions - // TacsScalar Q[NUM_Q_ENTRIES]; //, ABar[this->NUM_ABAR_ENTRIES]; - // for (int ii = 0; ii < NUM_Q_ENTRIES; ii++) { - // Q[ii] = 0.0; - // } - // for (int plyNum = 0; plyNum < numPlies; plyNum++) { - // const TacsScalar* Qply = &(QMats[plyNum * NUM_Q_ENTRIES]); - // for (int ii = 0; ii < NUM_Q_ENTRIES; ii++) { - // Q[ii] += plyFracs[plyNum] * Qply[ii]; - // } - // } - - // // Compute the effective elastic and shear moduli - // TacsScalar Q11 = Q[0]; - // TacsScalar Q12 = Q[1]; - // TacsScalar Q22 = Q[3]; - // TacsScalar Q66 = Q[5]; - // *E = Q11 - Q12 * Q12 / Q22; - // *G = Q66; + /* + TacsScalar Q[this->NUM_Q_ENTRIES], ABar[this->NUM_ABAR_ENTRIES]; + + this->computeSmearedStiffness(this->numStiffenerPlies, this->stiffenerQMats, + this->stiffenerAbarMats, + this->stiffenerPlyFracs, Q, ABar); + + // Compute the effective elastic and shear moduli + TacsScalar Q11 = Q[0]; + TacsScalar Q12 = Q[1]; + TacsScalar Q22 = Q[3]; + TacsScalar Q66 = Q[5]; + E = Q11 - Q12 * Q12 / Q22; + G = Q66; + */ } // Compute the failure criteria for the stiffener @@ -1745,32 +1890,21 @@ TacsScalar TACSBladeStiffenedShellConstitutive::computeStiffenerFailure( TACSOrthotropicPly* ply = this->stiffenerPly; // Compute the strain state at the tip of the stiffener - TacsScalar zCentroid = this->computeStiffenerCentroidHeight(); - TacsScalar zTipOffset = - -(this->stiffenerHeight + this->stiffenerThick) - zCentroid; - TacsScalar plyStrain[3]; - memset(plyStrain, 0, 3 * sizeof(TacsScalar)); - plyStrain[0] = stiffenerStrain[0] + zTipOffset * stiffenerStrain[2]; + TacsScalar zTipOffset = -(this->stiffenerHeight + this->stiffenerThick) - + this->computeStiffenerCentroidHeight(); + TacsScalar tipStrain[3]; + memset(tipStrain, 0, 3 * sizeof(TacsScalar)); + tipStrain[0] = stiffenerStrain[0] + zTipOffset * stiffenerStrain[2]; // Compute the failure criteria at this strain state for each ply angle for (int ii = 0; ii < this->numStiffenerPlies; ii++) { this->stiffenerPlyFailValues[ii] = - ply->failure(this->stiffenerPlyAngles[ii], plyStrain); - } - - // Now do the same thing for the bottom of the stiffener - zTipOffset = -zCentroid; - plyStrain[0] = stiffenerStrain[0] + zTipOffset * stiffenerStrain[2]; - - // Compute the failure criteria at this strain state for each ply angle - for (int ii = 0; ii < this->numStiffenerPlies; ii++) { - this->stiffenerPlyFailValues[this->numStiffenerPlies + ii] = - ply->failure(this->stiffenerPlyAngles[ii], plyStrain); + ply->failure(this->stiffenerPlyAngles[ii], tipStrain); } // Returned the aggregated value over all plies - return ksAggregation(this->stiffenerPlyFailValues, - 2 * this->numStiffenerPlies, this->ksWeight); + return ksAggregation(this->stiffenerPlyFailValues, this->numStiffenerPlies, + this->ksWeight); } TacsScalar TACSBladeStiffenedShellConstitutive::evalStiffenerFailureStrainSens( @@ -1785,15 +1919,11 @@ TacsScalar TACSBladeStiffenedShellConstitutive::evalStiffenerFailureStrainSens( memset(sens, 0, numStrain * sizeof(TacsScalar)); // Compute the strain state at the tip of the stiffener - TacsScalar zCentroid = this->computeStiffenerCentroidHeight(); - TacsScalar zTipOffset = - -(this->stiffenerHeight + this->stiffenerThick) - zCentroid; + TacsScalar zTipOffset = -(this->stiffenerHeight + this->stiffenerThick) - + this->computeStiffenerCentroidHeight(); TacsScalar tipStrain[3]; memset(tipStrain, 0, 3 * sizeof(TacsScalar)); tipStrain[0] = strain[0] + zTipOffset * strain[2]; - TacsScalar bottomStrain[3]; - memset(bottomStrain, 0, 3 * sizeof(TacsScalar)); - bottomStrain[0] = strain[0] - zCentroid * strain[2]; // Compute the failure criteria at this strain state for each ply angle and // the sensitivity of the failure criteria w.r.t the strains @@ -1801,22 +1931,15 @@ TacsScalar TACSBladeStiffenedShellConstitutive::evalStiffenerFailureStrainSens( TacsScalar plyFailStrainSens[3]; fails[ii] = ply->failureStrainSens(angles[ii], tipStrain, plyFailStrainSens); - // Convert the sensitivity w.r.t the tip strains to the sensitivity - // w.r.t beam strains + // Convert the sensitivity w.r.t the tip strain to the sensitivity w.r.t + // beam strains memset(dFaildStrain[ii], 0, numStrain * sizeof(TacsScalar)); dFaildStrain[ii][0] = plyFailStrainSens[0]; dFaildStrain[ii][2] = zTipOffset * plyFailStrainSens[0]; - - // Now do the same for the bottom of the stiffener - fails[numPlies + ii] = - ply->failureStrainSens(angles[ii], bottomStrain, plyFailStrainSens); - memset(dFaildStrain[numPlies + ii], 0, numStrain * sizeof(TacsScalar)); - dFaildStrain[numPlies + ii][0] = plyFailStrainSens[0]; - dFaildStrain[numPlies + ii][2] = -zCentroid * plyFailStrainSens[0]; } TacsScalar fail = ksAggregationSensProduct( - fails, 2 * numPlies, numStrain, this->ksWeight, dFaildStrain, sens); + fails, numPlies, numStrain, this->ksWeight, dFaildStrain, sens); return fail; } @@ -1826,18 +1949,16 @@ void TACSBladeStiffenedShellConstitutive::addStiffenerFailureDVSens( TACSOrthotropicPly* ply = this->stiffenerPly; TacsScalar* fails = this->stiffenerPlyFailValues; const TacsScalar* angles = this->stiffenerPlyAngles; - TacsScalar* dKSdFail = this->stiffenerPlyFailSens; + TacsScalar* dKSdFail = this->panelPlyFailSens; const int hNum = this->stiffenerHeightLocalNum - this->stiffenerDVStartNum; const int tNum = this->stiffenerThickLocalNum - this->stiffenerDVStartNum; - TacsScalar zCentroid = this->computeStiffenerCentroidHeight(); - TacsScalar zTipOffset = - -(this->stiffenerHeight + this->stiffenerThick) - zCentroid; - TacsScalar dZdh, dZdt; - this->computeStiffenerCentroidHeightSens(dZdt, dZdh); + TacsScalar zTipOffset = -(this->stiffenerHeight + this->stiffenerThick) - + this->computeStiffenerCentroidHeight(); TacsScalar dTipStraindh, dTipStraindt; - dTipStraindh = -dZdh - 1.0; - dTipStraindt = -dZdt - 1.0; + this->computeStiffenerCentroidHeightSens(dTipStraindt, dTipStraindh); + dTipStraindh = -dTipStraindh - 1.0; + dTipStraindt = -dTipStraindt - 1.0; // Compute the strain state at the tip of the stiffener TacsScalar tipStrain[3]; @@ -1851,36 +1972,19 @@ void TACSBladeStiffenedShellConstitutive::addStiffenerFailureDVSens( fails[ii] = ply->failure(angles[ii], tipStrain); } - // Now do the same thing for the bottom of the stiffener - TacsScalar bottomStrain[3]; - memset(bottomStrain, 0, 3 * sizeof(TacsScalar)); - zTipOffset = -zCentroid; - bottomStrain[0] = strain[0] + zTipOffset * strain[2]; - - for (int ii = 0; ii < this->numStiffenerPlies; ii++) { - fails[this->numStiffenerPlies + ii] = - ply->failure(this->stiffenerPlyAngles[ii], bottomStrain); - } - // Compute the sensitivity of the KS aggregation w.r.t the failure values - ksAggregationSens(fails, 2 * this->numStiffenerPlies, this->ksWeight, - dKSdFail); + ksAggregationSens(fails, this->numStiffenerPlies, this->ksWeight, dKSdFail); // Now go back through each ply, compute the strain sensitivity of it's // failure, then convert it to a DV sensitivity and add it to the dfdx array for (int ii = 0; ii < this->numStiffenerPlies; ii++) { - TacsScalar tipFailStrainSens[3], bottomFailStrainSens[3]; - ply->failureStrainSens(angles[ii], tipStrain, tipFailStrainSens); - ply->failureStrainSens(angles[ii], bottomStrain, bottomFailStrainSens); + TacsScalar plyFailStrainSens[3]; + ply->failureStrainSens(angles[ii], tipStrain, plyFailStrainSens); if (hNum >= 0) { - dfdx[hNum] += scale * dKSdFail[ii] * dTipStraindh * tipFailStrainSens[0]; - dfdx[hNum] += scale * dKSdFail[this->numStiffenerPlies + ii] * -dZdh * - strain[2] * bottomFailStrainSens[0]; + dfdx[hNum] += scale * dKSdFail[ii] * dTipStraindh * plyFailStrainSens[0]; } if (tNum >= 0) { - dfdx[tNum] += scale * dKSdFail[ii] * dTipStraindt * tipFailStrainSens[0]; - dfdx[tNum] += scale * dKSdFail[this->numStiffenerPlies + ii] * -dZdt * - strain[2] * bottomFailStrainSens[0]; + dfdx[tNum] += scale * dKSdFail[ii] * dTipStraindt * plyFailStrainSens[0]; } } } @@ -2016,24 +2120,9 @@ void TACSBladeStiffenedShellConstitutive::computeStiffenerMOISens( // Buckling functions // ============================================================================== -TacsScalar TACSBladeStiffenedShellConstitutive::evalGlobalPanelBuckling( - const TacsScalar e[]) { - TacsScalar stress[TACSShellConstitutive::NUM_STRESSES]; - TacsScalar D1, D2, D3; - this->computeCriticalGlobalBucklingStiffness(&D1, &D2, &D3); - const TacsScalar L = this->panelLength; - - const TacsScalar N1Crit = computeCriticalGlobalAxialLoad(D1, L); - const TacsScalar N12Crit = this->computeCriticalShearLoad(D1, D2, D3, L); - - this->evalStress(0, NULL, NULL, e, stress); - return this->bucklingEnvelope(-stress[0], N1Crit, stress[2], N12Crit); -} - void TACSBladeStiffenedShellConstitutive:: - computeCriticalGlobalBucklingStiffness(TacsScalar* const D1, - TacsScalar* const D2, - TacsScalar* const D3) { + computeCriticalGlobalBucklingStiffness(TacsScalar* D1, TacsScalar* D2, + TacsScalar* D3) { // Compute the Q matrices for the panel and stiffener TacsScalar QPanel[this->NUM_Q_ENTRIES], QStiffener[this->NUM_Q_ENTRIES]; this->computeSmearedStiffness(this->numPanelPlies, this->panelQMats, @@ -2066,15 +2155,14 @@ void TACSBladeStiffenedShellConstitutive:: // --- 2-direction bending stiffness --- // TacsScalar D2Panel = (tp3 * QPanel[3]) / 12.0; - // TacsScalar D2Stiff = ((tp + 0.5 * ts) * (tp + 0.5 * ts) * (tp + 0.5 * ts) - // * + // TacsScalar D2Stiff = ((tp + 0.5 * ts) * (tp + 0.5 * ts) * (tp + 0.5 * ts) * // (QPanel[3] + QStiffener[3])) / // 6.0; // *D2 = ps / ((ps - hs * kf) / D2Panel + (hs * kf) / D2Stiff); // NOTE: I am ignoring the contribution of the stiffener flange to the // 2-direction bending stiffness so that the stiffness used for the buckling - // calculation is consistent with the stiffness matrix. Not sure whether - // this is a good idea or not. + // calculation is consistent with the stiffness matrix. Not sure whether this + // is a good idea or not. *D2 = tp3 / 12.0 * QPanel[3]; // --- Twisting stiffness --- @@ -2084,158 +2172,15 @@ void TACSBladeStiffenedShellConstitutive:: *D3 = (QPanel[5] * (Jp + Ap * zg * zg) + 0.25 * QStiffener[5] * (Js + As * (zg - zs) * (zg - zs))) / ps; -} - -TacsScalar -TACSBladeStiffenedShellConstitutive::evalGlobalPanelBucklingStrainSens( - const TacsScalar e[], TacsScalar sens[]) { - TacsScalar stiffness[NUM_TANGENT_STIFFNESS_ENTRIES], stress[NUM_STRESSES]; - this->computeStiffness(stiffness); - const TacsScalar *A, *B, *D, *As; - TacsScalar drill; - this->extractTangentStiffness(stiffness, &A, &B, &D, &As, &drill); - this->computeStress(A, B, D, As, drill, e, stress); - TacsScalar N1GlobalSens, N1CritGlobalSens, N12GlobalSens, N12CritGlobalSens; - TacsScalar D1, D2, D3; - this->computeCriticalGlobalBucklingStiffness(&D1, &D2, &D3); - const TacsScalar L = this->panelLength; - TacsScalar N1CritGlobal = computeCriticalGlobalAxialLoad(D1, L); - TacsScalar N12CritGlobal = this->computeCriticalShearLoad(D1, D2, D3, L); - - const TacsScalar strengthRatio = this->bucklingEnvelopeSens( - -stress[0], N1CritGlobal, stress[2], N12CritGlobal, &N1GlobalSens, - &N1CritGlobalSens, &N12GlobalSens, &N12CritGlobalSens); - - sens[0] += N1GlobalSens * -A[0] + N12GlobalSens * A[2]; - sens[1] += N1GlobalSens * -A[1] + N12GlobalSens * A[4]; - sens[2] += N1GlobalSens * -A[2] + N12GlobalSens * A[5]; - sens[3] += N1GlobalSens * -B[0] + N12GlobalSens * B[2]; - sens[4] += N1GlobalSens * -B[1] + N12GlobalSens * B[4]; - sens[5] += N1GlobalSens * -B[2] + N12GlobalSens * B[5]; - - return strengthRatio; -} - -void TACSBladeStiffenedShellConstitutive::addGlobalPanelBucklingDVSens( - int elemIndex, TacsScalar scale, const double pt[], const TacsScalar X[], - const TacsScalar strain[], int dvLen, TacsScalar dfdx[]) { - TacsScalar stress[NUM_STRESSES]; - this->evalStress(0, NULL, NULL, strain, stress); - TacsScalar dfdN1Global, dfdN12Global, dfdN1CritGlobal, dfdN12CritGlobal; - TacsScalar D1, D2, D3; - TacsScalar N1Crit, N12Crit; - this->computeCriticalGlobalBucklingStiffness(&D1, &D2, &D3); - const TacsScalar L = this->panelLength; - N1Crit = computeCriticalGlobalAxialLoad(D1, L); - N12Crit = this->computeCriticalShearLoad(D1, D2, D3, L); - - this->bucklingEnvelopeSens(-stress[0], N1Crit, stress[2], N12Crit, - &dfdN1Global, &dfdN1CritGlobal, &dfdN12Global, - &dfdN12CritGlobal); - // dfdN1Global *= dKSdf[3]; - // dfdN12Global *= dKSdf[3]; - // dfdN1CritGlobal *= dKSdf[3]; - // dfdN12CritGlobal *= dKSdf[3]; - - // Add the sensitivity of the buckling failure criteria due to the - // dependence of the applied loads on the DVs - TacsScalar dfdStress[NUM_STRESSES]; - memset(dfdStress, 0, NUM_STRESSES * sizeof(TacsScalar)); - dfdStress[0] = -dfdN1Global; - dfdStress[2] = dfdN12Global; - this->addStressDVSens(elemIndex, scale, pt, X, strain, dfdStress, dvLen, - dfdx); - - // Propogate the sensitivity of the buckling failure criteria w.r.t the - // critical loads back to the DVs - TacsScalar dfdD1, dfdD2, dfdD3, dfdPanelLength; - this->computeCriticalShearLoadSens(D1, D2, D3, L, &dfdD1, &dfdD2, &dfdD3, - &dfdPanelLength); - dfdD1 *= dfdN12CritGlobal; - dfdD2 *= dfdN12CritGlobal; - dfdD3 *= dfdN12CritGlobal; - dfdPanelLength *= dfdN12CritGlobal; - // N1CritGlobal = M_PI * M_PI * D1 / (L * L); - dfdD1 += dfdN1CritGlobal * (M_PI * M_PI / (L * L)); - dfdPanelLength += dfdN1CritGlobal * (-2.0 * M_PI * M_PI * D1 / (L * L * L)); - - if (this->panelLengthNum >= 0) { - dfdx[this->panelLengthLocalNum] += scale * dfdPanelLength; - } - - TacsScalar globalBucklingspSens, globalBucklingtpSens, globalBucklinghsSens, - globalBucklingtsSens, globalBucklingQstiffSens[NUM_Q_ENTRIES], - globalBucklingQpanelSens[NUM_Q_ENTRIES]; - this->computeCriticalGlobalBucklingStiffnessSens( - dfdD1, dfdD2, dfdD3, &globalBucklingspSens, &globalBucklingtpSens, - &globalBucklinghsSens, &globalBucklingtsSens, globalBucklingQpanelSens, - globalBucklingQstiffSens); - - // --- Panel thickness sensitivity --- - TacsScalar t = this->panelThick; - if (this->panelThickNum >= 0) { - const int dvNum = this->panelThickLocalNum; - dfdx[dvNum] += scale * globalBucklingtpSens; - } - - // --- Stiffener height contribution --- - if (this->stiffenerHeightNum >= 0) { - const int dvNum = this->stiffenerHeightLocalNum; - // --- Global buckling contribution --- - dfdx[dvNum] += scale * globalBucklinghsSens; - } - - // --- Stiffener thickness contribution --- - if (this->stiffenerThickNum >= 0) { - const int dvNum = this->stiffenerThickLocalNum; - // --- Global buckling contribution --- - dfdx[dvNum] += scale * globalBucklingtsSens; - } - - // --- Panel Ply fraction sensitivities --- - for (int ii = 0; ii < this->numPanelPlies; ii++) { - if (this->panelPlyFracNums[ii] >= 0) { - const int dvNum = this->panelPlyFracLocalNums[ii]; - const TacsScalar* const Q = &(this->panelQMats[ii * NUM_Q_ENTRIES]); - - // --- Global buckling contribution --- - dfdx[dvNum] += scale * (globalBucklingQpanelSens[0] * Q[0] + - globalBucklingQpanelSens[1] * Q[1] + - globalBucklingQpanelSens[2] * Q[2] + - globalBucklingQpanelSens[3] * Q[3] + - globalBucklingQpanelSens[4] * Q[4] + - globalBucklingQpanelSens[5] * Q[5]); - } - } - // --- Stiffener ply fraction contributions --- - for (int ii = 0; ii < this->numStiffenerPlies; ii++) { - if (this->stiffenerPlyFracNums[ii] >= 0) { - const int dvNum = this->stiffenerPlyFracLocalNums[ii]; - const TacsScalar* const Q = &(this->stiffenerQMats[ii * NUM_Q_ENTRIES]); - - // --- Global buckling contribution --- - dfdx[dvNum] += scale * (globalBucklingQstiffSens[0] * Q[0] + - globalBucklingQstiffSens[1] * Q[1] + - globalBucklingQstiffSens[2] * Q[2] + - globalBucklingQstiffSens[3] * Q[3] + - globalBucklingQstiffSens[4] * Q[4] + - globalBucklingQstiffSens[5] * Q[5]); - } - } - - // --- Stiffener pitch sensitivity --- - if (this->stiffenerPitchNum >= 0) { - dfdx[this->stiffenerPitchLocalNum] += scale * globalBucklingspSens; - } + TacsScalar C[NUM_TANGENT_STIFFNESS_ENTRIES]; } void TACSBladeStiffenedShellConstitutive:: computeCriticalGlobalBucklingStiffnessSens( const TacsScalar dfdD1, const TacsScalar dfdD2, const TacsScalar dfdD3, - TacsScalar* const psSens, TacsScalar* const tpSens, - TacsScalar* const hsSens, TacsScalar* const tsSens, - TacsScalar QpanelSens[], TacsScalar QstiffSens[]) { + TacsScalar* psSens, TacsScalar* tpSens, TacsScalar* hsSens, + TacsScalar* tsSens, TacsScalar QpanelSens[], TacsScalar QstiffSens[]) { // Zero the sensitivities *tpSens = 0.0; *psSens = 0.0; @@ -2275,8 +2220,8 @@ void TACSBladeStiffenedShellConstitutive:: TacsScalar zns2 = (zn - zs) * (zn - zs); TacsScalar pInv = 1.0 / ps; - // D1 = (E1p * (Ip + Ap * zn * zn) + E1s * (Is + As * (zn - zs) * (zn - - // zs))) / ps; + // D1 = (E1p * (Ip + Ap * zn * zn) + E1s * (Is + As * (zn - zs) * (zn - zs))) + // / ps; TacsScalar E1pSens = dfdD1 * ((Ap * zn2 + Ip) / ps); TacsScalar IpSens = dfdD1 * (E1p / ps); @@ -2400,11 +2345,7 @@ void TACSBladeStiffenedShellConstitutive:: } void TACSBladeStiffenedShellConstitutive::testGlobalBucklingStiffnessSens() { -#ifdef TACS_USE_COMPLEX - const double eps = 1e-200; -#else - const double eps = 1e-6; -#endif + const TacsScalar eps = 1e-6; // Get the number of DVs int num_dvs = this->getDesignVars(0, 0, NULL); TacsScalar* DV0 = new TacsScalar[num_dvs]; @@ -2439,7 +2380,7 @@ void TACSBladeStiffenedShellConstitutive::testGlobalBucklingStiffnessSens() { TacsRealPart(psSens[2])); TacsScalar D1Pert, D2Pert, D3Pert; #ifdef TACS_USE_COMPLEX - DVPert[this->stiffenerPitchLocalNum] += TacsScalar(0.0, eps); + DVPert[this->stiffenerPitchLocalNum] += TacsScalar(0.0, 1e-200); #else DVPert[this->stiffenerPitchLocalNum] += eps; #endif @@ -2447,7 +2388,7 @@ void TACSBladeStiffenedShellConstitutive::testGlobalBucklingStiffnessSens() { this->computeCriticalGlobalBucklingStiffness(&D1Pert, &D2Pert, &D3Pert); #ifdef TACS_USE_COMPLEX - TacsScalar psSensFD = TacsImagPart(D1Pert) / eps; + TacsScalar psSensFD = TacsImagPart(D1Pert) * 1e200; #else TacsScalar psSensFD = (D1Pert - D1) / eps; #endif @@ -2459,7 +2400,7 @@ void TACSBladeStiffenedShellConstitutive::testGlobalBucklingStiffnessSens() { printf(" FD: D1sens = % 011.7e, ", TacsRealPart(psSensFD)); #ifdef TACS_USE_COMPLEX - psSensFD = TacsImagPart(D2Pert) / eps; + psSensFD = TacsImagPart(D2Pert) * 1e200; #else psSensFD = (D2Pert - D2) / eps; #endif @@ -2663,26 +2604,23 @@ void TACSBladeStiffenedShellConstitutive::testGlobalBucklingStiffnessSens() { delete[] DVPert; } -TacsScalar TACSBladeStiffenedShellConstitutive::evalLocalPanelBuckling( - const TacsScalar e[]) { - // Compute panel stiffness matrix and loads - TacsScalar panelStiffness[NUM_TANGENT_STIFFNESS_ENTRIES], - stress[NUM_STRESSES]; - this->computePanelStiffness(panelStiffness); - const TacsScalar *A, *D; - this->extractTangentStiffness(panelStiffness, &A, NULL, &D, NULL, NULL); - this->computePanelStress(e, stress); +TacsScalar +TACSBladeStiffenedShellConstitutive::computeCriticalLocalAxialLoadSens( + const TacsScalar D11, const TacsScalar D22, const TacsScalar D12, + const TacsScalar D66, const TacsScalar L, TacsScalar* D11Sens, + TacsScalar* D22Sens, TacsScalar* D12Sens, TacsScalar* D66Sens, + TacsScalar* LSens) { + double pi2 = M_PI * M_PI; + TacsScalar L2 = L * L; + TacsScalar root = sqrt(D11 * D22); - // Compute the critical local loads - const TacsScalar D11 = D[0], D12 = D[1], D22 = D[3], D66 = D[5], - L = this->stiffenerPitch; - const TacsScalar N1Crit = - this->computeCriticalLocalAxialLoad(D11, D22, D12, D66, L); - const TacsScalar N12Crit = - this->computeCriticalShearLoad(D11, D22, D12 + 2.0 * D66, L); + *D11Sens = pi2 * root / (D11 * L2); + *D12Sens = 2.0 * pi2 / L2; + *D22Sens = pi2 * root / (D22 * L2); + *D66Sens = 4.0 * pi2 / L2; + *LSens = 4.0 * pi2 * (-D12 - 2.0 * D66 - root) / (L2 * L); - // Compute the buckling criteria - return this->bucklingEnvelope(-stress[0], N1Crit, stress[2], N12Crit); + return 2.0 * pi2 / L2 * (root + D12 + 2.0 * D66); } TacsScalar TACSBladeStiffenedShellConstitutive::computeCriticalShearLoad( @@ -2708,166 +2646,15 @@ TacsScalar TACSBladeStiffenedShellConstitutive::computeCriticalShearLoad( return N12_crit; } -TacsScalar -TACSBladeStiffenedShellConstitutive::evalLocalPanelBucklingStrainSens( - const TacsScalar e[], TacsScalar sens[]) { - // Compute panel stiffness matrix and loads - TacsScalar panelStiffness[NUM_TANGENT_STIFFNESS_ENTRIES], - panelStress[NUM_STRESSES]; - this->computePanelStiffness(panelStiffness); - const TacsScalar *APanel, *DPanel; - this->extractTangentStiffness(panelStiffness, &APanel, NULL, &DPanel, NULL, - NULL); - this->computePanelStress(e, panelStress); - - // Compute the critical local loads (no need to compute their - // sensitivities because they're not dependent on the strain)) - TacsScalar N1CritLocal, N12CritLocal; - TacsScalar D11 = DPanel[0], D12 = DPanel[1], D22 = DPanel[3], D66 = DPanel[5], - L = this->stiffenerPitch; - N1CritLocal = this->computeCriticalLocalAxialLoad(D11, D22, D12, D66, L); - N12CritLocal = this->computeCriticalShearLoad(D11, D22, D12 + 2.0 * D66, L); - - // Compute the buckling criteria and it's sensitivities - TacsScalar N1LocalSens, N12LocalSens, N1CritLocalSens, N12CritLocalSens; - const TacsScalar strengthRatio = this->bucklingEnvelopeSens( - -panelStress[0], N1CritLocal, panelStress[2], N12CritLocal, &N1LocalSens, - &N1CritLocalSens, &N12LocalSens, &N12CritLocalSens); - - sens[0] = N1LocalSens * -APanel[0] + N12LocalSens * APanel[2]; - sens[1] = N1LocalSens * -APanel[1] + N12LocalSens * APanel[4]; - sens[2] = N1LocalSens * -APanel[2] + N12LocalSens * APanel[5]; - for (int ii = 3; ii < NUM_STRESSES; ii++) { - sens[ii] = 0.0; - } - - return strengthRatio; -} - -void TACSBladeStiffenedShellConstitutive::addLocalPanelBucklingDVSens( - int elemIndex, TacsScalar scale, const double pt[], const TacsScalar X[], - const TacsScalar strain[], int dvLen, TacsScalar dfdx[]) { - // Compute panel stiffness matrix and loads - TacsScalar panelStiffness[NUM_TANGENT_STIFFNESS_ENTRIES], - panelStress[NUM_STRESSES]; - const TacsScalar t = this->panelThick; - this->computePanelStiffness(panelStiffness); - const TacsScalar *A, *D; - this->extractTangentStiffness(panelStiffness, &A, NULL, &D, NULL, NULL); - this->computePanelStress(strain, panelStress); - - // Compute the critical local loads and their sensitivities - TacsScalar N1Crit, N12Crit; - const TacsScalar D11 = D[0], D12 = D[1], D22 = D[3], D66 = D[5], - L = this->stiffenerPitch; - - // Create arrays for the sensitivities of the critical loads: - // [dN/dD11, dNdD22, dNdD12, dNdD66, dNdL] - TacsScalar N1CritSens[5], N12CritSens[5]; - N1Crit = this->computeCriticalLocalAxialLoadSens( - D11, D22, D12, D66, L, &N1CritSens[0], &N1CritSens[1], &N1CritSens[2], - &N1CritSens[3], &N1CritSens[4]); - N12Crit = this->computeCriticalShearLoadSens( - D11, D22, D12 + 2.0 * D66, L, &N12CritSens[0], &N12CritSens[1], - &N12CritSens[2], &N12CritSens[4]); - - // N12CritSens[2] is currently dN12Crit/d(D12 + 2D66) - N12CritSens[3] = 2.0 * N12CritSens[2]; - - // Compute the buckling criteria and it's sensitivities to the applied and - // critical loads - TacsScalar dfdN1Local, dfdN12Local, dfdN1CritLocal, dfdN12CritLocal; - this->bucklingEnvelopeSens(-panelStress[0], N1Crit, panelStress[2], N12Crit, - &dfdN1Local, &dfdN1CritLocal, &dfdN12Local, - &dfdN12CritLocal); - // dfdN1Local *= dKSdf[2]; - // dfdN12Local *= dKSdf[2]; - // dfdN1CritLocal *= dKSdf[2]; - // dfdN12CritLocal *= dKSdf[2]; - - // Convert sensitivity w.r.t applied loads into sensitivity w.r.t DVs - TacsScalar dfdPanelStress[NUM_STRESSES]; - memset(dfdPanelStress, 0, NUM_STRESSES * sizeof(TacsScalar)); - dfdPanelStress[0] = -dfdN1Local; - dfdPanelStress[2] = dfdN12Local; - this->addPanelStressDVSens(scale, strain, dfdPanelStress, - &dfdx[this->panelDVStartNum]); - - // Convert the sensitivities of the critical loads w.r.t the D matrix - // entries to sensitivities of the buckling failure criteria w.r.t the D - // matrix entries - TacsScalar localBucklingDMatSens[4]; - for (int ii = 0; ii < 4; ii++) { - localBucklingDMatSens[ii] = - N1CritSens[ii] * dfdN1CritLocal + N12CritSens[ii] * dfdN12CritLocal; - } - - // --- Panel thickness sensitivity --- - if (this->panelThickNum >= 0) { - const int dvNum = this->panelThickLocalNum; - // --- Local buckling contribution --- - const TacsScalar dMatdt = 0.25 * t * t; // d/dt(t^3/12) = t^2/4 - for (int ii = 0; ii < this->numPanelPlies; ii++) { - TacsScalar* Q = &(this->panelQMats[ii * NUM_Q_ENTRIES]); - dfdx[dvNum] += - scale * dMatdt * this->panelPlyFracs[ii] * - (localBucklingDMatSens[0] * Q[0] + localBucklingDMatSens[1] * Q[3] + - localBucklingDMatSens[2] * Q[1] + localBucklingDMatSens[3] * Q[5]); - } - } - - // --- Panel Ply fraction sensitivities --- - for (int ii = 0; ii < this->numPanelPlies; ii++) { - if (this->panelPlyFracNums[ii] >= 0) { - const int dvNum = this->panelPlyFracLocalNums[ii]; - // --- Local buckling contribution --- - const TacsScalar DFactor = t * t * t / 12.0; - const TacsScalar* const Q = &(this->panelQMats[ii * NUM_Q_ENTRIES]); - // Do df/dx += df/dD11 * dD11/dx + df/dD22 * dD22/dx + df/dD12 * dD12/dx - // + df/dD66 * dD66/dx - dfdx[dvNum] += - scale * DFactor * - (localBucklingDMatSens[0] * Q[0] + localBucklingDMatSens[1] * Q[3] + - localBucklingDMatSens[2] * Q[1] + localBucklingDMatSens[3] * Q[5]); - } - } - - // --- Stiffener pitch sensitivity --- - if (this->stiffenerPitchNum >= 0) { - dfdx[this->stiffenerPitchLocalNum] += - scale * - (N12CritSens[4] * dfdN12CritLocal + N1CritSens[4] * dfdN1CritLocal); - } -} - -TacsScalar -TACSBladeStiffenedShellConstitutive::computeCriticalLocalAxialLoadSens( - const TacsScalar D11, const TacsScalar D22, const TacsScalar D12, - const TacsScalar D66, const TacsScalar L, TacsScalar* const D11Sens, - TacsScalar* const D22Sens, TacsScalar* const D12Sens, - TacsScalar* const D66Sens, TacsScalar* const LSens) { - const double pi2 = M_PI * M_PI; - const TacsScalar L2 = L * L; - const TacsScalar root = sqrt(D11 * D22); - - *D11Sens = pi2 * root / (D11 * L2); - *D12Sens = 2.0 * pi2 / L2; - *D22Sens = pi2 * root / (D22 * L2); - *D66Sens = 4.0 * pi2 / L2; - *LSens = 4.0 * pi2 * (-D12 - 2.0 * D66 - root) / (L2 * L); - - return 2.0 * pi2 / L2 * (root + D12 + 2.0 * D66); -} - TacsScalar TACSBladeStiffenedShellConstitutive::computeCriticalShearLoadSens( const TacsScalar D1, const TacsScalar D2, const TacsScalar D3, - const TacsScalar L, TacsScalar* const D1Sens, TacsScalar* const D2Sens, - TacsScalar* const D3Sens, TacsScalar* const LSens) { - const TacsScalar L2 = L * L; - const TacsScalar L3 = L2 * L; - const TacsScalar D32 = D3 * D3; - const TacsScalar root = sqrt(D1 * D2); - const TacsScalar xi = root / D3; + const TacsScalar L, TacsScalar* D1Sens, TacsScalar* D2Sens, + TacsScalar* D3Sens, TacsScalar* LSens) { + TacsScalar L2 = L * L; + TacsScalar L3 = L2 * L; + TacsScalar D32 = D3 * D3; + TacsScalar root = sqrt(D1 * D2); + TacsScalar xi = root / D3; TacsScalar N12CritVals[2]; N12CritVals[0] = -(4.0 / L2) * sqrt(D3 * D1 * xi) * (8.125 + 5.045 / xi); @@ -2878,32 +2665,32 @@ TacsScalar TACSBladeStiffenedShellConstitutive::computeCriticalShearLoadSens( N12Crit = ksAggregationSens(N12CritVals, 2, 50., N12CritSens); N12Crit *= -1.0; - const TacsScalar dN12Crit1_dD1 = + TacsScalar dN12Crit1_dD1 = sqrt(D1 * root) * (5.045 * D3 + 24.375 * root) / (D1 * L2 * root); - const TacsScalar dN12Crit2_dD1 = + TacsScalar dN12Crit2_dD1 = sqrt(D1 * D3) * (5.628 * D1 * D2 + 23.4 * D32 + 2.128 * D3 * root) / (D1 * D32 * L2); *D1Sens = dN12Crit1_dD1 * N12CritSens[0] + dN12Crit2_dD1 * N12CritSens[1]; - const TacsScalar dN12Crit1_dD2 = + TacsScalar dN12Crit1_dD2 = sqrt(D1 * root) * (-5.045 * D3 + 8.125 * root) / (D2 * L2 * root); - const TacsScalar dN12Crit2_dD2 = + TacsScalar dN12Crit2_dD2 = sqrt(D1 * D3) * (3.752 * D1 * D2 + 1.064 * D3 * root) / (D2 * D32 * L2); *D2Sens = dN12Crit1_dD2 * N12CritSens[0] + dN12Crit2_dD2 * N12CritSens[1]; - const TacsScalar dN12Crit1_dD3 = 20.18 * sqrt(D1 * root) / (L2 * root); - const TacsScalar dN12Crit2_dD3 = + TacsScalar dN12Crit1_dD3 = 20.18 * sqrt(D1 * root) / (L2 * root); + TacsScalar dN12Crit2_dD3 = sqrt(D1 * D3) * (-5.628 * D1 * D2 + 23.4 * D32 - 1.064 * D3 * root) / (D32 * D3 * L2); *D3Sens = dN12Crit1_dD3 * N12CritSens[0] + dN12Crit2_dD3 * N12CritSens[1]; - const TacsScalar dN12Crit1_dL = + TacsScalar dN12Crit1_dL = sqrt(D1 * root) * (-40.36 * D3 - 65.0 * root) / (L3 * root); - const TacsScalar dN12Crit2_dL = + TacsScalar dN12Crit2_dL = sqrt(D1 * D3) * (-7.504 * D1 * D2 - 93.6 * D32 - 4.256 * D3 * root) / (D32 * L3); @@ -2964,8 +2751,7 @@ bool TACSBladeStiffenedShellConstitutive::testCriticalShearLoadSens( " LSens = % 011.7e | LSensFD = % 011.7e | LSensRelError = % 011.7e\n", TacsRealPart(LSens), TacsRealPart(LSensFD), TacsRealPart(LSensRelError)); printf( - "----------------------------------------------------------------------" - "--" + "------------------------------------------------------------------------" "\n"); return fabs(TacsRealPart(D1SensRelError)) < tol && @@ -2977,21 +2763,20 @@ bool TACSBladeStiffenedShellConstitutive::testCriticalShearLoadSens( TacsScalar TACSBladeStiffenedShellConstitutive::bucklingEnvelope( const TacsScalar N1, const TacsScalar N1Crit, const TacsScalar N12, const TacsScalar N12Crit) { - const TacsScalar N1Term = N1 / N1Crit; - const TacsScalar N12Term = N12 / N12Crit; - const TacsScalar root = sqrt(N1Term * N1Term + 4.0 * N12Term * N12Term); + TacsScalar N1Term = N1 / N1Crit; + TacsScalar N12Term = N12 / N12Crit; + TacsScalar root = sqrt(N1Term * N1Term + 4.0 * N12Term * N12Term); return 0.5 * (N1Term + root); } -// Compute the sensitivity of the buckling failure criterion w.r.t the loads -// and critical loads +// Compute the sensitivity of the buckling failure criterion w.r.t the loads and +// critical loads TacsScalar TACSBladeStiffenedShellConstitutive::bucklingEnvelopeSens( const TacsScalar N1, const TacsScalar N1Crit, const TacsScalar N12, - const TacsScalar N12Crit, TacsScalar* const dfdN1, - TacsScalar* const dfdN1Crit, TacsScalar* const dfdN12, - TacsScalar* const dfdN12Crit) { - const TacsScalar N1Term = N1 / N1Crit; - const TacsScalar N12Term = N12 / N12Crit; + const TacsScalar N12Crit, TacsScalar* dfdN1, TacsScalar* dfdN1Crit, + TacsScalar* dfdN12, TacsScalar* dfdN12Crit) { + TacsScalar N1Term = N1 / N1Crit; + TacsScalar N12Term = N12 / N12Crit; TacsScalar root = sqrt(N1Term * N1Term + 4.0 * N12Term * N12Term); if (TacsRealPart(root) == 0.0) { root = 1e-13; @@ -3047,8 +2832,7 @@ bool TACSBladeStiffenedShellConstitutive::testBucklingEnvelopeSens( printf("testBucklingEnvelopeSens results:"); printf( - "N1 = % 011.7e, N1Crit = % 011.7e, N12 = % 011.7e, N12Crit = % " - "011.7e\n", + "N1 = % 011.7e, N1Crit = % 011.7e, N12 = % 011.7e, N12Crit = % 011.7e\n", TacsRealPart(N1), TacsRealPart(N1Crit), TacsRealPart(N12), TacsRealPart(N12Crit)); printf("f: % 011.7e\n", f); @@ -3069,552 +2853,6 @@ bool TACSBladeStiffenedShellConstitutive::testBucklingEnvelopeSens( fabs(dfdN12RelError) < tol && fabs(dfdN12CritRelError) < tol); } -TacsScalar TACSBladeStiffenedShellConstitutive::evalStiffenerColumnBuckling( - const TacsScalar stiffenerStrain[]) { - TacsScalar stiffenerStress[TACSBeamConstitutive::NUM_STRESSES]; - this->computeStiffenerStress(stiffenerStrain, stiffenerStress); - // The first component of the stiffener stress is the axial force - TacsScalar fCrit = this->computeStiffenerColumnBucklingLoad(); - return -stiffenerStress[0] / fCrit; -} - -TacsScalar -TACSBladeStiffenedShellConstitutive::computeStiffenerColumnBucklingLoad() { - TacsScalar E, G; - this->computeEffectiveModulii(this->numStiffenerPlies, this->stiffenerQMats, - this->stiffenerPlyFracs, &E, &G); - const TacsScalar Izz = this->computeStiffenerIzz(); - const TacsScalar fCrit = - this->computeColumnBucklingLoad(E, Izz, this->panelLength); - return fCrit; -} - -void TACSBladeStiffenedShellConstitutive::addStiffenerColumnBucklingLoadDVSens( - const TacsScalar scale, TacsScalar dfdx[]) { - // First compute the sensitivity of the critical load to the E, I and L - // values of the stiffener - TacsScalar E, G; - this->computeEffectiveModulii(this->numStiffenerPlies, this->stiffenerQMats, - this->stiffenerPlyFracs, &E, &G); - const TacsScalar Izz = this->computeStiffenerIzz(); - TacsScalar dFdE, dFdI, dFdL; - this->computeColumnBucklingLoadSens(E, Izz, this->panelLength, dFdE, dFdI, - dFdL); - - // Now convert those sensitivities into sensitivities with respect to the - // DVs Beam length contributions (only relevant DV is panel length) - if (this->panelLengthNum >= 0) { - dfdx[this->panelLengthLocalNum] += scale * dFdL; - } - - // Moment of inertia contributions (only relevant DVs are stiffener height - // and thickness) - TacsScalar dIzzdt, dIzzdh; - this->computeStiffenerIzzSens(dIzzdt, dIzzdh); - if (this->stiffenerHeightNum >= 0) { - dfdx[this->stiffenerHeightLocalNum] += scale * dFdI * dIzzdh; - } - if (this->stiffenerThickNum >= 0) { - dfdx[this->stiffenerThickLocalNum] += scale * dFdI * dIzzdt; - } - - // Young's modulus contributions (only relevant DVs are stiffener ply - // fractions) - for (int ii = 0; ii < this->numStiffenerPlies; ii++) { - if (this->stiffenerPlyFracNums[ii] >= 0) { - const int index = this->stiffenerPlyFracLocalNums[ii]; - - const TacsScalar* const Q = &(this->stiffenerQMats[ii * NUM_Q_ENTRIES]); - - const TacsScalar dEdx = Q[0] - (Q[1] * Q[1]) / Q[3]; - dfdx[index] += scale * dFdE * dEdx; - } - } -} - -TacsScalar TACSBladeStiffenedShellConstitutive::computeColumnBucklingLoad( - const TacsScalar E, const TacsScalar I, const TacsScalar L) { - return M_PI * M_PI * E * I / (L * L); -} - -TacsScalar TACSBladeStiffenedShellConstitutive::computeColumnBucklingLoadSens( - const TacsScalar E, const TacsScalar I, const TacsScalar L, - TacsScalar& dFdE, TacsScalar& dFdI, TacsScalar& dFdL) { - const double pi2 = M_PI * M_PI; - const TacsScalar L2inv = 1.0 / (L * L); - const TacsScalar F = pi2 * E * I * L2inv; - dFdE = pi2 * I * L2inv; - dFdI = pi2 * E * L2inv; - dFdL = -2.0 * F / L; - return F; -} - -TacsScalar -TACSBladeStiffenedShellConstitutive::evalStiffenerColumnBucklingStrainSens( - const TacsScalar stiffenerStrain[], TacsScalar sens[]) { - const int numStiff = TACSBeamConstitutive::NUM_TANGENT_STIFFNESS_ENTRIES; - TacsScalar C[numStiff]; - this->computeStiffenerStiffness(C); - const int numStress = TACSBeamConstitutive::NUM_STRESSES; - TacsScalar stiffenerStress[numStress]; - TACSBeamConstitutive::computeStress(C, stiffenerStrain, stiffenerStress); - - const TacsScalar columnBucklingLoad = - this->computeStiffenerColumnBucklingLoad(); - const TacsScalar stiffenerAxialLoad = -stiffenerStress[0]; - const TacsScalar fCrit = this->computeStiffenerColumnBucklingLoad(); - memset(sens, 0, numStress * sizeof(TacsScalar)); - - for (int ii = 0; ii < numStress; ii++) { - sens[ii] = -C[ii] / fCrit; - } - return stiffenerAxialLoad / fCrit; -} - -void TACSBladeStiffenedShellConstitutive::addStiffenerColumnBucklingDVSens( - const TacsScalar scale, const TacsScalar shellStrain[], - const TacsScalar stiffenerStrain[], const TacsScalar stiffenerAxialLoad, - const TacsScalar fCrit, TacsScalar dfdx[]) { - // F = stiffenerAxialLoad / fCrit - // Using the Quotient rule: - // dF/dx = 1/fCrit * dS11/dx - stiffenerAxialLoad/fCrit^2 * dfCrit/dx - - // First we add the 1/fCrit * dS11/dx term - // We need to consider both the direct dependence of the stiffener stress on - // the DVs, and also the dependence of the stiffener stress on the stiffener - // centroid strains, which depend on the DVS: - // psi^T dS/dx = psi^T pS/px + psi^T pS/pe * dstiffenerStrain/dx - // = psi^T pS/px + C * psi * d/dx(Te) * shellStrain - TacsScalar psi[TACSBeamConstitutive::NUM_STRESSES]; - memset(psi, 0, TACSBeamConstitutive::NUM_STRESSES * sizeof(TacsScalar)); - psi[0] = 1.0; - this->addStiffenerStressDVSens(-scale / fCrit, stiffenerStrain, psi, - &dfdx[this->stiffenerDVStartNum]); - - TacsScalar psiStress[TACSBeamConstitutive::NUM_STRESSES]; - this->computeStiffenerStress(psi, psiStress); - this->addStrainTransformProductDVsens(psiStress, shellStrain, -scale / fCrit, - dfdx); - - // Now we add the - stiffenerAxialLoad/fCrit^2 * dfCrit/dx term - this->addStiffenerColumnBucklingLoadDVSens( - -scale * stiffenerAxialLoad / (fCrit * fCrit), dfdx); -} - -void TACSBladeStiffenedShellConstitutive::computeStiffenerCripplingValues( - const TacsScalar stiffenerStrain[], TacsScalar plyFailValues[]) { - // We use the semi-empirical formula for one-edge-free crippling described - // in the DoD composite materials handbook Volume 3F (a.k.a MIL-HDBK-17-3F) - // (available at - // http://everyspec.com/MIL-HDBK/MIL-HDBK-0001-0099/MIL_HDBK_17_3F_216/) - // - // The formula is: F_crippling / F_ult = 1.63*(b/t)^-0.717 - // Where: - // F_crippling is the crippling load - // F_ult is the ultimate compressive load - // b is the flange width - // t is the flange thickness - // - // Here we compute the ultimate compressive load using the same criteria - // used to compute the stiffener material failure value but using on the - // axial strain component. This gives: SR_Crippling = SR_comp/1.63 * - // (b/t)^0.717 Where SR_Crippling is the crippling strength ratio SR_comp is - // the strength ratio computed using only the compressive component of the - // strain - - const int numPlies = this->numStiffenerPlies; - TacsScalar* const plyAngles = this->stiffenerPlyAngles; - memset(plyFailValues, 0, 2 * numPlies * sizeof(TacsScalar)); - - // We compute the crippling criteria for both the web and - // the flange of the stiffener. - TACSOrthotropicPly* ply = this->stiffenerPly; - TacsScalar zCentroid = this->computeStiffenerCentroidHeight(); - TacsScalar zFlange = -0.5 * this->stiffenerThick - zCentroid; - - TacsScalar flangeCrippleFactor = computeCripplingFactor( - 0.5 * this->flangeFraction * this->stiffenerHeight, this->stiffenerThick); - - TacsScalar plyStrain[3]; - memset(plyStrain, 0, 3 * sizeof(TacsScalar)); - plyStrain[0] = stiffenerStrain[0] + zFlange * stiffenerStrain[2]; - - if (TacsRealPart(plyStrain[0]) < 0.0) { - for (int ii = 0; ii < numPlies; ii++) { - plyFailValues[ii] = - ply->failure(plyAngles[ii], plyStrain) * flangeCrippleFactor; - } - } - - // For the stiffener web, we will use the method described by Kassapoglou to - // account for the variation in tension/compression over the height of the - // web. (See section 8.5.3 of Kassapoglou's "Design and Analysis of - // Composite Structures with Applications to Aerospace Structures"). - const TacsScalar zWebTop = - -(this->stiffenerHeight + this->stiffenerThick) - zCentroid; - const TacsScalar zWebBottom = -this->stiffenerThick - zCentroid; - const TacsScalar axStrainTop = - stiffenerStrain[0] + zWebTop * stiffenerStrain[2]; - const TacsScalar axStrainBottom = - stiffenerStrain[0] + zWebBottom * stiffenerStrain[2]; - - // There are 3 cases we have to consider: - // 1. Both strains are tensile - // 2. Both strains are compressive - // 3. One strain is tensile and the other is compressive - // - // In case 1 the strength ratio is simply zero as there can be no crippling. - // In case 2 we compute the crippling strength ratio as usual using the - // average of the two strains and the full height of the stiffener as the b - // value - // In cases 3 we compute the crippling strength ratio considering that - // only the portion of the web that is in compression can cripplie, we - // therefore use the length of that region as b and the average strain in - // that region as the strain value. - - // In theory, the point at which the strain changes sign can be treated as a - // simple support, therefore we should either treat the web as a - // one-edge-free or no-edge-free flange depending whether it is the top or - // bottom of the web that is in compression. However, this would lead to a - // discontinuity in the strength ratio, so we use the one-edge-free formula - // for both cases. - if (TacsRealPart(axStrainTop) < 0.0 || TacsRealPart(axStrainBottom) < 0.0) { - TacsScalar compressiveLength = 0.0, averageStrain = 0.0, - webCrippleFactor = 0.0; - if (TacsRealPart(axStrainTop) < 0 && TacsRealPart(axStrainBottom) < 0) { - // --- Case 2 --- - compressiveLength = this->stiffenerHeight; - averageStrain = 0.5 * (axStrainTop + axStrainBottom); - } else if (TacsRealPart(axStrainTop) < 0 && - TacsRealPart(axStrainBottom) >= 0) { - // --- Case 3 top in compression --- - compressiveLength = - this->stiffenerHeight * axStrainTop / (axStrainTop - axStrainBottom); - averageStrain = 0.5 * axStrainTop; - } else if (TacsRealPart(axStrainTop) >= 0 && - TacsRealPart(axStrainBottom) < 0) { - // --- Case 3 bottom in compression --- - compressiveLength = this->stiffenerHeight * axStrainBottom / - (axStrainBottom - axStrainTop); - averageStrain = 0.5 * axStrainBottom; - } - - webCrippleFactor = - computeCripplingFactor(compressiveLength, this->stiffenerThick); - - plyStrain[0] = averageStrain; - for (int ii = 0; ii < numPlies; ii++) { - plyFailValues[ii + numPlies] = - ply->failure(plyAngles[ii], plyStrain) * webCrippleFactor; - } - } -} -TacsScalar TACSBladeStiffenedShellConstitutive::evalStiffenerCrippling( - const TacsScalar stiffenerStrain[]) { - const int numPlies = this->numStiffenerPlies; - this->computeStiffenerCripplingValues(stiffenerStrain, - this->stiffenerPlyFailValues); - TacsScalar fail = - ksAggregation(this->stiffenerPlyFailValues, 2 * numPlies, this->ksWeight); - - return fail; -} - -TacsScalar -TACSBladeStiffenedShellConstitutive::evalStiffenerCripplingStrainSens( - const TacsScalar stiffenerStrain[], TacsScalar sens[]) { - const int numPlies = this->numStiffenerPlies; - const int numStrain = TACSBeamConstitutive::NUM_STRESSES; - TacsScalar* const plyFailValues = this->stiffenerPlyFailValues; - TacsScalar** const dFaildStrain = this->stiffenerPlyFailStrainSens; - TacsScalar* const plyAngles = this->stiffenerPlyAngles; - memset(plyFailValues, 0, 2 * numPlies * sizeof(TacsScalar)); - memset(sens, 0, numStrain * sizeof(TacsScalar)); - - for (int ii = 0; ii < 2 * numPlies; ii++) { - memset(dFaildStrain[ii], 0, numStrain * sizeof(TacsScalar)); - } - - // We compute the crippling criteria for both the web and - // the flange of the stiffener. - - // --- Flange crippling --- - TACSOrthotropicPly* ply = this->stiffenerPly; - TacsScalar zCentroid = this->computeStiffenerCentroidHeight(); - TacsScalar zFlange = -0.5 * this->stiffenerThick - zCentroid; - - TacsScalar flangeCrippleFactor = computeCripplingFactor( - 0.5 * this->flangeFraction * this->stiffenerHeight, this->stiffenerThick); - - TacsScalar plyStrain[3]; - memset(plyStrain, 0, 3 * sizeof(TacsScalar)); - plyStrain[0] = stiffenerStrain[0] + zFlange * stiffenerStrain[2]; - - if (TacsRealPart(plyStrain[0]) < 0.0) { - for (int ii = 0; ii < numPlies; ii++) { - TacsScalar plyFailStrainSens[3]; - plyFailValues[ii] = - ply->failureStrainSens(plyAngles[ii], plyStrain, plyFailStrainSens) * - flangeCrippleFactor; - // Convert the sensitivity w.r.t the ply strains to the sensitivity - // w.r.t beam strains - dFaildStrain[ii][0] = flangeCrippleFactor * plyFailStrainSens[0]; - dFaildStrain[ii][2] = - flangeCrippleFactor * zFlange * plyFailStrainSens[0]; - } - } - - // --- Web crippling --- - TacsScalar zWebTop = - -(this->stiffenerHeight + this->stiffenerThick) - zCentroid; - TacsScalar zWebBottom = -this->stiffenerThick - zCentroid; - TacsScalar axStrainTop = stiffenerStrain[0] + zWebTop * stiffenerStrain[2]; - TacsScalar axStrainBottom = - stiffenerStrain[0] + zWebBottom * stiffenerStrain[2]; - - if (TacsRealPart(axStrainTop) < 0.0 || TacsRealPart(axStrainBottom) < 0.0) { - TacsScalar compressiveLength = 0.0, averageStrain = 0.0, - webCrippleFactor = 0.0; - TacsScalar averageStrainSens[2] = {0.0, 0.0}; - TacsScalar lengthStrainSens[2] = {0.0, 0.0}; - if (TacsRealPart(axStrainTop) < 0 && TacsRealPart(axStrainBottom) < 0) { - // --- Case 2 --- - compressiveLength = this->stiffenerHeight; - averageStrain = 0.5 * (axStrainTop + axStrainBottom); - averageStrainSens[0] = 0.5; - averageStrainSens[1] = 0.5; - } else if (TacsRealPart(axStrainTop) < 0 && - TacsRealPart(axStrainBottom) >= 0) { - // --- Case 3 top in compression --- - const TacsScalar strainDiff = axStrainTop - axStrainBottom; - compressiveLength = this->stiffenerHeight * axStrainTop / strainDiff; - lengthStrainSens[0] = - -this->stiffenerHeight * axStrainBottom / (strainDiff * strainDiff); - lengthStrainSens[1] = - this->stiffenerHeight * axStrainTop / (strainDiff * strainDiff); - - averageStrain = 0.5 * axStrainTop; - averageStrainSens[0] = 0.5; - } else if (TacsRealPart(axStrainTop) >= 0 && - TacsRealPart(axStrainBottom) < 0) { - // --- Case 3 bottom in compression --- - const TacsScalar strainDiff = axStrainBottom - axStrainTop; - compressiveLength = this->stiffenerHeight * axStrainBottom / strainDiff; - lengthStrainSens[0] = - this->stiffenerHeight * axStrainBottom / (strainDiff * strainDiff); - lengthStrainSens[1] = - -this->stiffenerHeight * axStrainTop / (strainDiff * strainDiff); - - averageStrain = 0.5 * axStrainBottom; - averageStrainSens[1] = 0.5; - } - - TacsScalar dFactordL, dFactordt; - webCrippleFactor = computeCripplingFactorSens( - compressiveLength, this->stiffenerThick, dFactordL, dFactordt); - - plyStrain[0] = averageStrain; - for (int ii = 0; ii < numPlies; ii++) { - TacsScalar plyFailStrainSens[3]; - const TacsScalar strengthRatio = - ply->failureStrainSens(plyAngles[ii], plyStrain, plyFailStrainSens); - plyFailValues[ii + numPlies] = strengthRatio * webCrippleFactor; - - // Convert sensitivity w.r.t the average strain to sensitivity w.r.t - // the top and bottom strains - TacsScalar topStrainSens = - strengthRatio * dFactordL * lengthStrainSens[0] + - webCrippleFactor * plyFailStrainSens[0] * averageStrainSens[0]; - - TacsScalar bottomStrainSens = - strengthRatio * dFactordL * lengthStrainSens[1] + - webCrippleFactor * plyFailStrainSens[0] * averageStrainSens[1]; - - // Convert the sensitivity w.r.t the top and bottom strains to the - // sensitivity w.r.t beam strains - dFaildStrain[ii + numPlies][0] = topStrainSens + bottomStrainSens; - dFaildStrain[ii + numPlies][2] = - zWebTop * topStrainSens + zWebBottom * bottomStrainSens; - } - } - - TacsScalar fail = - ksAggregationSensProduct(plyFailValues, 2 * numPlies, numStrain, - this->ksWeight, dFaildStrain, sens); - - return fail; -} - -void TACSBladeStiffenedShellConstitutive::addStiffenerCripplingDVSens( - const TacsScalar scale, const TacsScalar stiffenerStrain[], - TacsScalar dfdx[]) { - TACSOrthotropicPly* ply = this->stiffenerPly; - TacsScalar* plyFailValues = this->stiffenerPlyFailValues; - const TacsScalar* plyAngles = this->stiffenerPlyAngles; - TacsScalar* dKSdFail = this->stiffenerPlyFailSens; - const int hNum = this->stiffenerHeightLocalNum - this->stiffenerDVStartNum; - const int tNum = this->stiffenerThickLocalNum - this->stiffenerDVStartNum; - const int numPlies = this->numStiffenerPlies; - const int numStrain = TACSBeamConstitutive::NUM_STRESSES; - const bool computeThicknessSens = (this->stiffenerThickNum >= 0); - const bool computeHeightSens = (this->stiffenerHeightNum >= 0); - - this->computeStiffenerCripplingValues(stiffenerStrain, plyFailValues); - - const TacsScalar fail = - ksAggregationSens(plyFailValues, 2 * numPlies, this->ksWeight, dKSdFail); - - // --- Flange crippling --- - TacsScalar dzCentroiddh, dzCentroiddt; - TacsScalar zCentroid = this->computeStiffenerCentroidHeight(); - this->computeStiffenerCentroidHeightSens(dzCentroiddt, dzCentroiddh); - TacsScalar zFlange = -0.5 * this->stiffenerThick - zCentroid; - const TacsScalar dzFlangedt = -0.5 - dzCentroiddt; - const TacsScalar dzFlangedh = -dzCentroiddh; - - TacsScalar dFactordh, dFactordt; - TacsScalar flangeCrippleFactor = computeCripplingFactorSens( - 0.5 * this->flangeFraction * this->stiffenerHeight, this->stiffenerThick, - dFactordh, dFactordt); - dFactordh *= 0.5 * this->flangeFraction; - - TacsScalar plyStrain[3]; - memset(plyStrain, 0, 3 * sizeof(TacsScalar)); - plyStrain[0] = stiffenerStrain[0] + zFlange * stiffenerStrain[2]; - if (TacsRealPart(plyStrain[0]) < 0.0) { - for (int ii = 0; ii < numPlies; ii++) { - TacsScalar plyFailStrainSens[3]; - const TacsScalar strengthRatio = - ply->failureStrainSens(plyAngles[ii], plyStrain, plyFailStrainSens); - if (computeThicknessSens) { - dfdx[tNum] += dKSdFail[ii] * scale * - (strengthRatio * dFactordt + - plyFailStrainSens[0] * dzFlangedt * stiffenerStrain[2] * - flangeCrippleFactor); - } - if (computeHeightSens) { - dfdx[hNum] += dKSdFail[ii] * scale * - (strengthRatio * dFactordh + - plyFailStrainSens[0] * dzFlangedh * stiffenerStrain[2] * - flangeCrippleFactor); - } - } - } - - // --- Flange crippling --- - const TacsScalar zWebTop = - -(this->stiffenerHeight + this->stiffenerThick) - zCentroid; - const TacsScalar dzWebTopdh = -1.0 - dzCentroiddh; - const TacsScalar dzWebTopdt = -1.0 - dzCentroiddt; - - const TacsScalar zWebBottom = -this->stiffenerThick - zCentroid; - const TacsScalar dzWebBottomdt = -1.0 - dzCentroiddt; - const TacsScalar dzWebBottomdh = -dzCentroiddh; - - const TacsScalar axStrainTop = - stiffenerStrain[0] + zWebTop * stiffenerStrain[2]; - const TacsScalar dStrainTopdh = stiffenerStrain[2] * dzWebTopdh; - const TacsScalar dStrainTopdt = stiffenerStrain[2] * dzWebTopdt; - - const TacsScalar axStrainBottom = - stiffenerStrain[0] + zWebBottom * stiffenerStrain[2]; - const TacsScalar dStrainBottomdh = stiffenerStrain[2] * dzWebBottomdh; - const TacsScalar dStrainBottomdt = stiffenerStrain[2] * dzWebBottomdt; - - if (TacsRealPart(axStrainTop) < 0.0 || TacsRealPart(axStrainBottom) < 0.0) { - TacsScalar compressiveLength = 0.0, averageStrain = 0.0, - webCrippleFactor = 0.0; - TacsScalar dAverageStraindh = 0.0; - TacsScalar dAverageStraindt = 0.0; - TacsScalar dLdh = 0.0; - TacsScalar dLdt = 0.0; - if (TacsRealPart(axStrainTop) < 0 && TacsRealPart(axStrainBottom) < 0) { - // --- Case 2 --- - compressiveLength = this->stiffenerHeight; - dLdh = 1.0; - - averageStrain = 0.5 * (axStrainTop + axStrainBottom); - dAverageStraindh = 0.5 * (dStrainTopdh + dStrainBottomdh); - dAverageStraindt = 0.5 * (dStrainTopdt + dStrainBottomdt); - } else if (TacsRealPart(axStrainTop) < 0 && - TacsRealPart(axStrainBottom) >= 0) { - // --- Case 3 top in compression --- - const TacsScalar strainDiff = axStrainTop - axStrainBottom; - compressiveLength = this->stiffenerHeight * axStrainTop / strainDiff; - - const TacsScalar dLdstrainTop = - -this->stiffenerHeight * axStrainBottom / (strainDiff * strainDiff); - const TacsScalar dLdstrainBottom = - this->stiffenerHeight * axStrainTop / (strainDiff * strainDiff); - - // dLdx = pLpx + pLpet * detdx + pLpeb * debdx - dLdh = axStrainTop / strainDiff + dLdstrainTop * dStrainTopdh + - dLdstrainBottom * dStrainBottomdh; - dLdt = dLdstrainTop * dStrainTopdt + dLdstrainBottom * dStrainBottomdt; - - averageStrain = 0.5 * axStrainTop; - dAverageStraindh = 0.5 * dStrainTopdh; - dAverageStraindt = 0.5 * dStrainTopdt; - } else if (TacsRealPart(axStrainTop) >= 0 && - TacsRealPart(axStrainBottom) < 0) { - // --- Case 3 bottom in compression --- - const TacsScalar strainDiff = axStrainBottom - axStrainTop; - - compressiveLength = this->stiffenerHeight * axStrainBottom / strainDiff; - - const TacsScalar dLdstrainTop = - this->stiffenerHeight * axStrainBottom / (strainDiff * strainDiff); - const TacsScalar dLdstrainBottom = - -this->stiffenerHeight * axStrainTop / (strainDiff * strainDiff); - - // dLdx = pLpx + pLpet * detdx + pLpeb * debdx - dLdh = axStrainBottom / strainDiff + dLdstrainTop * dStrainTopdh + - dLdstrainBottom * dStrainBottomdh; - dLdt = dLdstrainTop * dStrainTopdt + dLdstrainBottom * dStrainBottomdt; - - averageStrain = 0.5 * axStrainBottom; - dAverageStraindh = 0.5 * dStrainBottomdh; - dAverageStraindt = 0.5 * dStrainBottomdt; - } - - TacsScalar dFactordL, dFactordt; - webCrippleFactor = computeCripplingFactorSens( - compressiveLength, this->stiffenerThick, dFactordL, dFactordt); - const TacsScalar dFactordh = dFactordL * dLdh; - dFactordt += dFactordL * dLdt; - - plyStrain[0] = averageStrain; - for (int ii = 0; ii < numPlies; ii++) { - TacsScalar plyFailStrainSens[3]; - const TacsScalar strengthRatio = - ply->failureStrainSens(plyAngles[ii], plyStrain, plyFailStrainSens); - if (computeThicknessSens) { - const TacsScalar dFaildt = - strengthRatio * dFactordt + - webCrippleFactor * plyFailStrainSens[0] * dAverageStraindt; - - dfdx[tNum] += dKSdFail[ii + numPlies] * scale * dFaildt; - } - if (computeHeightSens) { - const TacsScalar dFaildh = - strengthRatio * dFactordh + - webCrippleFactor * plyFailStrainSens[0] * dAverageStraindh; - - dfdx[hNum] += dKSdFail[ii + numPlies] * scale * dFaildh; - } - } - } -} - -TacsScalar TACSBladeStiffenedShellConstitutive::computeCripplingFactorSens( - const TacsScalar b, const TacsScalar t, TacsScalar& dkdb, - TacsScalar& dkdt) { - const TacsScalar factor = pow(b / t, 0.717) / 1.63; - dkdb = 0.717 * factor / b; - dkdt = -0.717 * factor / t; - return factor; -} - // ============================================================================== // Utility functions // ============================================================================== diff --git a/src/constitutive/TACSBladeStiffenedShellConstitutive.h b/src/constitutive/TACSBladeStiffenedShellConstitutive.h index aa5317cd9..0e1d93148 100644 --- a/src/constitutive/TACSBladeStiffenedShellConstitutive.h +++ b/src/constitutive/TACSBladeStiffenedShellConstitutive.h @@ -131,14 +131,6 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { * @param _stiffenerPlyFracNums Stiffener ply fraction design variable numbers * @param _flangeFraction Stiffener base width as a fraction of the stiffener * height - * @param _includeGlobalBuckling Whether to include global panel buckling in - * the aggregated failure criteria - * @param _includeLocalBuckling Whether to include local inter-stringer - * buckling in the aggregated failure criteria - * @param _includeStiffenerColumnBuckling Whether to include stiffener column - * buckling in the aggregated failure criteria - * @param _includeStiffenerCrippling Whether to include stiffener crippling in - * the aggregated failure criteria */ TACSBladeStiffenedShellConstitutive( TACSOrthotropicPly* _panelPly, TACSOrthotropicPly* _stiffenerPly, @@ -150,69 +142,13 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { int _stiffenerHeightNum, TacsScalar _stiffenerThick, int _stiffenerThickNum, int _numStiffenerPlies, TacsScalar _stiffenerPlyAngles[], TacsScalar _stiffenerPlyFracs[], - int _stiffenerPlyFracNums[], TacsScalar _flangeFraction = 1.0, - bool _includeGlobalBuckling = true, bool _includeLocalBuckling = true, - bool _includeStiffenerColumnBuckling = true, - bool _includeStiffenerCrippling = true); + int _stiffenerPlyFracNums[], TacsScalar _flangeFraction = 1.0); ~TACSBladeStiffenedShellConstitutive(); // ============================================================================== // Set non-default values // ============================================================================== - - /** - * @brief Enable or disable the material failure mode in the aggregated - * failure criteria - * - * @param includeMaterialFailure Whether to include material failure - */ - void setIncludeMaterialFailure(bool _includMaterialFailure) { - this->includeMaterialFailure = _includMaterialFailure; - } - - /** - * @brief Enable or disable the global panel buckling mode in the aggregated - * failure criteria - * - * @param includeGlobalBuckling Whether to include global panel buckling - */ - void setIncludeGlobalBuckling(bool _includeGlobalBuckling) { - this->includeGlobalBuckling = _includeGlobalBuckling; - } - - /** - * @brief Enable or disable the local inter-stringer buckling mode in the - * aggregated failure criteria - * - * @param _includeLocalBuckling Whether to include local inter-stringer - * buckling - */ - void setIncludeLocalBuckling(bool _includeLocalBuckling) { - this->includeLocalBuckling = _includeLocalBuckling; - } - - /** - * @brief Enable or disable the stiffener column buckling mode in the - * aggregated failure criteria - * - * @param _includeStiffenerColumnBuckling Whether to include stiffener column - * buckling - */ - void setIncludeStiffenerColumnBuckling(bool _includeStiffenerColumnBuckling) { - this->includeStiffenerColumnBuckling = _includeStiffenerColumnBuckling; - } - - /** - * @brief Enable or disable the stiffener crippling mode in the aggregated - * failure criteria - * - * @param _includeStiffenerCrippling Whether to include stiffener crippling - */ - void setIncludeStiffenerCrippling(bool _includeStiffenerCrippling) { - this->includeStiffenerCrippling = _includeStiffenerCrippling; - } - /** * @brief Set the Stiffener Pitch DV Bounds * @@ -717,7 +653,7 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { static void computeEffectiveModulii(const int numPlies, const TacsScalar QMats[], const TacsScalar plyFracs[], - TacsScalar* const E, TacsScalar* const G); + TacsScalar* E, TacsScalar* G); /** * @brief Compute the failure criterion for the stiffener in a given strain @@ -851,14 +787,6 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { // Buckling functions // ============================================================================== - /** - * @brief Compute the strength ratio for the global buckling of the panel - * - * @param e Shell strains - * @return TacsScalar Strength ratio - */ - TacsScalar evalGlobalPanelBuckling(const TacsScalar e[]); - /** * @brief Compute the panel + stiffener stiffness values used to compute the * global buckling loads @@ -867,37 +795,8 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { * @output D2 2-direction bending stiffness * @output D3 Twisting stiffness */ - void computeCriticalGlobalBucklingStiffness(TacsScalar* const D1, - TacsScalar* const D2, - TacsScalar* const D3); - - /** - * @brief Compute the sensitivity of the global buckling strength ratio w.r.t - * the shell strains - * - * @param e Shell strains - * @param sens Sensitivity of the output w.r.t the shell strains - * @return TacsScalar Strength Ratio - */ - TacsScalar evalGlobalPanelBucklingStrainSens(const TacsScalar e[], - TacsScalar sens[]); - - /** - * @brief Add the derivative of the global panel buckling strength ratio w.r.t - * the design variables - - @param elemIndex The local element index (not used) - @param scale Value by which to scale the derivatives - @param pt The parametric point (not used) - @param X The physical node location (not used) - @param strain The shell strains - @param dvLen The length of the design vector (not used) - @param dfdx The DV sensitivity array to add to - */ - void addGlobalPanelBucklingDVSens(int elemIndex, TacsScalar scale, - const double pt[], const TacsScalar X[], - const TacsScalar strain[], int dvLen, - TacsScalar dfdx[]); + void computeCriticalGlobalBucklingStiffness(TacsScalar* D1, TacsScalar* D2, + TacsScalar* D3); /** * @brief Compute the sensitivities of the panel + stiffener stiffness values @@ -917,9 +816,8 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { */ void computeCriticalGlobalBucklingStiffnessSens( const TacsScalar dfdD1, const TacsScalar dfdD2, const TacsScalar dfdD3, - TacsScalar* const spSens, TacsScalar* const tpSens, - TacsScalar* const hsSens, TacsScalar* const tsSens, - TacsScalar QstiffSens[], TacsScalar QpanelSens[]); + TacsScalar* spSens, TacsScalar* tpSens, TacsScalar* hsSens, + TacsScalar* tsSens, TacsScalar QstiffSens[], TacsScalar QpanelSens[]); /** * @brief Compute the critical axial load for the global buckling of the @@ -935,15 +833,6 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { return M_PI * M_PI * D1 / (L * L); } - /** - * @brief Compute the strength ratio for the local buckling of the panel skin - * between stiffeners - * - * @param e Shell strains - * @return TacsScalar Strength ratio - */ - TacsScalar evalLocalPanelBuckling(const TacsScalar e[]); - /** * @brief Compute the critical axial load for local buckling of the panel * @@ -962,33 +851,6 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { return 2.0 * M_PI * M_PI / (L * L) * (sqrt(D11 * D22) + D12 + 2.0 * D66); } - /** - * @brief Compute the sensitivity of the local panel buckling strength ratio - * - * @param e Shell strains - * @param sens Sensitivity of the output w.r.t the shell strains - * @return TacsScalar Strength Ratio - */ - TacsScalar evalLocalPanelBucklingStrainSens(const TacsScalar e[], - TacsScalar sens[]); - - /** - * @brief Add the derivative of the local panel buckling strength ratio w.r.t - * the design variables - - @param elemIndex The local element index (not used) - @param scale Value by which to scale the derivatives - @param pt The parametric point (not used) - @param X The physical node location (not used) - @param strain The shell strains - @param dvLen The length of the design vector (not used) - @param dfdx The DV sensitivity array to add to - */ - void addLocalPanelBucklingDVSens(int elemIndex, TacsScalar scale, - const double pt[], const TacsScalar X[], - const TacsScalar strain[], int dvLen, - TacsScalar dfdx[]); - /** * @brief Compute the sensitivity of the critical axial load for local * buckling of the panel @@ -1007,9 +869,9 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { */ static TacsScalar computeCriticalLocalAxialLoadSens( const TacsScalar D11, const TacsScalar D22, const TacsScalar D12, - const TacsScalar D66, const TacsScalar L, TacsScalar* const D11Sens, - TacsScalar* const D22Sens, TacsScalar* const D12Sens, - TacsScalar* const D66Sens, TacsScalar* const LSens); + const TacsScalar D66, const TacsScalar L, TacsScalar* D11Sens, + TacsScalar* D22Sens, TacsScalar* D12Sens, TacsScalar* D66Sens, + TacsScalar* LSens); /** * @brief Compute the critical shear load for either local or global buckling @@ -1080,8 +942,8 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { */ static TacsScalar computeCriticalShearLoadSens( const TacsScalar D1, const TacsScalar D2, const TacsScalar D3, - const TacsScalar L, TacsScalar* const D1Sens, TacsScalar* const D2Sens, - TacsScalar* const D3Sens, TacsScalar* const LSens); + const TacsScalar L, TacsScalar* D1Sens, TacsScalar* D2Sens, + TacsScalar* D3Sens, TacsScalar* LSens); /** * @brief Compute the buckling failure criterion @@ -1116,141 +978,8 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { */ static TacsScalar bucklingEnvelopeSens( const TacsScalar N1, const TacsScalar N1Crit, const TacsScalar N12, - const TacsScalar N12Crit, TacsScalar* const N1Sens, - TacsScalar* const N1CritSens, TacsScalar* const N12Sens, - TacsScalar* const N12CritSens); - - /** - * @brief Compute the strength ratio for the stiffener column buckling failure - * mode - * - * @param stiffenerStrain Stiffener centroid strains - * @return TacsScalar Strength ratio - */ - TacsScalar evalStiffenerColumnBuckling(const TacsScalar stiffenerStrain[]); - - /** - * @brief Compute the critical buckling load of the panel stiffeners - * - * @return TacsScalar - */ - TacsScalar computeStiffenerColumnBucklingLoad(); - - /** - * @brief Add the sensitivity of the stiffener column buckling failure with - * respect to the DVs - * - * @param scale Scaling factor to apply to the sensitivities (e.g the - * sensitivity of the final output to the stiffener column buckling load) - * @param dfdx Array of DV sensitivities - */ - void addStiffenerColumnBucklingLoadDVSens(const TacsScalar scale, - TacsScalar dfdx[]); - - /** - * @brief Compute the critical column buckling load of a beam - * - * @param E Elastic modulus - * @param I Bending moment of inertia - * @param L Beam length - * @return TacsScalar Critical column buckling load - */ - static TacsScalar computeColumnBucklingLoad(const TacsScalar E, - const TacsScalar I, - const TacsScalar L); - - /** - * @brief Compute the sensitivity of the critical column buckling load of a - * beam w.r.t the it's stiffness and length - * - * @param E Elastic modulus - * @param I Bending moment of inertia - * @param L Beam length - * @param dFdE Sensitivity of the critical load w.r.t the elastic modulus - * @param dFdI Sensitivity of the critical load w.r.t the bending moment of - * inertia - * @param dFdL Sensitivity of the critical load w.r.t the beam length - * @return TacsScalar Critical column buckling load - */ - static TacsScalar computeColumnBucklingLoadSens( - const TacsScalar E, const TacsScalar I, const TacsScalar L, - TacsScalar& dFdE, TacsScalar& dFdI, TacsScalar& dFdL); - - /** - * @brief Compute the sensitivity of the stiffener column buckling failure - * criteria with respect to the beam strains - * - * @param stiffenerStrain Stiffener centroid beam strains - * @param sens Sensitivity of the column failure criteria w.r.t the beam - * strains - * @return TacsScalar The column buckling failure criteria value - */ - TacsScalar evalStiffenerColumnBucklingStrainSens( - const TacsScalar stiffenerStrain[], TacsScalar sens[]); - - /** - * @brief Add the sensitivity of the stiffener column buckling failure - * criteria with respect to the design variables to an existing strain - * sensitivity vector - * - * @param scale - * @param shellStrain - * @param stiffenerAxialLoad - * @param fCrit - * @param dfdx - */ - void addStiffenerColumnBucklingDVSens(const TacsScalar scale, - const TacsScalar shellStrain[], - const TacsScalar stiffenerStrain[], - const TacsScalar stiffenerAxialLoad, - const TacsScalar fCrit, - TacsScalar dfdx[]); - - /** - * @brief Compute the empirical crippling factor for a stiffener flange - * - * @param b Flange length - * @param t Flange thickness - * @return TacsScalar Crippling factor by which the compressive strength - * should be scaled down or, conversely, the strength ratio should be scaled - * up - */ - static TacsScalar computeCripplingFactor(const TacsScalar b, - const TacsScalar t) { - return pow(b / t, 0.717) / 1.63; - } - - /** - * @brief Compute the sensitivity of the flange crippling factor w.r.t the - * flange dimensions - * - * @param b Flange length - * @param t Flange thickness - * @return dkdb Sensitivity of the crippling factor w.r.t the flange length - * @return dkdt Sensitivity of the crippling factor w.r.t the flange thickness - * @return TacsScalar The crippling factor - */ - static TacsScalar computeCripplingFactorSens(const TacsScalar b, - const TacsScalar t, - TacsScalar& dkdb, - TacsScalar& dkdt); - - /** - * @brief Compute the strength ratio with respect to stiffener crippling - * - * @param stiffenerStrain Stiffener centroid beam strains - * @return TacsScalar Strength ratio - */ - TacsScalar evalStiffenerCrippling(const TacsScalar stiffenerStrain[]); - void computeStiffenerCripplingValues(const TacsScalar stiffenerStrain[], - TacsScalar plyFailValues[]); - - TacsScalar evalStiffenerCripplingStrainSens( - const TacsScalar stiffenerStrain[], TacsScalar sens[]); - - void addStiffenerCripplingDVSens(const TacsScalar scale, - const TacsScalar stiffenerStrain[], - TacsScalar dfdx[]); + const TacsScalar N12Crit, TacsScalar* N1Sens, TacsScalar* N1CritSens, + TacsScalar* N12Sens, TacsScalar* N12CritSens); // ============================================================================== // Attributes @@ -1264,16 +993,6 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { int numGeneralDV; ///< Number of general DVs int numPanelDV; ///< Number of panel DVs int numStiffenerDV; ///< Number of stiffener DVs - bool includeMaterialFailure = - true; ///< Whether to consider panel material failure - bool includeGlobalBuckling = - true; ///< Whether to consider global panel buckling failure - bool includeLocalBuckling = true; ///< Whether to consider local buckling - ///< failure of the skin between stiffeners - bool includeStiffenerColumnBuckling = - true; ///< Whether to consider stiffener column buckling failure - bool includeStiffenerCrippling = - true; ///< Whether to consider stiffener crippling failure // --- Material properties --- TACSOrthotropicPly* panelPly; ///< Orthotropic ply for the panel @@ -1347,18 +1066,9 @@ class TACSBladeStiffenedShellConstitutive : public TACSShellConstitutive { TacsScalar* panelPlyFailSens; TacsScalar* stiffenerPlyFailSens; - static const char* const constName; ///< Constitutive model name + static const char* constName; ///< Constitutive model name static const int NUM_Q_ENTRIES = 6; ///< Number of entries in the Q matrix static const int NUM_ABAR_ENTRIES = - 3; ///< Number of entries in the ABar matrix - static const int NUM_FAILURES = - 6; ///< Number of failure modes, we have: - ///< 1. Panel material failure - ///< 2. Panel material failure - ///< 3. Local panel buckling (between stiffeners) - ///< 4. Global panel buckling - ///< 5. Stiffener column buckling - ///< 6. Stiffener crippling - static constexpr TacsScalar DUMMY_FAIL_VALUE = - -1e200; ///< Dummy failure value used for failure modes that are disabled + 3; ///< Number of entries in the ABar matrix + static const int NUM_FAILURES = 4; ///< Number of failure modes }; diff --git a/tacs/constitutive.pyx b/tacs/constitutive.pyx index c2cc91260..fb332bfd3 100644 --- a/tacs/constitutive.pyx +++ b/tacs/constitutive.pyx @@ -910,11 +910,7 @@ cdef class BladeStiffenedShellConstitutive(ShellConstitutive): np.ndarray[int, ndim=1, mode='c'] panelPlyFracNums = None, int stiffenerHeightNum = -1, int stiffenerThickNum = -1, - np.ndarray[int, ndim=1, mode='c'] stiffenerPlyFracNums = None, - bool includeGlobalBuckling=True, - bool includeLocalBuckling=True, - bool includeStiffenerColumnBuckling=True, - bool includeStiffenerCrippling=True + np.ndarray[int, ndim=1, mode='c'] stiffenerPlyFracNums = None ): numPanelPlies = len(panelPlyAngles) @@ -963,36 +959,11 @@ cdef class BladeStiffenedShellConstitutive(ShellConstitutive): stiffenerPlyAngles.data, stiffenerPlyFracs.data, stiffenerPlyFracNums.data, - flangeFraction, - includeGlobalBuckling, - includeLocalBuckling, - includeStiffenerColumnBuckling, - includeStiffenerCrippling + flangeFraction ) self.ptr = self.cptr = self.blade_ptr self.ptr.incref() - def setFailureModes( - self, - includeMaterialFailure=None, - includeGlobalBuckling=None, - includeLocalBuckling=None, - includeStiffenerColumnBuckling=None, - includeStiffenerCrippling=None, - ): - - if self.blade_ptr: - if includeMaterialFailure is not None: - self.blade_ptr.setIncludeMaterialFailure(includeMaterialFailure) - if includeGlobalBuckling is not None: - self.blade_ptr.setIncludeGlobalBuckling(includeGlobalBuckling) - if includeLocalBuckling is not None: - self.blade_ptr.setIncludeLocalBuckling(includeLocalBuckling) - if includeStiffenerColumnBuckling is not None: - self.blade_ptr.setIncludeStiffenerColumnBuckling(includeStiffenerColumnBuckling) - if includeStiffenerCrippling is not None: - self.blade_ptr.setIncludeStiffenerCrippling(includeStiffenerCrippling) - def setKSWeight(self, double ksWeight): """ Update the ks weight used for aggregating the different failure modes diff --git a/tacs/cpp_headers/constitutive.pxd b/tacs/cpp_headers/constitutive.pxd index 326cf123b..d1d3ceb71 100644 --- a/tacs/cpp_headers/constitutive.pxd +++ b/tacs/cpp_headers/constitutive.pxd @@ -1,5 +1,7 @@ # This file is part of TACS: The Toolkit for the Analysis of Composite # Structures, a parallel finite-element code for structural and +# multidisciplinary de# This file is part of TACS: The Toolkit for the Analysis of Composite +# Structures, a parallel finite-element code for structural and # multidisciplinary design optimization. # # Copyright (C) 2014 Georgia Tech Research Corporation @@ -152,19 +154,46 @@ cdef extern from "TACSBladeStiffenedShellConstitutive.h": TacsScalar[], # stiffenerPlyAngles TacsScalar[], # stiffenerPlyFracs int[], # stiffenerPlyFracNums - TacsScalar, # flangeFraction - bool, # includeGlobalBuckling, - bool, # includeLocalBuckling, - bool, # includeStiffenerColumnBuckling, - bool, # includeStiffenerCrippling + TacsScalar # flangeFraction + ) + int getNumPanelPlies() + int getNumStiffenerPlies() + void setKSWeight(double ksWeight) + void setStiffenerPitchBounds(TacsScalar lowerBound, TacsScalar upperBound) + void setStiffenerHeightBounds(TacsScalar lowerBound, TacsScalar upperBound) + void setStiffenerThicknessBounds(TacsScalar lowerBound, TacsScalar upperBound) + void setPanelThicknessBounds(TacsScalar lowerBound, TacsScalar upperBound) + void setStiffenerPlyFractionBounds(TacsScalar[] lowerBound, TacsScalar[] upperBound) + void setPanelPlyFractionBounds(TacsScalar[] lowerBound, TacsScalar[] upperBound) + +cdef extern from "TACSBladeStiffenedShellConstitutive.h": + cdef cppclass TACSBladeStiffenedShellConstitutive(TACSShellConstitutive): + TACSBladeStiffenedShellConstitutive( + TACSOrthotropicPly*, # panelPly + TACSOrthotropicPly*, # stiffenerPly + TacsScalar, # kcorr + TacsScalar, # panelLength + int, # panelLengthNum + TacsScalar, # stiffenerPitch + int, # stiffenerPitchNum + TacsScalar, # panelThick + int, # panelThickNum + int, # numPanelPlies + TacsScalar[], # panelPlyAngles + TacsScalar[], # panelPlyFracs + int[], # panelPlyFracNums + TacsScalar, # stiffenerHeight + int, # stiffenerHeightNum + TacsScalar, # stiffenerThick + int, # stiffenerThickNum + int, # numStiffenerPlies + TacsScalar[], # stiffenerPlyAngles + TacsScalar[], # stiffenerPlyFracs + int[], # stiffenerPlyFracNums + TacsScalar # flangeFraction ) int getNumPanelPlies() int getNumStiffenerPlies() - void setIncludeMaterialFailure(bool _includeMaterialFailure) - void setIncludeGlobalBuckling(bool _includeGlobalBuckling) - void setIncludeLocalBuckling(bool _includeLocalBuckling) - void setIncludeStiffenerColumnBuckling(bool _includeStiffenerColumnBuckling) - void setIncludeStiffenerCrippling(bool _includeStiffenerCrippling) void setKSWeight(double ksWeight) void setStiffenerPitchBounds(TacsScalar lowerBound, TacsScalar upperBound) void setStiffenerHeightBounds(TacsScalar lowerBound, TacsScalar upperBound) diff --git a/tests/constitutive_tests/test_blade_sitffened_shell_constitutive.py b/tests/constitutive_tests/test_blade_sitffened_shell_constitutive.py index bd335f006..4f1d11b93 100644 --- a/tests/constitutive_tests/test_blade_sitffened_shell_constitutive.py +++ b/tests/constitutive_tests/test_blade_sitffened_shell_constitutive.py @@ -34,15 +34,15 @@ def setUp(self): self.x = np.ones(3, dtype=self.dtype) self.pt = np.zeros(3) - self.panelLength = 2.1 + self.panelLength = 2.0 self.panelLengthNum = 0 - self.stiffenerPitch = 0.178 + self.stiffenerPitch = 0.2 self.stiffenerPitchNum = 1 - self.stiffenerHeight = 0.314 + self.stiffenerHeight = 0.075 self.stiffenerHeightNum = 2 - self.stiffenerThickness = 1.23e-2 + self.stiffenerThickness = 1e-2 self.stiffenerThicknessNum = 3 - self.panelThickness = 1.586e-2 + self.panelThickness = 1.5e-2 self.panelThicknessNum = 4 self.numPanelPlies = 3 @@ -69,13 +69,11 @@ def setUp(self): [ self.panelLength, self.stiffenerPitch, - self.panelThickness, - ] - + list(self.panelPlyFracs) - + [ self.stiffenerHeight, self.stiffenerThickness, + self.panelThickness, ] + + list(self.panelPlyFracs) + list(self.stiffenerPlyFracs) ) self.dvs = np.array(self.dvs, dtype=self.dtype) @@ -138,16 +136,6 @@ def setUp(self): self.ply_list = [iso_ply, ortho_ply] - # These are the individual failure modes that can be enabled/disabled in the model. We will run the failure sensitivity tests with each enabled individually and then again with all enabled. - self.failure_modes = [ - "MaterialFailure", - "GlobalBuckling", - "LocalBuckling", - "StiffenerColumnBuckling", - "StiffenerCrippling", - ] - self.failure_modes_to_test = self.failure_modes + ["all"] - # Seed random number generator in tacs for consistent test results elements.SeedRandomGenerator(0) @@ -172,11 +160,11 @@ def get_con(self, ply): self.panelPlyFracNums, self.stiffenerHeightNum, self.stiffenerThicknessNum, - self.stiffenerPlyFracNums, + self.stiffenerPlyFracNums ) # Set the KS weight really low so that all failure modes make a # significant contribution to the failure function derivatives - con.setKSWeight(1.0) + con.setKSWeight(0.1) return con def test_constitutive_density(self): @@ -269,64 +257,43 @@ def test_constitutive_stress(self): # ) # self.assertFalse(fail) - def test_constitutive_failure_dv_sens(self): + def test_constitutive_failure(self): # Test failure dv sensitivity for ply in self.ply_list: with self.subTest(ply=ply): - for enabled_failure_mode in self.failure_modes_to_test: - with self.subTest(failure_mode=enabled_failure_mode): - includeFailureModes = {} - for failure_mode in self.failure_modes: - includeFailureModes[f"include{failure_mode}"] = ( - failure_mode == enabled_failure_mode - or enabled_failure_mode == "all" - ) - con = self.get_con(ply) - con.setFailureModes(**includeFailureModes) - fail = False - for _ in range(self.numFailureTests): - fail = fail or constitutive.TestConstitutiveFailure( - con, - self.elem_index, - self.pt, - self.x, - self.dvs, - self.dh, - self.print_level, - self.atol, - self.rtol, - ) - self.assertFalse(fail) + con = self.get_con(ply) + fail = False + for _ in range(self.numFailureTests): + fail = fail or constitutive.TestConstitutiveFailure( + con, + self.elem_index, + self.pt, + self.x, + self.dvs, + self.dh, + self.print_level, + self.atol, + self.rtol, + ) + self.assertFalse(fail) def test_constitutive_failure_strain_sens(self): for ply in self.ply_list: with self.subTest(ply=ply): - for enabled_failure_mode in self.failure_modes_to_test: - with self.subTest(failure_mode=enabled_failure_mode): - includeFailureModes = {} - for failure_mode in self.failure_modes: - includeFailureModes[f"include{failure_mode}"] = ( - failure_mode == enabled_failure_mode - or enabled_failure_mode == "all" - ) - con = self.get_con(ply) - con.setFailureModes(**includeFailureModes) - fail = False - for _ in range(self.numFailureTests): - fail = ( - fail - or constitutive.TestConstitutiveFailureStrainSens( - con, - self.elem_index, - self.pt, - self.x, - self.dh, - self.print_level, - self.rtol, - self.atol, - ) - ) - self.assertFalse(fail) + con = self.get_con(ply) + fail = False + for _ in range(self.numFailureTests): + fail = fail or constitutive.TestConstitutiveFailureStrainSens( + con, + self.elem_index, + self.pt, + self.x, + self.dh, + self.print_level, + self.rtol, + self.atol, + ) + self.assertFalse(fail) if __name__ == "__main__": diff --git a/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py b/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py index 42e49c506..9d3239a1a 100644 --- a/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py +++ b/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py @@ -42,7 +42,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest): "point_load_cgy": 0.500000000000004, "point_load_cgz": -0.0035714285714285718, "point_load_compliance": 71.87161795633577, - "point_load_ks_failure": 1.55459927291223, + "point_load_ks_vmfailure": 0.3304517634314118, "point_load_mass": 17.5, "pressure_Ixx": 1.4589304315475928, "pressure_Ixy": 3.907985046680551e-14, @@ -54,7 +54,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest): "pressure_cgy": 0.500000000000004, "pressure_cgz": -0.0035714285714285718, "pressure_compliance": 377.11604579572935, - "pressure_ks_failure": 2.071155983783409, + "pressure_ks_vmfailure": 0.9085085771277308, "pressure_mass": 17.5, "gravity_Ixx": 1.4589304315475928, "gravity_Ixy": 3.907985046680551e-14, @@ -66,7 +66,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest): "gravity_cgy": 0.500000000000004, "gravity_cgz": -0.0035714285714285718, "gravity_compliance": 11.114479357783475, - "gravity_ks_failure": 0.36091429055815866, + "gravity_ks_vmfailure": 0.0908244150089233, "gravity_mass": 17.5, "modal_eigsm.0": 728895.1077101853, "modal_eigsm.1": 1591857.6791554866,