diff --git a/src/constitutive/TACSMaterialProperties.cpp b/src/constitutive/TACSMaterialProperties.cpp index 74a8c4399..504b1d787 100644 --- a/src/constitutive/TACSMaterialProperties.cpp +++ b/src/constitutive/TACSMaterialProperties.cpp @@ -522,8 +522,9 @@ TACSOrthotropicPly::TACSOrthotropicPly(TacsScalar _plyThickness, eYc = Yc / E2; eS12 = S12 / G12; - // By default, use Tsai-Wu - useTsaiWuCriterion = 1; + // By default, use Tsai-Wu modified to return a strength ratio + useTsaiWuCriterion = true; + useModifiedTsaiWu = true; // Compute the coefficients for the Tsai-Wu failure criteria F1 = (Xc - Xt) / (Xt * Xc); @@ -567,9 +568,15 @@ void TACSOrthotropicPly::setKSWeight(TacsScalar _ksWeight) { /* Set the failure criteria to use */ -void TACSOrthotropicPly::setUseMaxStrainCriterion() { useTsaiWuCriterion = 0; } +void TACSOrthotropicPly::setUseMaxStrainCriterion() { + useTsaiWuCriterion = false; +} + +void TACSOrthotropicPly::setUseTsaiWuCriterion() { useTsaiWuCriterion = true; } -void TACSOrthotropicPly::setUseTsaiWuCriterion() { useTsaiWuCriterion = 1; } +void TACSOrthotropicPly::setUseModifiedTsaiWu(bool _useModifiedTsaiWu) { + useModifiedTsaiWu = _useModifiedTsaiWu; +} /* Get the density of the material @@ -789,7 +796,13 @@ void TACSOrthotropicPly::calculateStress(TacsScalar angle, F11*s[0]**2 + F22*s[1]**2 + 2.0*F12*s[0]*s[1] + F66*s[2]**2 <= 1.0 - The value returned is actually 1/2 * (b + sqrt(b^2 + 4a)) where "b" + To enable this method, call `setUseTsaiWuCriterion()` and + `setUseModifiedTsaiWu(false)` + + 2. The Tsai-Wu strength ratio: + + This is similar to the Tsai-Wu criterion, except that the value + returned is actually 1/2 * (b + sqrt(b^2 + 4a)) where "b" is the linear part of the Tsai-Wu criterion and "a" is the quadratic part. This form is equivalent to the original form at the failure boundary but has the advantage that it scales lineary when the stress @@ -798,7 +811,9 @@ void TACSOrthotropicPly::calculateStress(TacsScalar angle, form. It also means that the failure criterion becomes equivalent to the von-Mises criterion when the material properties are isotropic. - 2. The maximum strain failure criteria: + This is the default method + + 3. The maximum strain failure criteria: KS(e_{i}/e_max^{+/-}, ksWeight) <= 1.0 @@ -820,7 +835,11 @@ TacsScalar TACSOrthotropicPly::failure(TacsScalar angle, linTerm = F1 * s[0] + F2 * s[1]; quadTerm = F11 * s[0] * s[0] + F22 * s[1] * s[1] + 2.0 * F12 * s[0] * s[1] + F66 * s[2] * s[2]; - fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); + if (useModifiedTsaiWu) { + fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); + } else { + fail = linTerm + quadTerm; + } } else { // Calculate the values of each of the failure criteria TacsScalar f[6]; @@ -866,22 +885,32 @@ TacsScalar TACSOrthotropicPly::failureStrainSens(TacsScalar angle, linTerm = F1 * s[0] + F2 * s[1]; quadTerm = F11 * s[0] * s[0] + F22 * s[1] * s[1] + 2.0 * F12 * s[0] * s[1] + F66 * s[2] * s[2]; - fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); - TacsScalar tmp = (F1 * s[0] + F2 * s[1]) * (F1 * s[0] + F2 * s[1]); - TacsScalar tmp2 = sqrt(4.0 * F11 * s2[0] + 8.0 * F12 * s[0] * s[1] + - 4.0 * F22 * s2[1] + 4.0 * F66 * s2[2] + tmp); + if (useModifiedTsaiWu) { + fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); + + TacsScalar tmp = (F1 * s[0] + F2 * s[1]) * (F1 * s[0] + F2 * s[1]); + TacsScalar tmp2 = sqrt(4.0 * F11 * s2[0] + 8.0 * F12 * s[0] * s[1] + + 4.0 * F22 * s2[1] + 4.0 * F66 * s2[2] + tmp); - // Calculate the sensitivity of the failure criteria w.r.t the 3 stresses - sens[0] = (F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + 4.0 * F11 * s[0] + - 4.0 * F12 * s[1]) / - (2.0 * tmp2); + // Calculate the sensitivity of the failure criteria w.r.t the 3 stresses + sens[0] = (F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + 4.0 * F11 * s[0] + + 4.0 * F12 * s[1]) / + (2.0 * tmp2); - sens[1] = (4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + F2 * tmp2 + - 4.0 * F22 * s[1]) / - (2.0 * tmp2); + sens[1] = (4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + F2 * tmp2 + + 4.0 * F22 * s[1]) / + (2.0 * tmp2); - sens[2] = 2.0 * F66 * s[2] / tmp2; + sens[2] = 2.0 * F66 * s[2] / tmp2; + } + + else { + fail = linTerm + quadTerm; + sens[0] = F1 + 2.0 * F11 * s[0] + 2.0 * F12 * s[1]; + sens[1] = F2 + 2.0 * F22 * s[1] + 2.0 * F12 * s[0]; + sens[2] = 2.0 * F66 * s[2]; + } TacsScalar sSens[3]; getPlyStress(sens, sSens); @@ -929,40 +958,48 @@ TacsScalar TACSOrthotropicPly::failureAngleSens(TacsScalar angle, TacsScalar e[3], se[3]; // The ply strain transformStrainGlobal2Ply(angle, strain, e); + // Compute the sensitivity of the transformation (se = d(e)/d(angle)) + transformStrainGlobal2PlyAngleSens(angle, strain, se); + TacsScalar fail = 0.0; if (useTsaiWuCriterion) { TacsScalar s[3], s2[3]; // The ply stress getPlyStress(e, s); - for (int ii = 0; ii < 3; ii++) { - s2[ii] = s[ii] * s[ii]; - } + // Compute the sensitivity of the stress + TacsScalar ss[3]; + getPlyStress(se, ss); TacsScalar linTerm, quadTerm; linTerm = F1 * s[0] + F2 * s[1]; quadTerm = F11 * s[0] * s[0] + F22 * s[1] * s[1] + 2.0 * F12 * s[0] * s[1] + F66 * s[2] * s[2]; - fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); - // Compute the sensitivity of the transformation (se = d(e)/d(angle)) - transformStrainGlobal2PlyAngleSens(angle, strain, se); + if (useModifiedTsaiWu) { + fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); + // ss = d(s)/d(e) * d(e)/d(angle) = d(s)/d(angle) - // Compute the sensitivity of the stress - TacsScalar ss[3]; - getPlyStress(se, ss); - // ss = d(s)/d(e) * d(e)/d(angle) = d(s)/d(angle) - - TacsScalar tmp = (F1 * s[0] + F2 * s[1]) * (F1 * s[0] + F2 * s[1]); - TacsScalar tmp2 = sqrt(4.0 * F11 * s2[0] + 8.0 * F12 * s[0] * s[1] + - 4.0 * F22 * s2[1] + 4.0 * F66 * s2[2] + tmp); - - *failSens = ss[0] * ((F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + - 4.0 * F11 * s[0] + 4.0 * F12 * s[1]) / - (2.0 * tmp2)); - *failSens += ss[1] * ((4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + - F2 * tmp2 + 4.0 * F22 * s[1]) / - (2.0 * tmp2)); - *failSens += ss[2] * (2.0 * F66 * s[2] / tmp2); + for (int ii = 0; ii < 3; ii++) { + s2[ii] = s[ii] * s[ii]; + } + TacsScalar tmp = (F1 * s[0] + F2 * s[1]) * (F1 * s[0] + F2 * s[1]); + TacsScalar tmp2 = sqrt(4.0 * F11 * s2[0] + 8.0 * F12 * s[0] * s[1] + + 4.0 * F22 * s2[1] + 4.0 * F66 * s2[2] + tmp); + + *failSens = ss[0] * ((F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + + 4.0 * F11 * s[0] + 4.0 * F12 * s[1]) / + (2.0 * tmp2)); + *failSens += ss[1] * ((4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + + F2 * tmp2 + 4.0 * F22 * s[1]) / + (2.0 * tmp2)); + *failSens += ss[2] * (2.0 * F66 * s[2] / tmp2); + } else { + fail = linTerm + quadTerm; + *failSens = + (2.0 * (F11 * s[0] * ss[0] + F22 * s[1] * ss[1] + + F12 * (ss[0] * s[1] + s[0] * ss[1]) + F66 * s[2] * ss[2]) + + F1 * ss[0] + F2 * ss[1]); + } } else { // Calculate the values of each of the failure criteria @@ -975,7 +1012,6 @@ TacsScalar TACSOrthotropicPly::failureAngleSens(TacsScalar angle, f[5] = -e[2] / eS12; // Compute the sensitivity of the transformation - transformStrainGlobal2PlyAngleSens(angle, strain, se); fs[0] = se[0] / eXt; fs[1] = -se[0] / eXc; fs[2] = se[1] / eYt; diff --git a/src/constitutive/TACSMaterialProperties.h b/src/constitutive/TACSMaterialProperties.h index 5b3f49b8c..4eb55d053 100644 --- a/src/constitutive/TACSMaterialProperties.h +++ b/src/constitutive/TACSMaterialProperties.h @@ -172,6 +172,9 @@ class TACSOrthotropicPly : public TACSObject { void setKSWeight(TacsScalar _ksWeight); void setUseMaxStrainCriterion(); void setUseTsaiWuCriterion(); + // Set whether to use the standard Tsai-Wu failure index or the modified + // version which returns a strength ratio + void setUseModifiedTsaiWu(bool _useModifiedTsaiWu); // Retrieve the material properties // -------------------------------- @@ -275,7 +278,8 @@ class TACSOrthotropicPly : public TACSObject { TacsScalar G12, G23, G13; // Keep track of which failure criterion to use - int useTsaiWuCriterion; + bool useTsaiWuCriterion; + bool useModifiedTsaiWu; // The stress-based strength properties TacsScalar Xt, Xc;