Skip to content

Commit

Permalink
Add back old Tsai-Wu failure index, add method to switch between the …
Browse files Browse the repository at this point in the history 
…two forms
  • Loading branch information
@A-CGray
A-CGray committed May 15, 2023
1 parent 188db1c commit 1760b99
Show file tree
Hide file tree
Showing 2 changed files with 83 additions and 43 deletions.
120 changes: 78 additions & 42 deletions src/constitutive/TACSMaterialProperties.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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];
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand All @@ -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;
Expand Down
6 changes: 5 additions & 1 deletion src/constitutive/TACSMaterialProperties.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
// --------------------------------
Expand Down Expand Up @@ -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;
Expand Down

0 comments on commit 1760b99

Please sign in to comment.