Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/damage #33

Merged
merged 4 commits into from
Oct 8, 2020
Merged

Feature/damage #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6299659
Initial implementation of damage extension to Collins material strength.
jmikeowen Sep 30, 2020
7d27973
Adding a "strengthInDamage" option to the solid SPH hydro package.
jmikeowen Sep 30, 2020
25fe694
Typo!
jmikeowen Oct 6, 2020
bc4b2d1
Ensuring yield strength is always positive (limiting negative pressur…
jmikeowen Oct 8, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion src/Pybind11Wraps/SPH/SolidSPHHydroBase.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ def pyinit(smoothingScaleMethod = "const SmoothingScaleBase<%(Dimension)s>&",
nTensile = "const double",
damageRelieveRubble = "const bool",
negativePressureInDamage = "const bool",
strengthInDamage = "const bool",
xmin = "const Vector&",
xmax = "const Vector&"):
"SolidSPHHydroBase constructor"
Expand Down Expand Up @@ -92,7 +93,9 @@ def enforceBoundaries(state = "State<%(Dimension)s>&",
damageRelieveRubble = PYB11property("bool", "damageRelieveRubble", "damageRelieveRubble",
doc="Control whether allow damaged material to have stress relieved.")
negativePressureInDamage = PYB11property("bool", "negativePressureInDamage", "negativePressureInDamage",
doc="Should we allow damage material to support negative pressures?")
doc="Should we allow damaged material to support negative pressures?")
strengthInDamage = PYB11property("bool", "strengthInDamage", "strengthInDamage",
doc="Should we allow damaged material to support strength?")

DdeviatoricStressDt = PYB11property("const FieldList<%(Dimension)s, SymTensor>&", "DdeviatoricStressDt", returnpolicy="reference_internal")
bulkModulus = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "bulkModulus", returnpolicy="reference_internal")
Expand Down
1 change: 1 addition & 0 deletions src/Pybind11Wraps/SPH/SolidSPHHydroBaseRZ.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ def pyinit(smoothingScaleMethod = "const SmoothingScaleBase<%(Dimension)s>&",
nTensile = "const double",
damageRelieveRubble = "const bool",
negativePressureInDamage = "const bool",
strengthInDamage = "const bool",
xmin = "const Vector&",
xmax = "const Vector&"):
"SolidSPHHydroBaseRZ constructor"
Expand Down
16 changes: 16 additions & 0 deletions src/Pybind11Wraps/SolidMaterial/CollinsStrength.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,25 @@ class CollinsStrength(StrengthModel):
def pyinit(self,
shearModulusModel = "const StrengthModel<%(Dimension)s>&",
mui = "const double",
mud = "const double",
Y0 = "const double",
Ym = "const double"):
"""CollinsStrength constructor

mui : Coefficient of internal friction (intact material)
mud : Coefficient of internal friction (damaged material)
Y0 : Shear strength at zero pressure
yM : von Mises plastic limit"""

def pyinit1(self,
shearModulusModel = "const StrengthModel<%(Dimension)s>&",
mui = "const double",
Y0 = "const double",
Ym = "const double"):
"""CollinsStrength constructor

DEPRECATION WARNING: this version without mud is deprecated

mui : Coefficient of internal friction
Y0 : Shear strength at zero pressure
yM : von Mises plastic limit"""
Expand All @@ -55,6 +70,7 @@ def soundSpeed(self,
#...........................................................................
# Properties
mui = PYB11property("double")
mud = PYB11property("double")
Y0 = PYB11property("double")
Ym = PYB11property("double")

Expand Down
16 changes: 14 additions & 2 deletions src/SPH/SPHHydros.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,7 @@ def __init__(self,
nTensile = 4.0,
damageRelieveRubble = False,
negativePressureInDamage = False,
strengthInDamage = False,
xmin = Vector%(dim)s(-1e100, -1e100, -1e100),
xmax = Vector%(dim)s( 1e100, 1e100, 1e100)):
self._smoothingScaleMethod = %(smoothingScaleMethod)s%(dim)s()
Expand Down Expand Up @@ -112,6 +113,7 @@ def __init__(self,
nTensile,
damageRelieveRubble,
negativePressureInDamage,
strengthInDamage,
xmin,
xmax)
return
Expand Down Expand Up @@ -256,6 +258,7 @@ def __init__(self,
nTensile = 4.0,
damageRelieveRubble = False,
negativePressureInDamage = False,
strengthInDamage = False,
xmin = Vector2d(-1e100, -1e100),
xmax = Vector2d( 1e100, 1e100),
etaMinAxis = 0.1):
Expand Down Expand Up @@ -286,6 +289,7 @@ def __init__(self,
nTensile,
damageRelieveRubble,
negativePressureInDamage,
strengthInDamage,
xmin,
xmax)
self.zaxisBC = AxisBoundaryRZ(etaMinAxis)
Expand Down Expand Up @@ -349,6 +353,7 @@ def SPH(dataBase,
nTensile = 4.0,
damageRelieveRubble = False,
negativePressureInDamage = False,
strengthInDamage = False,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100),
ASPH = False,
Expand Down Expand Up @@ -424,8 +429,9 @@ def SPH(dataBase,
"xmax" : eval("Vector%id(%g, %g, %g)" % xmax)}

if nsolid > 0:
kwargs.update({"damageRelieveRubble" : damageRelieveRubble,
"negativePressureInDamage" : negativePressureInDamage})
kwargs.update({"damageRelieveRubble" : damageRelieveRubble,
"negativePressureInDamage" : negativePressureInDamage,
"strengthInDamage" : strengthInDamage})

# Build and return the thing.
result = Constructor(**kwargs)
Expand Down Expand Up @@ -454,6 +460,7 @@ def ASPH(dataBase,
nTensile = 4.0,
damageRelieveRubble = False,
negativePressureInDamage = False,
strengthInDamage = False,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100)):
return SPH(dataBase = dataBase,
Expand All @@ -475,6 +482,7 @@ def ASPH(dataBase,
nTensile = nTensile,
damageRelieveRubble = damageRelieveRubble,
negativePressureInDamage = negativePressureInDamage,
strengthInDamage = strengthInDamage,
xmin = xmin,
xmax = xmax,
ASPH = True)
Expand All @@ -501,6 +509,7 @@ def SPHRZ(dataBase,
nTensile = 4.0,
damageRelieveRubble = False,
negativePressureInDamage = False,
strengthInDamage = False,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100),
etaMinAxis = 0.1,
Expand All @@ -524,6 +533,7 @@ def SPHRZ(dataBase,
nTensile = nTensile,
damageRelieveRubble = damageRelieveRubble,
negativePressureInDamage = negativePressureInDamage,
strengthInDamage = strengthInDamage,
xmin = xmin,
xmax = xmax,
ASPH = ASPH,
Expand Down Expand Up @@ -551,6 +561,7 @@ def ASPHRZ(dataBase,
nTensile = 4.0,
damageRelieveRubble = False,
negativePressureInDamage = False,
strengthInDamage = False,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100),
etaMinAxis = 0.1):
Expand All @@ -573,6 +584,7 @@ def ASPHRZ(dataBase,
nTensile = nTensile,
damageRelieveRubble = damageRelieveRubble,
negativePressureInDamage = negativePressureInDamage,
strengthInDamage = strengthInDamage,
xmin = xmin,
xmax = xmax,
ASPH = True,
Expand Down
14 changes: 9 additions & 5 deletions src/SPH/SolidSPHEvaluateDerivatives.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -307,12 +307,16 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/,
const auto Peffj = (mNegativePressureInDamage or Pj > 0.0 ? Pj : fDeffij*Pj);

// Compute the stress tensors.
sigmai = -Peffi*SymTensor::one;
sigmaj = -Peffj*SymTensor::one;
if (sameMatij) {
sigmai = fDeffij*Si - Peffi*SymTensor::one;
sigmaj = fDeffij*Sj - Peffj*SymTensor::one;
} else {
sigmai = -Peffi*SymTensor::one;
sigmaj = -Peffj*SymTensor::one;
if (mStrengthInDamage) {
sigmai += Si;
sigmaj += Sj;
} else {
sigmai += fDeffij*Si;
sigmaj += fDeffij*Sj;
}
}

// Compute the tensile correction to add to the stress as described in
Expand Down
2 changes: 2 additions & 0 deletions src/SPH/SolidSPHHydroBase.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ SolidSPHHydroBase(const SmoothingScaleBase<Dimension>& smoothingScaleMethod,
const double nTensile,
const bool damageRelieveRubble,
const bool negativePressureInDamage,
const bool strengthInDamage,
const Vector& xmin,
const Vector& xmax):
SPHHydroBase<Dimension>(smoothingScaleMethod,
Expand All @@ -143,6 +144,7 @@ SolidSPHHydroBase(const SmoothingScaleBase<Dimension>& smoothingScaleMethod,
xmax),
mDamageRelieveRubble(damageRelieveRubble),
mNegativePressureInDamage(negativePressureInDamage),
mStrengthInDamage(strengthInDamage),
mGradKernel(WGrad),
mDdeviatoricStressDt(FieldStorageType::CopyFields),
mBulkModulus(FieldStorageType::CopyFields),
Expand Down
7 changes: 6 additions & 1 deletion src/SPH/SolidSPHHydroBase.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ public:
const double nTensile,
const bool damageRelieveRubble,
const bool negativePressureInDamage,
const bool strengthInDamage,
const Vector& xmin,
const Vector& xmax);

Expand Down Expand Up @@ -115,6 +116,10 @@ public:
bool negativePressureInDamage() const;
void negativePressureInDamage(bool x);

// Do we allow damaged material to have strength?
bool strengthInDamage() const;
void strengthInDamage(bool x);

//****************************************************************************
// Methods required for restarting.
virtual std::string label() const override { return "SolidSPHHydroBase"; }
Expand All @@ -124,7 +129,7 @@ public:

protected:
//--------------------------- Protected Interface ---------------------------//
bool mDamageRelieveRubble, mNegativePressureInDamage;
bool mDamageRelieveRubble, mNegativePressureInDamage, mStrengthInDamage;

private:
//--------------------------- Private Interface ---------------------------//
Expand Down
19 changes: 19 additions & 0 deletions src/SPH/SolidSPHHydroBaseInline.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,25 @@ negativePressureInDamage(bool x) {
mNegativePressureInDamage = x;
}

//------------------------------------------------------------------------------
// Do we allow damaged material to have strength?
//------------------------------------------------------------------------------
template<typename Dimension>
inline
bool
SolidSPHHydroBase<Dimension>::
strengthInDamage() const {
return mStrengthInDamage;
}

template<typename Dimension>
inline
void
SolidSPHHydroBase<Dimension>::
strengthInDamage(bool x) {
mStrengthInDamage = x;
}

//------------------------------------------------------------------------------
// Kernel for velocity gradient.
//------------------------------------------------------------------------------
Expand Down
16 changes: 11 additions & 5 deletions src/SPH/SolidSPHHydroBaseRZ.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,7 @@ SolidSPHHydroBaseRZ(const SmoothingScaleBase<Dim<2> >& smoothingScaleMethod,
const double nTensile,
const bool damageRelieveRubble,
const bool negativePressureInDamage,
const bool strengthInDamage,
const Vector& xmin,
const Vector& xmax):
SolidSPHHydroBase<Dim<2> >(smoothingScaleMethod,
Expand All @@ -119,6 +120,7 @@ SolidSPHHydroBaseRZ(const SmoothingScaleBase<Dim<2> >& smoothingScaleMethod,
nTensile,
damageRelieveRubble,
negativePressureInDamage,
strengthInDamage,
xmin,
xmax) {
}
Expand Down Expand Up @@ -540,12 +542,16 @@ evaluateDerivatives(const Dim<2>::Scalar /*time*/,
const auto Peffj = (mNegativePressureInDamage or Pj > 0.0 ? Pj : fDeffij*Pj);

// Compute the stress tensors.
sigmai = -Peffi*SymTensor::one;
sigmaj = -Peffj*SymTensor::one;
if (sameMatij) {
sigmai = fDeffij*Si - Peffi*SymTensor::one;
sigmaj = fDeffij*Sj - Peffj*SymTensor::one;
} else {
sigmai = -Peffi*SymTensor::one;
sigmaj = -Peffj*SymTensor::one;
if (mStrengthInDamage) {
sigmai += Si;
sigmaj += Sj;
} else {
sigmai += fDeffij*Si;
sigmaj += fDeffij*Sj;
}
}

// Compute the tensile correction to add to the stress as described in
Expand Down
1 change: 1 addition & 0 deletions src/SPH/SolidSPHHydroBaseRZ.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@ public:
const double nTensile,
const bool damageRelieveRubble,
const bool negativePressureInDamage,
const bool strengthInDamage,
const Vector& xmin,
const Vector& xmax);

Expand Down
47 changes: 44 additions & 3 deletions src/SolidMaterial/CollinsStrength.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "Utilities/DBC.hh"
#include "SolidEquationOfState.hh"
#include "Field/Field.hh"
#include "NodeList/SolidNodeList.hh"

namespace Spheral {

Expand All @@ -30,15 +31,35 @@ template<typename Dimension>
CollinsStrength<Dimension>::
CollinsStrength(const StrengthModel<Dimension>& shearModulusModel,
const double mui,
const double mud,
const double Y0,
const double Ym):
StrengthModel<Dimension>(),
mShearModulusModel(shearModulusModel),
mmui(mui),
mmud(mud),
mY0(Y0),
mYm(Ym) {
}

//------------------------------------------------------------------------------
// Constructor.
//------------------------------------------------------------------------------
template<typename Dimension>
CollinsStrength<Dimension>::
CollinsStrength(const StrengthModel<Dimension>& shearModulusModel,
const double mui,
const double Y0,
const double Ym):
StrengthModel<Dimension>(),
mShearModulusModel(shearModulusModel),
mmui(mui),
mmud(0.0),
mY0(Y0),
mYm(Ym) {
if (Process::getRank() == 0) printf("Deprecation WARNING: specifying the Collins strength model without the coefficient of friction in damage (mud) is deprecated.\n");
}

//------------------------------------------------------------------------------
// Destructor.
//------------------------------------------------------------------------------
Expand Down Expand Up @@ -72,10 +93,23 @@ yieldStrength(Field<Dimension, Scalar>& yieldStrength,
const Field<Dimension, Scalar>& pressure,
const Field<Dimension, Scalar>& /*plasticStrain*/,
const Field<Dimension, Scalar>& /*plasticStrainRate*/) const {
const unsigned n = density.numInternalElements();
const Scalar YdiffInv = safeInvVar(mYm - mY0);

// Do a janky extraction of the damage (unknown time level) from the NodeList. Better
// be a SolidNodeList since we're using strength!
const auto& nodes = dynamic_cast<const SolidNodeList<Dimension>&>(yieldStrength.nodeList());
const auto& Dfield = nodes.damage();

const auto n = density.numInternalElements();
const auto YdiffInv = safeInvVar(mYm - mY0);
for (unsigned i = 0; i != n; ++i) {
yieldStrength(i) = mY0 + mmui*pressure(i)/(1.0 + mmui*pressure(i)*YdiffInv);
const auto Pi = std::max(0.0, pressure(i));
const auto Yi = mY0 + mmui*Pi/(1.0 + mmui*Pi*YdiffInv);
const auto Yd = std::min(Yi, mmud*Pi);
CHECK(Yi >= 0.0 and Yd >= 0.0);
const auto Di = Dfield(i).Trace()/Dimension::nDim;
CHECK(Di >= 0.0 and Di <= 1.0);
yieldStrength(i) = (1.0 - Di)*Yi + Di*Yd;
CHECK(yieldStrength(i) >= 0.0);
}
}

Expand Down Expand Up @@ -110,6 +144,13 @@ mui() const {
return mmui;
}

template<typename Dimension>
double
CollinsStrength<Dimension>::
mud() const {
return mmud;
}

template<typename Dimension>
double
CollinsStrength<Dimension>::
Expand Down
Loading