🔁 Merge pull request #705 from cgeudeker/material

Add shear and compression wave velocities to Linear Elastic Material Properties for #692
cb-geo · Mar 3, 2021 · 86ba10e · 86ba10e
2 parents 27086f7 + 30ef52e
commit 86ba10e
Show file tree

Hide file tree

Showing 3 changed files with 46 additions and 2 deletions.
diff --git a/include/materials/linear_elastic.h b/include/materials/linear_elastic.h
@@ -77,6 +77,10 @@ class LinearElastic : public Material<Tdim> {
   double poisson_ratio_{std::numeric_limits<double>::max()};
   //! Bulk modulus
   double bulk_modulus_{std::numeric_limits<double>::max()};
+  //! Compressional Wave Velocity
+  double vp_{std::numeric_limits<double>::max()};
+  //! Shear Wave Velocity
+  double vs_{std::numeric_limits<double>::max()};
 };  // LinearElastic class
 }  // namespace mpm
 

diff --git a/include/materials/linear_elastic.tcc b/include/materials/linear_elastic.tcc
@@ -9,10 +9,24 @@ mpm::LinearElastic<Tdim>::LinearElastic(unsigned id,
         material_properties.at("youngs_modulus").template get<double>();
     poisson_ratio_ =
         material_properties.at("poisson_ratio").template get<double>();
+
     // Calculate bulk modulus
     bulk_modulus_ = youngs_modulus_ / (3.0 * (1. - 2. * poisson_ratio_));
 
+    // Calculate constrained and shear modulus
+    double constrained_modulus =
+        youngs_modulus_ * (1. - poisson_ratio_) /
+        ((1. + poisson_ratio_) * (1. - 2. * poisson_ratio_));
+    double shear_modulus = youngs_modulus_ / (2.0 * (1. + poisson_ratio_));
+
+    // Calculate wave velocities
+    vp_ = sqrt(constrained_modulus / density_);
+    vs_ = sqrt(shear_modulus / density_);
+
     properties_ = material_properties;
+    properties_["pwave_velocity"] = vp_;
+    properties_["swave_velocity"] = vs_;
+
     // Set elastic tensor
     this->compute_elastic_tensor();
   } catch (Json::exception& except) {

diff --git a/tests/materials/linear_elastic_test.cc b/tests/materials/linear_elastic_test.cc
@@ -71,8 +71,6 @@ TEST_CASE("LinearElastic is checked in 2D", "[material][linear_elastic][2D]") {
             Approx(jmaterial["density"]).epsilon(Tolerance));
     REQUIRE(material->template property<double>("youngs_modulus") ==
             Approx(jmaterial["youngs_modulus"]).epsilon(Tolerance));
-    REQUIRE(material->template property<double>("poisson_ratio") ==
-            Approx(jmaterial["poisson_ratio"]).epsilon(Tolerance));
 
     // Check if state variable is initialised
     SECTION("State variable is initialised") {
@@ -147,6 +145,20 @@ TEST_CASE("LinearElastic is checked in 2D", "[material][linear_elastic][2D]") {
     REQUIRE(stress(4) == Approx(0.00000000000000e+00).epsilon(Tolerance));
     REQUIRE(stress(5) == Approx(0.00000000000000e+00).epsilon(Tolerance));
   }
+
+  SECTION("LinearElastic check properties earthquake") {
+    unsigned id = 0;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic2D", std::move(id), jmaterial);
+
+    // Get P-Wave and S-Wave Velocities
+    REQUIRE(material->template property<double>("pwave_velocity") ==
+            Approx(116.023870223).epsilon(Tolerance));
+    REQUIRE(material->template property<double>("swave_velocity") ==
+            Approx(62.0173672946).epsilon(Tolerance));
+  }
 }
 
 //! Check linearelastic class in 3D
@@ -289,4 +301,18 @@ TEST_CASE("LinearElastic is checked in 3D", "[material][linear_elastic][3D]") {
     REQUIRE(stress(4) == Approx(7.69230769230769e+01).epsilon(Tolerance));
     REQUIRE(stress(5) == Approx(1.15384615384615e+02).epsilon(Tolerance));
   }
+
+  SECTION("LinearElastic check properties earthquake") {
+    unsigned id = 0;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic3D", std::move(id), jmaterial);
+
+    // Get P-Wave and S-Wave Velocities
+    REQUIRE(material->template property<double>("pwave_velocity") ==
+            Approx(116.023870223).epsilon(Tolerance));
+    REQUIRE(material->template property<double>("swave_velocity") ==
+            Approx(62.0173672946).epsilon(Tolerance));
+  }
 }