#123 ported today's changes from #42 which fixed the problem

pou036 · May 8, 2017 · 238d00a · 238d00a
1 parent 5d5364c
commit 238d00a
Show file tree

Hide file tree

Showing 9 changed files with 37 additions and 19 deletions.
diff --git a/include/materials/RedbackMechMaterialDP.h b/include/materials/RedbackMechMaterialDP.h
@@ -35,7 +35,7 @@ class RedbackMechMaterialDP : public RedbackMechMaterial
   void getFlowTensor(const RankTwoTensor &, Real, Real, Real, RankTwoTensor &);
   Real getFlowIncrement(Real, Real, Real, Real, Real);
   void get_py_qy(Real, Real, Real &, Real &, Real);
-  Real getDerivativeFlowIncrement(const RankTwoTensor &, Real, Real, Real, Real);
+  Real getDerivativeFlowIncrement(const RankTwoTensor &, Real, Real, Real, Real, Real);
   Real getPressureProjection(Real, Real, Real);
 
   // virtual void form_damage_kernels(Real);

diff --git a/src/materials/RedbackMechMaterialCC.C b/src/materials/RedbackMechMaterialCC.C
@@ -64,8 +64,14 @@ RedbackMechMaterialCC::getFlowIncrement(Real sig_eqv, Real pressure, Real q_yiel
                                        /*x=*/pressure,
                                        /*y=*/sig_eqv))
   {
-    Real flow_incr_vol = _ref_pe_rate * _dt * std::pow(std::fabs(pressure - p_yield_stress), _exponent) * _exponential;
-    Real flow_incr_dev = _ref_pe_rate * _dt * std::pow(std::fabs(sig_eqv - q_yield_stress), _exponent) * _exponential;
+    Real flow_incr_vol =
+      _ref_pe_rate * _dt * std::pow(std::fabs((pressure - p_yield_stress) / pc), _exponent) * _exponential;
+    Real flow_incr_dev =
+      _ref_pe_rate * _dt * std::pow(std::fabs((sig_eqv - q_yield_stress) / pc), _exponent) * _exponential;
+    // Real flow_incr_vol = _ref_pe_rate * _dt * std::pow(std::fabs(pressure - p_yield_stress), _exponent) *
+    // _exponential;
+    // Real flow_incr_dev = _ref_pe_rate * _dt * std::pow(std::fabs(sig_eqv - q_yield_stress), _exponent) *
+    // _exponential;
     return std::pow(flow_incr_vol * flow_incr_vol + flow_incr_dev * flow_incr_dev, 0.5);
   }
   else
@@ -85,10 +91,14 @@ RedbackMechMaterialCC::getDerivativeFlowIncrement(
       _ref_pe_rate * _dt * std::pow(std::fabs(pressure - p_yield_stress), _exponent) * _exponential;
     Real delta_lambda_q = _ref_pe_rate * _dt * std::pow(std::fabs(sig_eqv - q_yield_stress), _exponent) * _exponential;
     Real delta_lambda = (std::pow(delta_lambda_p * delta_lambda_p + delta_lambda_q * delta_lambda_q, 0.5));
-    Real der_flow_incr_dev =
-      _ref_pe_rate * _dt * _exponent * std::pow(std::fabs(sig_eqv - q_yield_stress), _exponent - 1.0) * _exponential;
-    Real der_flow_incr_vol =
-      _ref_pe_rate * _dt * _exponent * std::pow(std::fabs(pressure - p_yield_stress), _exponent - 1.0) * _exponential;
+    Real der_flow_incr_dev = _ref_pe_rate * _dt * _exponent *
+                             std::pow(std::fabs((sig_eqv - q_yield_stress) / pc), _exponent - 1.0) * _exponential /
+                             std::fabs(pc);
+    // _ref_pe_rate * _dt * _exponent * std::pow(std::fabs(sig_eqv - q_yield_stress), _exponent - 1.0) * _exponential;
+    Real der_flow_incr_vol = _ref_pe_rate * _dt * _exponent *
+                             std::pow(std::fabs((pressure - p_yield_stress) / pc), _exponent - 1.0) * _exponential /
+                             std::fabs(pc);
+    // _ref_pe_rate * _dt * _exponent * std::pow(std::fabs(pressure - p_yield_stress), _exponent - 1.0) * _exponential;
     return (delta_lambda_q * der_flow_incr_dev + delta_lambda_p * der_flow_incr_vol) / delta_lambda;
   }
   else

diff --git a/src/materials/RedbackMechMaterialDP.C b/src/materials/RedbackMechMaterialDP.C
@@ -79,39 +79,47 @@ RedbackMechMaterialDP::getFlowTensor(
  */
 Real
 RedbackMechMaterialDP::getFlowIncrement(
-  Real sig_eqv, Real pressure, Real q_yield_stress, Real p_yield_stress, Real /*yield_stress*/)
+  Real sig_eqv, Real pressure, Real q_yield_stress, Real p_yield_stress, Real yield_stress)
 {
   Real flow_incr_vol =
     _ref_pe_rate * _dt *
-    std::pow(macaulayBracket((pressure - p_yield_stress) * (_slope_yield_surface < 0 ? 1 : -1)), _exponent) *
+    std::pow(macaulayBracket(((pressure - p_yield_stress) / yield_stress) * (_slope_yield_surface < 0 ? 1 : -1)),
+             _exponent) *
     _exponential;
   // TODO: q_yield_stress can be 0, we should handle that case properly...
   Real flow_incr_dev =
     _ref_pe_rate * _dt *
-    std::pow(macaulayBracket((q_yield_stress > 0 ? 1 : -1) * (sig_eqv / q_yield_stress - 1.0)), _exponent) *
+    std::pow(macaulayBracket((q_yield_stress > 0 ? 1 : -1) * ((sig_eqv - q_yield_stress) / yield_stress)), _exponent) *
     _exponential;
   //(q_yield_stress > 0 ? 1:-1) is the sign function
   return std::pow(flow_incr_vol * flow_incr_vol + flow_incr_dev * flow_incr_dev, 0.5);
   // TODO: change the formula to use dist_pq^m
 }
 
 Real
-RedbackMechMaterialDP::getDerivativeFlowIncrement(
-  const RankTwoTensor & /*sig*/, Real pressure, Real sig_eqv, Real q_yield_stress, Real p_yield_stress)
+RedbackMechMaterialDP::getDerivativeFlowIncrement(const RankTwoTensor & /*sig*/,
+                                                  Real pressure,
+                                                  Real sig_eqv,
+                                                  Real yield_stress,
+                                                  Real q_yield_stress,
+                                                  Real p_yield_stress)
 {
-  Real delta_lambda_p =
-    _ref_pe_rate * _dt * std::pow(macaulayBracket(pressure - p_yield_stress), _exponent) * _exponential;
+  Real delta_lambda_p = _ref_pe_rate * _dt *
+                        std::pow(macaulayBracket((pressure - p_yield_stress) / yield_stress), _exponent) *
+                        _exponential;
   Real delta_lambda_q =
     _ref_pe_rate * _dt *
-    std::pow(macaulayBracket((q_yield_stress > 0 ? 1 : -1) * (sig_eqv / q_yield_stress - 1.0)), _exponent) *
+    std::pow(macaulayBracket((q_yield_stress > 0 ? 1 : -1) * ((sig_eqv - q_yield_stress) / yield_stress)), _exponent) *
     _exponential;
   Real delta_lambda = (std::pow(delta_lambda_p * delta_lambda_p + delta_lambda_q * delta_lambda_q, 0.5));
   Real der_flow_incr_dev =
     _ref_pe_rate * _dt * _exponent *
-    std::pow(macaulayBracket((q_yield_stress > 0 ? 1 : -1) * (sig_eqv / q_yield_stress - 1.0)), _exponent - 1.0) *
-    _exponential / q_yield_stress;
+    std::pow(macaulayBracket((q_yield_stress > 0 ? 1 : -1) * ((sig_eqv - q_yield_stress) / yield_stress)),
+             _exponent - 1.0) *
+    _exponential / yield_stress;
   Real der_flow_incr_vol = _ref_pe_rate * _dt * _exponent *
-                           std::pow(macaulayBracket(pressure - p_yield_stress), _exponent - 1.0) * _exponential;
+                           std::pow(macaulayBracket((pressure - p_yield_stress) / yield_stress), _exponent - 1.0) *
+                           _exponential / yield_stress;
   return (delta_lambda_q * der_flow_incr_dev + delta_lambda_p * der_flow_incr_vol) / delta_lambda;
 }
 
@@ -135,7 +143,7 @@ RedbackMechMaterialDP::getJac(const RankTwoTensor & sig,
 
   sig_dev = sig.deviatoric();
 
-  dfi_dseqv = getDerivativeFlowIncrement(sig, pressure, sig_eqv, q_yield_stress, p_yield_stress);
+  dfi_dseqv = getDerivativeFlowIncrement(sig, pressure, sig_eqv, yield_stress, q_yield_stress, p_yield_stress);
   getFlowTensor(sig, sig_eqv, pressure, yield_stress, flow_dirn);
 
   /* The following calculates the tensorial derivative (Jacobian) of the

diff --git a/tests/benchmark_10_TMC/gold/bench_TMC_DP_out.e b/tests/benchmark_10_TMC/gold/bench_TMC_DP_out.e
diff --git a/tests/benchmark_11_THMC/gold/bench_THMC_DP_out.e b/tests/benchmark_11_THMC/gold/bench_THMC_DP_out.e
diff --git a/tests/benchmark_3_M/gold/bench_DP_out.e b/tests/benchmark_3_M/gold/bench_DP_out.e
diff --git a/tests/benchmark_6_TM/gold/bench_TM_DP_out.e b/tests/benchmark_6_TM/gold/bench_TM_DP_out.e
diff --git a/tests/benchmark_7_HM/gold/bench_HM_DP_out.e b/tests/benchmark_7_HM/gold/bench_HM_DP_out.e
diff --git a/tests/benchmark_8_THM/gold/bench_THM_DP_out.e b/tests/benchmark_8_THM/gold/bench_THM_DP_out.e