diff --git a/src/simulation/plasticity.cpp b/src/simulation/plasticity.cpp index b207156..89e019e 100644 --- a/src/simulation/plasticity.cpp +++ b/src/simulation/plasticity.cpp @@ -33,7 +33,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F // If hardening law: // q_yield = f(q_max, hardening variable) .... - T delta_gamma = hencky_deviatoric_norm - q_yield / e_mu_prefac; // this is eps_pl_dev_instant + T delta_gamma = hencky_deviatoric_norm - q_yield / e_mu_prefac; // this is eps_pl_dev_inst if (delta_gamma > 0){ // project to yield surface plastic_count++; @@ -281,7 +281,12 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F T delta_gamma; if (perzyna_exp == 1){ - delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt); + if (use_duvaut_lions_formulation){ + delta_gamma = (q_trial-q_yield) / (f_mu_prefac*(1+perzyna_visc/dt)); + } + else{ + delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt); + } if (delta_gamma < 0){ debug("DPVisc: FATAL negative delta_gamma = ", delta_gamma); exit = 1; @@ -298,16 +303,30 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F exit = 1; } - T tm = perzyna_visc * delta_gamma + dt; - T tmp = dt / tm; - T tmp1 = std::pow(tmp, perzyna_exp); + T residual; + T residual_diff; + if (use_duvaut_lions_formulation){ + + residual = q_trial - q_yield - f_mu_prefac * ( std::pow(perzyna_visc/dt*delta_gamma, perzyna_exp) + delta_gamma ); + if (std::abs(residual) < 1e-2) { + break; + } + residual_diff = -f_mu_prefac * ( std::pow(perzyna_visc/dt, perzyna_exp) * perzyna_exp * std::pow(delta_gamma, perzyna_exp-1) + 1 ); - T residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield; - if (std::abs(residual) < 1e-2) { - break; } + else{ - T residual_diff = -f_mu_prefac * tmp1 + (q_trial - f_mu_prefac * delta_gamma) * perzyna_exp * std::pow(tmp, perzyna_exp - 1) * (-perzyna_visc * dt) / (tm * tm); + T tm = perzyna_visc * delta_gamma + dt; + T tmp = dt / tm; + T tmp1 = std::pow(tmp, perzyna_exp); + + residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield; + if (std::abs(residual) < 1e-2) { + break; + } + residual_diff = -f_mu_prefac * tmp1 + (q_trial - f_mu_prefac * delta_gamma) * perzyna_exp * std::pow(tmp, perzyna_exp - 1) * (-perzyna_visc * dt) / (tm * tm); + + } if (std::abs(residual_diff) < 1e-14){ // otherwise division by zero debug("DPVisc: FATAL residual_diff too small in abs value = ", residual_diff); @@ -493,10 +512,9 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F T p_trial = p_stress; T q_trial = q_stress; - if (plastic_model == PlasticModel::MCCVisc){ - if (use_mibf) - particles.muI[p] = M; - } + if (use_mibf) + particles.muI[p] = M; + bool perform_rma; T p_c; @@ -526,98 +544,115 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F exit = 1; } - if (perform_rma) { // returns true if it performs a return mapping plastic_count++; - T ep = p_stress / (K*dim); - T eps_pl_vol_inst = hencky_trace + dim * ep; - particles.eps_pl_vol[p] += eps_pl_vol_inst; + T b = 1; // must be equal to one unless using duvaut lion formulation - T delta_gamma; //////////////////////////////////////////////////////////////// if (plastic_model == PlasticModel::MCCVisc){ - - T q_yield = q_stress; - - if (perzyna_exp == 1){ - delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt); - if (delta_gamma < 0){ - if (delta_gamma > -1e-14){ - delta_gamma = 0; - } - else{ - debug("MCCVisc: FATAL negative delta_gamma = ", delta_gamma); - exit = 1; - } + if (use_duvaut_lions_formulation){ + b = 1.0 / (1.0 + perzyna_visc / dt); // b defined such that "stress = (1-b) * strain_trial + b * stress_yield" + if (perzyna_exp != 1){ + debug("MCCVisc: FATAL viscous exponent must be 1 when use_duvaut_lions_formulation = true"); + exit = 1; } } - else{ // persyna_exp is not one + else{ // perzyna + T delta_gamma; + T q_yield = q_stress; + + if (perzyna_exp == 1){ + delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt); + if (delta_gamma < 0){ + if (delta_gamma > -1e-14){ + delta_gamma = 0; + } + else{ + debug("MCCVisc: FATAL negative delta_gamma = ", delta_gamma); + exit = 1; + } + } + } + else{ // perzyna_exp is not one - delta_gamma = 0.01 * (q_trial - q_yield) / f_mu_prefac; // initial guess + delta_gamma = 0.01 * (q_trial - q_yield) / f_mu_prefac; // initial guess - int max_iter = 60; - for (int iter = 0; iter < max_iter; iter++) { - if (iter == max_iter - 1){ // did not break loop - debug("MCCVisc: FATAL did not exit loop at iter = ", iter); - exit = 1; - } + int max_iter = 60; + for (int iter = 0; iter < max_iter; iter++) { + if (iter == max_iter - 1){ // did not break loop + debug("MCCVisc: FATAL did not exit loop at iter = ", iter); + exit = 1; + } - T tm = perzyna_visc * delta_gamma + dt; - T tmp = dt / tm; - T tmp1 = std::pow(tmp, perzyna_exp); + T tm = perzyna_visc * delta_gamma + dt; + T tmp = dt / tm; + T tmp1 = std::pow(tmp, perzyna_exp); - T residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield; - if (std::abs(residual) < 1e-2) { - break; - } + T residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield; + if (std::abs(residual) < 1e-2) { + break; + } - T residual_diff = -f_mu_prefac * tmp1 + (q_trial - f_mu_prefac * delta_gamma) * perzyna_exp * std::pow(tmp, perzyna_exp - 1) * (-perzyna_visc * dt) / (tm * tm); + T residual_diff = -f_mu_prefac * tmp1 + (q_trial - f_mu_prefac * delta_gamma) * perzyna_exp * std::pow(tmp, perzyna_exp - 1) * (-perzyna_visc * dt) / (tm * tm); - if (std::abs(residual_diff) < 1e-14){ // otherwise division by zero - debug("MCCVisc: FATAL residual_diff too small in abs value = ", residual_diff); - exit = 1; - } + if (std::abs(residual_diff) < 1e-14){ // otherwise division by zero + debug("MCCVisc: FATAL residual_diff too small in abs value = ", residual_diff); + exit = 1; + } - delta_gamma -= residual / residual_diff; + delta_gamma -= residual / residual_diff; - if (delta_gamma < 0) // not possible and can also lead to division by zero - delta_gamma = 1e-10; + if (delta_gamma < 0) // not possible and can also lead to division by zero + delta_gamma = 1e-10; - } // end N-R iterations + } // end N-R iterations - } // end if perzyna_exp == 1 + } // end if perzyna_exp == 1 - q_stress = std::max(q_yield, q_trial - f_mu_prefac * delta_gamma); // delta_gamma = dt * gamma_dot_S + q_stress = std::max(q_yield, q_trial - f_mu_prefac * delta_gamma); // delta_gamma = dt * gamma_dot_S - if (use_mibf){ - if (hardening_law == HardeningLaw::ExpoImpl){ - p_c = std::max(T(0), p0 * std::exp(-xi*particles.eps_pl_vol[p])); - } - else if (hardening_law == HardeningLaw::SinhImpl){ - p_c = std::max(T(0), K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K))); - } + } // end if use_duvaut_lions_formulation + + } // end if MCCVisc + //////////////////////////////////////////////////////////////// + + T eps_pl_vol_inst = - b * (p_trial - p_stress) / K; + T eps_pl_dev_inst = b * (q_trial - q_stress) / e_mu_prefac; + + particles.eps_pl_vol[p] += eps_pl_vol_inst; + particles.eps_pl_dev[p] += eps_pl_dev_inst; - if ( (p_stress < -beta * p_c + stress_tolerance) || (p_stress > p_c - stress_tolerance) ){ - particles.muI[p] = M; + particles.delta_gamma[p] = d_prefac * eps_pl_dev_inst / dt; + + T h_vol = - p_trial / K - eps_pl_vol_inst; + T h_dev = q_trial / e_mu_prefac - eps_pl_dev_inst; + hencky = h_dev * hencky_deviatoric + (h_vol/dim) * TV::Ones(); + particles.F[p] = svd.matrixU() * hencky.array().exp().matrix().asDiagonal() * svd.matrixV().transpose(); + + if (use_mibf){ + if (hardening_law == HardeningLaw::ExpoImpl){ + p_c = std::max(T(0), p0 * std::exp(-xi*particles.eps_pl_vol[p])); + } + else if (hardening_law == HardeningLaw::SinhImpl){ + p_c = std::max(T(0), K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K))); + } + + if ( (p_stress < -beta * p_c + stress_tolerance) || (p_stress > p_c - stress_tolerance) ){ + particles.muI[p] = M; + } + else{ + if (use_duvaut_lions_formulation){ + particles.muI[p] = ( (e_mu_prefac*h_dev-q_trial*(1-b))/b ) / std::sqrt((p_stress + beta * p_c) * (p_c - p_stress)); } else{ particles.muI[p] = q_stress / std::sqrt((p_stress + beta * p_c) * (p_c - p_stress)); } } + } // end if use_mibf - } // end if MCCVisc - //////////////////////////////////////////////////////////////// - - - delta_gamma = (q_trial - q_stress) / f_mu_prefac; - particles.eps_pl_dev[p] += (1.0/d_prefac) * delta_gamma; - particles.delta_gamma[p] = delta_gamma / dt; - - hencky = q_stress / e_mu_prefac * hencky_deviatoric - ep*TV::Ones(); - particles.F[p] = svd.matrixU() * hencky.array().exp().matrix().asDiagonal() * svd.matrixV().transpose(); - } // if perform_rma - } // end MCC / MCCVisc + } // end if perform_rma + } // end MCC or MCCVisc } // end plastic_model type diff --git a/src/simulation/simulation.hpp b/src/simulation/simulation.hpp index 0f0ee0e..25e4f5c 100644 --- a/src/simulation/simulation.hpp +++ b/src/simulation/simulation.hpp @@ -93,6 +93,7 @@ class Simulation{ // Perzyna T perzyna_exp = 1; T perzyna_visc = 0; + bool use_duvaut_lions_formulation = false; // MCC T beta = 0;