diff --git a/src/Pybind11Wraps/SPH/SolidSPHHydroBase.py b/src/Pybind11Wraps/SPH/SolidSPHHydroBase.py index 007fa24ff..11dc9dbec 100644 --- a/src/Pybind11Wraps/SPH/SolidSPHHydroBase.py +++ b/src/Pybind11Wraps/SPH/SolidSPHHydroBase.py @@ -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" @@ -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") diff --git a/src/Pybind11Wraps/SPH/SolidSPHHydroBaseRZ.py b/src/Pybind11Wraps/SPH/SolidSPHHydroBaseRZ.py index f887ae078..0af0e0ed0 100644 --- a/src/Pybind11Wraps/SPH/SolidSPHHydroBaseRZ.py +++ b/src/Pybind11Wraps/SPH/SolidSPHHydroBaseRZ.py @@ -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" diff --git a/src/Pybind11Wraps/SolidMaterial/CollinsStrength.py b/src/Pybind11Wraps/SolidMaterial/CollinsStrength.py index c2d1fc43d..894f5041e 100644 --- a/src/Pybind11Wraps/SolidMaterial/CollinsStrength.py +++ b/src/Pybind11Wraps/SolidMaterial/CollinsStrength.py @@ -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""" @@ -55,6 +70,7 @@ def soundSpeed(self, #........................................................................... # Properties mui = PYB11property("double") + mud = PYB11property("double") Y0 = PYB11property("double") Ym = PYB11property("double") diff --git a/src/SPH/SPHHydros.py b/src/SPH/SPHHydros.py index 6acd3edea..be99ff697 100644 --- a/src/SPH/SPHHydros.py +++ b/src/SPH/SPHHydros.py @@ -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() @@ -112,6 +113,7 @@ def __init__(self, nTensile, damageRelieveRubble, negativePressureInDamage, + strengthInDamage, xmin, xmax) return @@ -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): @@ -286,6 +289,7 @@ def __init__(self, nTensile, damageRelieveRubble, negativePressureInDamage, + strengthInDamage, xmin, xmax) self.zaxisBC = AxisBoundaryRZ(etaMinAxis) @@ -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, @@ -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) @@ -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, @@ -475,6 +482,7 @@ def ASPH(dataBase, nTensile = nTensile, damageRelieveRubble = damageRelieveRubble, negativePressureInDamage = negativePressureInDamage, + strengthInDamage = strengthInDamage, xmin = xmin, xmax = xmax, ASPH = True) @@ -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, @@ -524,6 +533,7 @@ def SPHRZ(dataBase, nTensile = nTensile, damageRelieveRubble = damageRelieveRubble, negativePressureInDamage = negativePressureInDamage, + strengthInDamage = strengthInDamage, xmin = xmin, xmax = xmax, ASPH = ASPH, @@ -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): @@ -573,6 +584,7 @@ def ASPHRZ(dataBase, nTensile = nTensile, damageRelieveRubble = damageRelieveRubble, negativePressureInDamage = negativePressureInDamage, + strengthInDamage = strengthInDamage, xmin = xmin, xmax = xmax, ASPH = True, diff --git a/src/SPH/SolidSPHEvaluateDerivatives.cc b/src/SPH/SolidSPHEvaluateDerivatives.cc index 8ee561323..1a0b1f23d 100644 --- a/src/SPH/SolidSPHEvaluateDerivatives.cc +++ b/src/SPH/SolidSPHEvaluateDerivatives.cc @@ -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 diff --git a/src/SPH/SolidSPHHydroBase.cc b/src/SPH/SolidSPHHydroBase.cc index 8af1e4c8c..b752ebc45 100644 --- a/src/SPH/SolidSPHHydroBase.cc +++ b/src/SPH/SolidSPHHydroBase.cc @@ -119,6 +119,7 @@ SolidSPHHydroBase(const SmoothingScaleBase& smoothingScaleMethod, const double nTensile, const bool damageRelieveRubble, const bool negativePressureInDamage, + const bool strengthInDamage, const Vector& xmin, const Vector& xmax): SPHHydroBase(smoothingScaleMethod, @@ -143,6 +144,7 @@ SolidSPHHydroBase(const SmoothingScaleBase& smoothingScaleMethod, xmax), mDamageRelieveRubble(damageRelieveRubble), mNegativePressureInDamage(negativePressureInDamage), + mStrengthInDamage(strengthInDamage), mGradKernel(WGrad), mDdeviatoricStressDt(FieldStorageType::CopyFields), mBulkModulus(FieldStorageType::CopyFields), diff --git a/src/SPH/SolidSPHHydroBase.hh b/src/SPH/SolidSPHHydroBase.hh index 8853b18ab..941f1c552 100644 --- a/src/SPH/SolidSPHHydroBase.hh +++ b/src/SPH/SolidSPHHydroBase.hh @@ -57,6 +57,7 @@ public: const double nTensile, const bool damageRelieveRubble, const bool negativePressureInDamage, + const bool strengthInDamage, const Vector& xmin, const Vector& xmax); @@ -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"; } @@ -124,7 +129,7 @@ public: protected: //--------------------------- Protected Interface ---------------------------// - bool mDamageRelieveRubble, mNegativePressureInDamage; + bool mDamageRelieveRubble, mNegativePressureInDamage, mStrengthInDamage; private: //--------------------------- Private Interface ---------------------------// diff --git a/src/SPH/SolidSPHHydroBaseInline.hh b/src/SPH/SolidSPHHydroBaseInline.hh index 5917374fc..b49e186ef 100644 --- a/src/SPH/SolidSPHHydroBaseInline.hh +++ b/src/SPH/SolidSPHHydroBaseInline.hh @@ -38,6 +38,25 @@ negativePressureInDamage(bool x) { mNegativePressureInDamage = x; } +//------------------------------------------------------------------------------ +// Do we allow damaged material to have strength? +//------------------------------------------------------------------------------ +template +inline +bool +SolidSPHHydroBase:: +strengthInDamage() const { + return mStrengthInDamage; +} + +template +inline +void +SolidSPHHydroBase:: +strengthInDamage(bool x) { + mStrengthInDamage = x; +} + //------------------------------------------------------------------------------ // Kernel for velocity gradient. //------------------------------------------------------------------------------ diff --git a/src/SPH/SolidSPHHydroBaseRZ.cc b/src/SPH/SolidSPHHydroBaseRZ.cc index 3b3813258..aa7fd15a8 100644 --- a/src/SPH/SolidSPHHydroBaseRZ.cc +++ b/src/SPH/SolidSPHHydroBaseRZ.cc @@ -96,6 +96,7 @@ SolidSPHHydroBaseRZ(const SmoothingScaleBase >& smoothingScaleMethod, const double nTensile, const bool damageRelieveRubble, const bool negativePressureInDamage, + const bool strengthInDamage, const Vector& xmin, const Vector& xmax): SolidSPHHydroBase >(smoothingScaleMethod, @@ -119,6 +120,7 @@ SolidSPHHydroBaseRZ(const SmoothingScaleBase >& smoothingScaleMethod, nTensile, damageRelieveRubble, negativePressureInDamage, + strengthInDamage, xmin, xmax) { } @@ -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 diff --git a/src/SPH/SolidSPHHydroBaseRZ.hh b/src/SPH/SolidSPHHydroBaseRZ.hh index ec1210ece..e1fa70871 100644 --- a/src/SPH/SolidSPHHydroBaseRZ.hh +++ b/src/SPH/SolidSPHHydroBaseRZ.hh @@ -65,6 +65,7 @@ public: const double nTensile, const bool damageRelieveRubble, const bool negativePressureInDamage, + const bool strengthInDamage, const Vector& xmin, const Vector& xmax); diff --git a/src/SolidMaterial/CollinsStrength.cc b/src/SolidMaterial/CollinsStrength.cc index 3950d451d..cdba556fd 100644 --- a/src/SolidMaterial/CollinsStrength.cc +++ b/src/SolidMaterial/CollinsStrength.cc @@ -16,6 +16,7 @@ #include "Utilities/DBC.hh" #include "SolidEquationOfState.hh" #include "Field/Field.hh" +#include "NodeList/SolidNodeList.hh" namespace Spheral { @@ -30,15 +31,35 @@ template CollinsStrength:: CollinsStrength(const StrengthModel& shearModulusModel, const double mui, + const double mud, const double Y0, const double Ym): StrengthModel(), mShearModulusModel(shearModulusModel), mmui(mui), + mmud(mud), mY0(Y0), mYm(Ym) { } +//------------------------------------------------------------------------------ +// Constructor. +//------------------------------------------------------------------------------ +template +CollinsStrength:: +CollinsStrength(const StrengthModel& shearModulusModel, + const double mui, + const double Y0, + const double Ym): + StrengthModel(), + 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. //------------------------------------------------------------------------------ @@ -72,10 +93,23 @@ yieldStrength(Field& yieldStrength, const Field& pressure, const Field& /*plasticStrain*/, const Field& /*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&>(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); } } @@ -110,6 +144,13 @@ mui() const { return mmui; } +template +double +CollinsStrength:: +mud() const { + return mmud; +} + template double CollinsStrength:: diff --git a/src/SolidMaterial/CollinsStrength.hh b/src/SolidMaterial/CollinsStrength.hh index 8bb3bc63a..809bc903a 100644 --- a/src/SolidMaterial/CollinsStrength.hh +++ b/src/SolidMaterial/CollinsStrength.hh @@ -26,12 +26,19 @@ public: typedef typename Dimension::Scalar Scalar; // Constructors, destructor. + CollinsStrength(const StrengthModel& shearModulusModel, // Used to compute the shear modulus + const double mui, // Coefficient of internal friction (intact material) + const double mud, // Coefficient of internal friction (damaged material) + const double Y0, // Shear strength at zero pressure + const double Ym); // von Mises plastic limit + // Backwards compatible constructor without the mud term. To be deprecated... CollinsStrength(const StrengthModel& shearModulusModel, // Used to compute the shear modulus const double mui, // Coefficient of internal friction const double Y0, // Shear strength at zero pressure const double Ym); // von Mises plastic limit virtual ~CollinsStrength(); + // Override the required generic interface. virtual bool providesSoundSpeed() const override { return mShearModulusModel.providesSoundSpeed(); } virtual void shearModulus(Field& shearModulus, @@ -55,13 +62,14 @@ public: // Access the strength parameters. const StrengthModel& shearModulusModel() const; double mui() const; + double mud() const; double Y0() const; double Ym() const; private: //--------------------------- Private Interface ---------------------------// const StrengthModel& mShearModulusModel; - double mmui, mY0, mYm; + double mmui, mmud, mY0, mYm; // No copying or assignment. CollinsStrength(const CollinsStrength&);