search

luchete80 / WeldForm

Weakly-Compressible Smoothed Particle Hydrodynamics Parallel Solver for Elasto Plastic and thermal coupled Solid Mechanics
GNU General Public License v3.0
9 stars 5 forks source link

Plastic Strain in greater than Absolute in Verlet and ModVerlet #84

Open luchete80 opened 3 years ago

luchete80 commented 3 years ago

image

image

luchete80 commented 3 years ago 
    if (Fail == 1) {
        double J2   = 0.5*(ShearStress(0,0)*ShearStress(0,0) + 2.0*ShearStress(0,1)*ShearStress(1,0) +
                        2.0*ShearStress(0,2)*ShearStress(2,0) + ShearStress(1,1)*ShearStress(1,1) +
                        2.0*ShearStress(1,2)*ShearStress(2,1) + ShearStress(2,2)*ShearStress(2,2));
        //Scale back, Fraser Eqn 3-53
        double sig_trial = sqrt(3.0*J2);
        ShearStress = std::min((Sigmay/sig_trial),1.0)*ShearStress;
        if ( sig_trial > Sigmay) {
            double dep=( sig_trial - Sigmay)/ (3.*G);   //Fraser, Eq 3-49 TODO: MODIFY FOR TANGENT MODULUS = 0
            pl_strain+=dep;
        }
    }
luchete80 commented 3 years ago

image

image

luchete80 commented 3 years ago

ORIGINALLY IT WAS:

    // Elastic prediction step (ShearStress_e n+1)
    Stress          = ShearStress;
    if (ct == 30)
        ShearStress = dt*(2.0*G*(StrainRate-1.0/3.0*(StrainRate(0,0)+StrainRate(1,1)+StrainRate(2,2))*OrthoSys::I)+SRT+RS) + ShearStress;
    else
        ShearStress = 2.0*dt*(2.0*G*(StrainRate-1.0/3.0*(StrainRate(0,0)+StrainRate(1,1)+StrainRate(2,2))*OrthoSys::I)+SRT+RS) + ShearStressb;
    ShearStressb    = Stress;

    if (Fail == 1)
    {
        double J2   = 0.5*(ShearStress(0,0)*ShearStress(0,0) + 2.0*ShearStress(0,1)*ShearStress(1,0) +
                        2.0*ShearStress(0,2)*ShearStress(2,0) + ShearStress(1,1)*ShearStress(1,1) +
                        2.0*ShearStress(1,2)*ShearStress(2,1) + ShearStress(2,2)*ShearStress(2,2));
        //Scale back
        ShearStress = std::min((Sigmay/sqrt(3.0*J2)),1.0)*ShearStress;
    }

    Sigma           = -Pressure * OrthoSys::I + ShearStress;

    Stress  = Strain;
    if (ct == 30)
        Strain  = dt*StrainRate + Strain;
    else
        Strain  = 2.0*dt*StrainRate + Strainb;
    Strainb = Stress;