Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/TsaiWuMod' into bladeStiffenedShell
Browse files Browse the repository at this point in the history
  • Loading branch information
@A-CGray
A-CGray committed May 19, 2023
2 parents 9c20101 + 1760b99 commit 69dbcbb
Show file tree
Hide file tree
Showing 4 changed files with 99 additions and 64 deletions.
141 changes: 86 additions & 55 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,24 @@ 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
2. The maximum strain failure criteria:
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
state is uniformly scaled. This means that the failure criterion can be
used to compute safety factors which is not true of the original quadratic
form. It also means that the failure criterion becomes equivalent to the
von-Mises criterion when the material properties are isotropic.
This is the default method
3. The maximum strain failure criteria:
KS(e_{i}/e_max^{+/-}, ksWeight) <= 1.0
Expand All @@ -807,14 +831,15 @@ TacsScalar TACSOrthotropicPly::failure(TacsScalar angle,
TacsScalar s[3]; // Ply stress
getPlyStress(e, s);

// fail = (F11 * s[0] * s[0] + F22 * s[1] * s[1] + 2.0 * F12 * s[0] * s[1] +
// F66 * s[2] * s[2] + F1 * s[0] + F2 * s[1]);

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));
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 @@ -856,30 +881,36 @@ TacsScalar TACSOrthotropicPly::failureStrainSens(TacsScalar angle,
s2[ii] = s[ii] * s[ii];
}

// fail = (F11 * s[0] * s[0] + F22 * s[1] * s[1] + 2.0 * F12 * s[0] * s[1] +
// F66 * s[2] * s[2] + F1 * s[0] + F2 * s[1]);
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));

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));

// sens[0] = F1 + 2.0 * F11 * s[0] + 2.0 * F12 * s[1];
sens[0] = (F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + 4.0 * F11 * s[0] +
4.0 * F12 * s[1]) /
(2.0 * tmp2);
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);

// sens[1] = F2 + 2.0 * F22 * s[1] + 2.0 * F12 * s[0];
sens[1] = (4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + F2 * tmp2 +
4.0 * F22 * 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[2] = 2.0 * F66 * s[2];
sens[2] = 2.0 * F66 * s[2] / 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;
}

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 @@ -927,47 +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);

// fail = (F11 * s[0] * s[0] + F22 * s[1] * s[1] + 2.0 * F12 * s[0] * s[1] +
// F66 * s[2] * s[2] + F1 * s[0] + F2 * s[1]);
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 =
// (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]);

*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 @@ -980,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
6 changes: 2 additions & 4 deletions tests/integration_tests/test_shell_beam_comp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,12 @@
from tacs import TACS, elements, constitutive, functions

"""
Create a cantilevered composite beam of linear quad shells with an
Create a cantilevered composite beam of linear quad shells with an
unbalanced tip shear load on the right corner
and test KSFailure, StructuralMass, and Compliance functions and sensitivities
"""

FUNC_REFS = np.array(
[5812.5, 63907185.558059536, 12.21799417804536, 174.71401901274177]
)
FUNC_REFS = np.array([5812.5, 63907185.558059536, 12.21799417804536, 23.51274876481848])

# Length of plate in x/y direction
Lx = 10.0
Expand Down
10 changes: 6 additions & 4 deletions tests/integration_tests/test_shell_comp_unbalanced.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
Tests an unbalanced laminate shell model with the following layup: [0, 45, 30]s.
The laminate information is read in from a PCOMP card in the BDF file.
Two load cases are tested: an in-plane tension and out-of-plane shear.
tests KSDisplacement, KSFailure, StructuralMass, and Compliance functions
tests KSDisplacement, KSFailure, StructuralMass, and Compliance functions
and sensitivities.
"""

Expand All @@ -22,13 +22,13 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):

FUNC_REFS = {
"Tension_compliance": 15047.4827204001,
"Tension_ks_vmfailure": 34.54666371379912,
"Tension_ks_TsaiWufailure": 7.875796240395447,
"Tension_mass": 1.1625,
"Tension_x_disp": 0.08602975190111069,
"Tension_y_disp": 0.01957454912511978,
"Tension_z_disp": 5.484526868562824e-17,
"VertShear_compliance": 0.00020961674292023337,
"VertShear_ks_vmfailure": 0.0007498886916845651,
"VertShear_ks_TsaiWufailure": 0.0013092603829596005,
"VertShear_mass": 1.1625,
"VertShear_x_disp": 2.6216565960935596e-23,
"VertShear_y_disp": 5.97480777083504e-23,
Expand Down Expand Up @@ -71,7 +71,9 @@ def setup_tacs_problems(self, comm):
for problem in tacs_probs:
problem.addFunction("mass", functions.StructuralMass)
problem.addFunction("compliance", functions.Compliance)
problem.addFunction("ks_vmfailure", functions.KSFailure, ksWeight=ksweight)
problem.addFunction(
"ks_TsaiWufailure", functions.KSFailure, ksWeight=ksweight
)
problem.addFunction(
"x_disp",
functions.KSDisplacement,
Expand Down

0 comments on commit 69dbcbb

Please sign in to comment.