#123 refined the delta_p used for numerical derivative (tried to get …

…semi-analytical solution but the maths are really too ugly...). Note that too coarse delta_p can lead to non-convergence, but we can't use too fine a value either as "s" is itself computed numerically... That value of delta_p should ideally be dynamic (depending on the overstress)
pou036 · Apr 11, 2018 · ca8d838 · ca8d838
1 parent d4665f5
commit ca8d838
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 6 deletions.
diff --git a/src/materials/RedbackMechMaterialCC.C b/src/materials/RedbackMechMaterialCC.C
@@ -79,8 +79,8 @@ void
 RedbackMechMaterialCC::getDerivativeFlowIncrement(Real & dfi_dp,
                                                   Real & dfi_dq,
                                                   const RankTwoTensor & /*sig*/,
-                                                  Real pressure,
-                                                  Real sig_eqv,
+                                                  Real p,
+                                                  Real q,
                                                   Real pc,
                                                   Real /*q_yield_stress*/,
                                                   Real /*p_yield_stress*/,
@@ -97,11 +97,11 @@ RedbackMechMaterialCC::getDerivativeFlowIncrement(Real & dfi_dp,
   // Compute numerically derivatives of s with respect to p and q
   Real p_y2, q_y2, s2, ds_dp, ds_dq;
   bool is_plastic;
-  Real delta_p = sigma_0 / 100.;
+  Real delta_p = sigma_0 / 1000.; // TODO: this value influences the convergence and should be dynamic!
   Real delta_q = _M*delta_p;
-  get_py_qy(pressure + delta_p, sig_eqv, p_y2, q_y2, -pc, is_plastic, s2);
+  get_py_qy(p + delta_p, q, p_y2, q_y2, -pc, is_plastic, s2);
   ds_dp = (s2 - s) / delta_p;
-  get_py_qy(pressure, sig_eqv + delta_q, p_y2, q_y2, -pc, is_plastic, s2);
+  get_py_qy(p, q + delta_q, p_y2, q_y2, -pc, is_plastic, s2);
   ds_dq = (s2 - s) / delta_q;
   Real tmp = _ref_pe_rate * _dt * _exponent * _exponential * std::pow(s/sigma_0, _exponent - 1) / sigma_0;
   dfi_dp = tmp * ds_dp;

diff --git a/src/utils/Ellipse.C b/src/utils/Ellipse.C
@@ -231,7 +231,7 @@ Ellipse::getYieldPointCC(Real const m, Real const p_c, Real const y0, Real const
         (std::exp(-2*t)*(2*std::sqrt(1+alpha*std::exp(beta*2*t)) - beta*hyp2f1(0.5, -1/beta, (beta-1.0)/beta, -alpha*std::exp(beta*2*t))) \
          -2*std::sqrt(1+alpha) + beta*hyp2f1(0.5, -1/beta, (beta-1.0)/beta, -alpha)))
   */
-  int n_iter = 100;
+  int n_iter = 100; // TODO: is this value good enough even very far from the ellipse?
   s = 0;
   Real t_old = 0;
   Real t_new = 0;