Merge pull request #15 from larsblatny/Associative_Viscoplasticity

Implemented Associative viscoplastic flow rule with MCC model

Author: larsblatny

src/simulation/plasticity.cpp

@@ -50,7 +50,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++;
@@ -298,7 +298,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;
@@ -315,16 +320,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);
@@ -510,10 +529,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;
@@ -543,98 +561,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
 
         // Get new stress
         TM Fe = particles.F[p];

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;