Skip to content

Commit

Permalink
Modify Tsai-Wu to work with safety factors (#203)
Browse files Browse the repository at this point in the history 
* Add __init__ method to BaseUI

* Modify problem classes to use __init__ from BaseUI

* Modify Tsai-Wu failure to work with safety factors

* Update integration tests with new failure values

* black .

* Remove commented out code and add comments explaining new Tsai-Wu calculation

* Add back old Tsai-Wu failure index, add method to switch between the two forms
  • Loading branch information
@A-CGray
A-CGray authored Jun 2, 2023
1 parent 3aea51c commit f3ddea0
Show file tree
Hide file tree
Showing 4 changed files with 119 additions and 34 deletions.
131 changes: 106 additions & 25 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,8 +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];
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 @@ -844,15 +875,42 @@ TacsScalar TACSOrthotropicPly::failureStrainSens(TacsScalar angle,

TacsScalar fail = 0.0;
if (useTsaiWuCriterion) {
TacsScalar s[3]; // Ply stress
TacsScalar s[3], s2[3]; // Ply stress
getPlyStress(e, s);
for (int ii = 0; ii < 3; ii++) {
s2[ii] = s[ii] * s[ii];
}

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

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

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]);
// 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[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];
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 @@ -900,25 +958,49 @@ 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]; // The ply stress
TacsScalar s[3], s2[3]; // The 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]);

// Compute the sensitivity of the transformation
transformStrainGlobal2PlyAngleSens(angle, strain, se);

// Compute the sensitivity of the stress
TacsScalar ss[3];
getPlyStress(se, ss);

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

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)

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
TacsScalar f[6], fs[6];
Expand All @@ -930,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 f3ddea0

Please sign in to comment.