Skip to content

Commit 4373c6e

Browse files
larsblatnylarsblatny
committed
added boolean option use_duvaut_lion_formulation, merged MCCVisc and MCCVisc2
1 parent eaa3b14 commit 4373c6e

File tree

3 files changed

+109
-131
lines changed

3 files changed

+109
-131
lines changed

‎src/simulation/plasticity.cpp‎

Lines changed: 107 additions & 130 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
1313
// Do nothing
1414
}
1515

16-
else if (plastic_model == PlasticModel::VM || plastic_model == PlasticModel::DP || plastic_model == PlasticModel::DPSoft || plastic_model == PlasticModel::MCC || plastic_model == PlasticModel::VMVisc || plastic_model == PlasticModel::DPVisc || plastic_model == PlasticModel::MCCVisc || plastic_model == PlasticModel::DPMui || plastic_model == PlasticModel::MCCMui || plastic_model == PlasticModel::MCCVisc2){
16+
else if (plastic_model == PlasticModel::VM || plastic_model == PlasticModel::DP || plastic_model == PlasticModel::DPSoft || plastic_model == PlasticModel::MCC || plastic_model == PlasticModel::VMVisc || plastic_model == PlasticModel::DPVisc || plastic_model == PlasticModel::MCCVisc || plastic_model == PlasticModel::DPMui || plastic_model == PlasticModel::MCCMui){
1717

1818
Eigen::JacobiSVD<TM> svd(Fe_trial, Eigen::ComputeFullU | Eigen::ComputeFullV);
1919
// TV hencky = svd.singularValues().array().log();
@@ -281,7 +281,12 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
281281
T delta_gamma;
282282

283283
if (perzyna_exp == 1){
284-
delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt);
284+
if (use_duvaut_lions_formulation){
285+
delta_gamma = (q_trial-q_yield) / (f_mu_prefac*(1+perzyna_visc/dt));
286+
}
287+
else{
288+
delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt);
289+
}
285290
if (delta_gamma < 0){
286291
debug("DPVisc: FATAL negative delta_gamma = ", delta_gamma);
287292
exit = 1;
@@ -298,16 +303,30 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
298303
exit = 1;
299304
}
300305

301-
T tm = perzyna_visc * delta_gamma + dt;
302-
T tmp = dt / tm;
303-
T tmp1 = std::pow(tmp, perzyna_exp);
306+
T residual;
307+
T residual_diff;
308+
if (use_duvaut_lions_formulation){
309+
310+
residual = q_trial - q_yield - f_mu_prefac * ( std::pow(perzyna_visc/dt*delta_gamma, perzyna_exp) + delta_gamma );
311+
if (std::abs(residual) < 1e-2) {
312+
break;
313+
}
314+
residual_diff = -f_mu_prefac * ( std::pow(perzyna_visc/dt, perzyna_exp) * perzyna_exp * std::pow(delta_gamma, perzyna_exp-1) + 1 );
304315

305-
T residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield;
306-
if (std::abs(residual) < 1e-2) {
307-
break;
308316
}
317+
else{
309318

310-
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);
319+
T tm = perzyna_visc * delta_gamma + dt;
320+
T tmp = dt / tm;
321+
T tmp1 = std::pow(tmp, perzyna_exp);
322+
323+
residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield;
324+
if (std::abs(residual) < 1e-2) {
325+
break;
326+
}
327+
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);
328+
329+
}
311330

312331
if (std::abs(residual_diff) < 1e-14){ // otherwise division by zero
313332
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
493512
T p_trial = p_stress;
494513
T q_trial = q_stress;
495514

496-
if (plastic_model == PlasticModel::MCCVisc){
497-
if (use_mibf)
498-
particles.muI[p] = M;
499-
}
515+
if (use_mibf)
516+
particles.muI[p] = M;
517+
500518

501519
bool perform_rma;
502520
T p_c;
@@ -526,142 +544,81 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
526544
exit = 1;
527545
}
528546

529-
530547
if (perform_rma) { // returns true if it performs a return mapping
531548
plastic_count++;
532549

533-
T ep = p_stress / (K*dim);
534-
T eps_pl_vol_inst = hencky_trace + dim * ep;
535-
particles.eps_pl_vol[p] += eps_pl_vol_inst;
550+
T b = 1; // must be equal to one unless using duvaut lion formulation
536551

537-
T delta_gamma;
538552
////////////////////////////////////////////////////////////////
539553
if (plastic_model == PlasticModel::MCCVisc){
540-
541-
T q_yield = q_stress;
542-
543-
if (perzyna_exp == 1){
544-
delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt);
545-
if (delta_gamma < 0){
546-
if (delta_gamma > -1e-14){
547-
delta_gamma = 0;
548-
}
549-
else{
550-
debug("MCCVisc: FATAL negative delta_gamma = ", delta_gamma);
551-
exit = 1;
552-
}
554+
if (use_duvaut_lions_formulation){
555+
b = 1.0 / (1.0 + perzyna_visc / dt); // b defined such that "stress = (1-b) * strain_trial + b * stress_yield"
556+
if (perzyna_exp != 1){
557+
debug("MCCVisc: FATAL viscous exponent must be 1 when use_duvaut_lions_formulation = true");
558+
exit = 1;
553559
}
554560
}
555-
else{ // persyna_exp is not one
556-
557-
delta_gamma = 0.01 * (q_trial - q_yield) / f_mu_prefac; // initial guess
558-
559-
int max_iter = 60;
560-
for (int iter = 0; iter < max_iter; iter++) {
561-
if (iter == max_iter - 1){ // did not break loop
562-
debug("MCCVisc: FATAL did not exit loop at iter = ", iter);
563-
exit = 1;
561+
else{ // perzyna
562+
T delta_gamma;
563+
T q_yield = q_stress;
564+
565+
if (perzyna_exp == 1){
566+
delta_gamma = (q_trial-q_yield) / (f_mu_prefac + q_yield*perzyna_visc/dt);
567+
if (delta_gamma < 0){
568+
if (delta_gamma > -1e-14){
569+
delta_gamma = 0;
570+
}
571+
else{
572+
debug("MCCVisc: FATAL negative delta_gamma = ", delta_gamma);
573+
exit = 1;
574+
}
564575
}
576+
}
577+
else{ // perzyna_exp is not one
565578

566-
T tm = perzyna_visc * delta_gamma + dt;
567-
T tmp = dt / tm;
568-
T tmp1 = std::pow(tmp, perzyna_exp);
569-
570-
T residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield;
571-
if (std::abs(residual) < 1e-2) {
572-
break;
573-
}
574-
575-
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);
576-
577-
if (std::abs(residual_diff) < 1e-14){ // otherwise division by zero
578-
debug("MCCVisc: FATAL residual_diff too small in abs value = ", residual_diff);
579-
exit = 1;
580-
}
581-
582-
delta_gamma -= residual / residual_diff;
583-
584-
if (delta_gamma < 0) // not possible and can also lead to division by zero
585-
delta_gamma = 1e-10;
579+
delta_gamma = 0.01 * (q_trial - q_yield) / f_mu_prefac; // initial guess
586580

587-
} // end N-R iterations
581+
int max_iter = 60;
582+
for (int iter = 0; iter < max_iter; iter++) {
583+
if (iter == max_iter - 1){ // did not break loop
584+
debug("MCCVisc: FATAL did not exit loop at iter = ", iter);
585+
exit = 1;
586+
}
588587

589-
} // end if perzyna_exp == 1
588+
T tm = perzyna_visc * delta_gamma + dt;
589+
T tmp = dt / tm;
590+
T tmp1 = std::pow(tmp, perzyna_exp);
590591

591-
q_stress = std::max(q_yield, q_trial - f_mu_prefac * delta_gamma); // delta_gamma = dt * gamma_dot_S
592+
T residual = (q_trial - f_mu_prefac * delta_gamma) * tmp1 - q_yield;
593+
if (std::abs(residual) < 1e-2) {
594+
break;
595+
}
592596

593-
if (use_mibf){
594-
if (hardening_law == HardeningLaw::ExpoImpl){
595-
p_c = std::max(T(0), p0 * std::exp(-xi*particles.eps_pl_vol[p]));
596-
}
597-
else if (hardening_law == HardeningLaw::SinhImpl){
598-
p_c = std::max(T(0), K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K)));
599-
}
600-
601-
if ( (p_stress < -beta * p_c + stress_tolerance) || (p_stress > p_c - stress_tolerance) ){
602-
particles.muI[p] = M;
603-
}
604-
else{
605-
particles.muI[p] = q_stress / std::sqrt((p_stress + beta * p_c) * (p_c - p_stress));
606-
}
607-
}
597+
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);
608598

609-
} // end if MCCVisc
610-
////////////////////////////////////////////////////////////////
599+
if (std::abs(residual_diff) < 1e-14){ // otherwise division by zero
600+
debug("MCCVisc: FATAL residual_diff too small in abs value = ", residual_diff);
601+
exit = 1;
602+
}
611603

604+
delta_gamma -= residual / residual_diff;
612605

613-
delta_gamma = (q_trial - q_stress) / f_mu_prefac;
614-
particles.eps_pl_dev[p] += (1.0/d_prefac) * delta_gamma;
615-
particles.delta_gamma[p] = delta_gamma / dt;
606+
if (delta_gamma < 0) // not possible and can also lead to division by zero
607+
delta_gamma = 1e-10;
616608

617-
hencky = q_stress / e_mu_prefac * hencky_deviatoric - ep*TV::Ones();
618-
particles.F[p] = svd.matrixU() * hencky.array().exp().matrix().asDiagonal() * svd.matrixV().transpose();
619-
} // if perform_rma
620-
} // end MCC / MCCVisc
621-
622-
else if (plastic_model == PlasticModel::MCCVisc2){
609+
} // end N-R iterations
623610

624-
// the trial stress states
625-
T p_stress = -K * hencky_trace;
626-
T q_stress = e_mu_prefac * hencky_deviatoric_norm;
611+
} // end if perzyna_exp == 1
627612

628-
// make copies
629-
T p_trial = p_stress;
630-
T q_trial = q_stress;
613+
q_stress = std::max(q_yield, q_trial - f_mu_prefac * delta_gamma); // delta_gamma = dt * gamma_dot_S
631614

632-
bool perform_rma;
633-
T p_c;
634-
if (hardening_law == HardeningLaw::NoHard){ // Exponential Explicit Hardening
635-
perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p0, beta, mu, K, f_mu_prefac);
636-
// perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K);
637-
}
638-
else if (hardening_law == HardeningLaw::ExpoExpl){ // Exponential Explicit Hardening
639-
p_c = std::max(stress_tolerance, p0*std::exp(-xi*particles.eps_pl_vol[p]));
640-
perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, f_mu_prefac);
641-
// perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K);
642-
}
643-
else if (hardening_law == HardeningLaw::SinhExpl){ // Sinh Explicit Hardening
644-
p_c = std::max(stress_tolerance, K*std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K)));
645-
perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, f_mu_prefac);
646-
// perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K);
647-
}
648-
else if (hardening_law == HardeningLaw::ExpoImpl){ // Exponential Implicit Hardening
649-
perform_rma = MCCRMAImplicitExponentialOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, particles.eps_pl_vol[p]);
650-
// perform_rma = MCCRMAImplicitExponential(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, f_mu_prefac, particles.eps_pl_vol[p]);
651-
}
652-
else if (hardening_law == HardeningLaw::SinhImpl) {
653-
perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, particles.eps_pl_vol[p]);
654-
}
655-
else{
656-
debug("You specified an invalid HARDENING LAW!");
657-
exit = 1;
658-
}
615+
} // end if use_duvaut_lions_formulation
659616

660-
if (perform_rma) { // returns true if it performs a return mapping
617+
} // end if MCCVisc
618+
////////////////////////////////////////////////////////////////
661619

662-
T tmp = 1.0 / (1.0 + perzyna_visc / dt);
663-
T eps_pl_vol_inst = - tmp * (p_trial - p_stress) / K;
664-
T eps_pl_dev_inst = tmp * (q_trial - q_stress) / e_mu_prefac;
620+
T eps_pl_vol_inst = - b * (p_trial - p_stress) / K;
621+
T eps_pl_dev_inst = b * (q_trial - q_stress) / e_mu_prefac;
665622

666623
particles.eps_pl_vol[p] += eps_pl_vol_inst;
667624
particles.eps_pl_dev[p] += eps_pl_dev_inst;
@@ -672,10 +629,30 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
672629
T h_dev = q_trial / e_mu_prefac - eps_pl_dev_inst;
673630
hencky = h_dev * hencky_deviatoric + (h_vol/dim) * TV::Ones();
674631
particles.F[p] = svd.matrixU() * hencky.array().exp().matrix().asDiagonal() * svd.matrixV().transpose();
675-
676-
} // if perform_rma
677-
678-
} // end if MCCVisc2
632+
633+
if (use_mibf){
634+
if (hardening_law == HardeningLaw::ExpoImpl){
635+
p_c = std::max(T(0), p0 * std::exp(-xi*particles.eps_pl_vol[p]));
636+
}
637+
else if (hardening_law == HardeningLaw::SinhImpl){
638+
p_c = std::max(T(0), K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K)));
639+
}
640+
641+
if ( (p_stress < -beta * p_c + stress_tolerance) || (p_stress > p_c - stress_tolerance) ){
642+
particles.muI[p] = M;
643+
}
644+
else{
645+
if (use_duvaut_lions_formulation){
646+
particles.muI[p] = ( (e_mu_prefac*h_dev-q_trial*(1-b))/b ) / std::sqrt((p_stress + beta * p_c) * (p_c - p_stress));
647+
}
648+
else{
649+
particles.muI[p] = q_stress / std::sqrt((p_stress + beta * p_c) * (p_c - p_stress));
650+
}
651+
}
652+
} // end if use_mibf
653+
654+
} // end if perform_rma
655+
} // end MCC or MCCVisc
679656

680657
} // end plastic_model type
681658

‎src/simulation/simulation.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,7 @@ class Simulation{
9393
// Perzyna
9494
T perzyna_exp = 1;
9595
T perzyna_visc = 0;
96+
bool use_duvaut_lions_formulation = false;
9697

9798
// MCC
9899
T beta = 0;

‎src/tools.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@ typedef double T; // float or double
5252

5353
enum class PlateType { top, bottom, left, right, front, back };
5454
enum class ElasticModel { Hencky, NeoHookean };
55-
enum class PlasticModel { NoPlasticity, VM, DP, DPSoft, MCC, VMVisc, DPVisc, MCCVisc, MCCVisc2, DPMui, MCCMui};
55+
enum class PlasticModel { NoPlasticity, VM, DP, DPSoft, MCC, VMVisc, DPVisc, MCCVisc, DPMui, MCCMui};
5656
enum class HardeningLaw { NoHard, ExpoExpl, ExpoImpl, SinhExpl, SinhImpl };
5757
enum class BC { NoSlip, SlipStick, SlipFree };
5858

0 commit comments

Comments
 (0)