diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh index e3f1d27c92caec76d4d0065fac194e3d9f0614e1..3ac28d2ade2a9c4ced7213872daa665099a4ac8a 100644 --- a/dune/tectonic/frictionpotential.hh +++ b/dune/tectonic/frictionpotential.hh @@ -69,6 +69,38 @@ class TruncatedRateState : public FrictionPotential { double Vmin; }; +class RegularisedRateState : public FrictionPotential { +public: + RegularisedRateState(double coefficient, double _normalStress, + FrictionData _fd) + : fd(_fd), weight(coefficient), normalStress(_normalStress) {} + + double coefficientOfFriction(double V) const { + return fd.a * std::asinh(0.5 * V / Vmin); + } + + double differential(double V) const { + return weight * (fd.C - normalStress * coefficientOfFriction(V)); + } + + double second_deriv(double V) const { + return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V); + } + + double regularity(double V) const { return std::abs(second_deriv(V)); } + + void updateAlpha(double alpha) { + double const logrest = (fd.mu0 + fd.b * alpha) / fd.a; + Vmin = fd.V0 / std::exp(logrest); + } + +private: + FrictionData const fd; + double const weight; + double const normalStress; + double Vmin; +}; + class TrivialFunction : public FrictionPotential { public: double evaluate(double) const { return 0; } diff --git a/src/assemblers.cc b/src/assemblers.cc index 2360258342a29cec388d8107053d7978a1353abc..c2a9fae984b1aa9154ab6f17c001b36ed4346a45 100644 --- a/src/assemblers.cc +++ b/src/assemblers.cc @@ -122,6 +122,7 @@ void MyAssembler<GridView, dimension>::assembleNormalStress( template <class GridView, int dimension> auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( + Config::FrictionModel frictionModel, BoundaryPatch<GridView> const &frictionalBoundary, GlobalFrictionData<dimension> const &frictionInfo, ScalarVector const &normalStress) @@ -136,9 +137,18 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler, weights, frictionalBoundary); } - return std::make_shared< - GlobalRateStateFriction<Matrix, Vector, TruncatedRateState, GridView>>( - frictionalBoundary, gridView, frictionInfo, weights, normalStress); + switch (frictionModel) { + case Config::Truncated: + return std::make_shared<GlobalRateStateFriction< + Matrix, Vector, TruncatedRateState, GridView>>( + frictionalBoundary, gridView, frictionInfo, weights, normalStress); + case Config::Regularised: + return std::make_shared<GlobalRateStateFriction< + Matrix, Vector, RegularisedRateState, GridView>>( + frictionalBoundary, gridView, frictionInfo, weights, normalStress); + default: + assert(false); + } } template <class GridView, int dimension> diff --git a/src/assemblers.hh b/src/assemblers.hh index 7b0c26c92a079f237c343b0dfc86e9f4e03d59fc..207bbde220327e999d418f2906de7e09f467a99f 100644 --- a/src/assemblers.hh +++ b/src/assemblers.hh @@ -16,6 +16,8 @@ #include <dune/tectonic/globalfriction.hh> #include <dune/tectonic/globalfrictiondata.hh> +#include "enums.hh" + template <class GridView, int dimension> class MyAssembler { public: using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>; @@ -76,6 +78,7 @@ template <class GridView, int dimension> class MyAssembler { double poissonRatio, Vector const &displacement); std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity( + Config::FrictionModel frictionModel, BoundaryPatch<GridView> const &frictionalBoundary, GlobalFrictionData<dimension> const &frictionInfo, ScalarVector const &normalStress); diff --git a/src/enum_friction_model.cc b/src/enum_friction_model.cc new file mode 100644 index 0000000000000000000000000000000000000000..e8a7d066a44a61790ee2c798bd2eb2facb3c3149 --- /dev/null +++ b/src/enum_friction_model.cc @@ -0,0 +1,15 @@ +#include <dune/common/exceptions.hh> + +#include "enums.hh" + +template <> struct StringToEnum<Config::FrictionModel> { + static Config::FrictionModel convert(std::string const &s) { + if (s == "Truncated") + return Config::Truncated; + + if (s == "Regularised") + return Config::Regularised; + + DUNE_THROW(Dune::Exception, "failed to parse enum"); + } +}; diff --git a/src/enums.hh b/src/enums.hh index 7dbd9b95997e0076f4bf078112991b31971c8fd4..5c04e4a41fc3ed756434e65cbcf81f309177d141 100644 --- a/src/enums.hh +++ b/src/enums.hh @@ -2,6 +2,10 @@ #define SRC_ENUMS_HH struct Config { + enum FrictionModel { + Truncated, + Regularised + }; enum stateModel { AgeingLaw, SlipLaw diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg index 8a2b3bd9d2c27e1af76e1a27d9953619038edbce..9425327f7712cf2319759f8d5910e982c81c5918 100644 --- a/src/sand-wedge-data/parset.cfg +++ b/src/sand-wedge-data/parset.cfg @@ -27,6 +27,7 @@ V0 = 5e-5 # [m/s] L = 2e-5 # [m] ? initialAlpha = 0.916290731874155 # [ ] ? stateModel = AgeingLaw +frictionModel = Truncated [boundary.friction.weakening] a = 0.015 # [1] ? b = 0.030 # [1] ? diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc index 8c562f87413ec94e0c2ab062d3d71c5470d083df..d4053ad24b8eb4f83b2ca4184d5069e4f2c82dc8 100644 --- a/src/sand-wedge.cc +++ b/src/sand-wedge.cc @@ -68,6 +68,7 @@ #include "assemblers.hh" #include "tobool.hh" #include "enum_parser.cc" +#include "enum_friction_model.cc" #include "enum_scheme.cc" #include "enum_state_model.cc" #include "enum_verbosity.cc" @@ -248,6 +249,7 @@ int main(int argc, char *argv[]) { body.getPoissonRatio(), u_initial); MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction")); auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity( + parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundary, frictionInfo, normalStress); myGlobalFriction->updateAlpha(alpha_initial);