Hi guys, I am currently working on building my own model of radial return mapping algorithm for plasticity of a beam, but encountered some issue (which I don’t know why, the values just doesn’t match with the one I got from commercial software like Abaqus or Comsol)…… If any of you are willing to enlighten me / share about how your ‘successful code’ looks like, please do so……..
For helping you, we need a description (in terms of equations) of the problem you want to solve and/or a code try.
for(int step=1; step<=nSteps; ++step){
real alpha = real(step) / real(nSteps);
Wh epbarold = epbar; // i want to save the old epbar
Wh ep11old = ep11, ep22old = ep22, ep12old = ep12, ep33old = ep33;
Wh eel11old = eel11, eel22old = eel22, eel12old = eel12;
Wh e11old = e11, e22old = e22, e12old = e12; // old total strain… used for computing the trial elastic strain (incremental)
Qh s11old = s11, s22old = s22, s12old = s12, s33old = s33;
Wh de11old = de11, de22old = de22, de33old = de33, de12old = de12;
cout << "STEP " << alpha << ", epbar = " << epbar[].max << endl;
for(int it=0; it<maxIt; ++it){
u1old[] = u1[];
u2old[] = u2[]; // displacement solved by the solver
e11 = dx(u1); e22 = dy(u2); e12 = 0.5*(dy(u1)+dx(u2)); // tensorial shear (no sqrt2 here) engineering strain no 1/2
eel11 = e11 - ep11;
eel22 = e22 - ep22;
eel12 = e12 - ep12;
eel33 = - ep33; // e33 - ep33old , tr(ep) = 0 OR ep11old + ep22old
Wh deel11 = eel11 - eel11old;
Wh deel22 = eel22 - eel22old;
Wh deel12 = eel12 - eel12old;
Wh tre = e11 + e22; // (nice cancellation thanks to plane strain)
cout << "Checking if dgam is non zero = " << dgam << endl;
// epbar = epbarold + dgam; // initialise at first at the top loop? does it forget?
// deviatoric strain
de11 = eel11 - tre/3.; // trying purely elastic - assumed all increase is done through elastic strain
de22 = eel22 - tre/3.;
de33 = eel33 - tre/3.;
de12 = eel12; // of current
Wh dde11 = de11 - de11old; //computing the deviatoric strain increment
Wh dde22 = de22 - de22old;
Wh dde33 = de33 - de33old;
Wh dde12 = de12 - de12old;
// trial deviatoric stress (by increment)
s11tr = s11old + 2.*mu0*dde11;
s22tr = s22old + 2.*mu0*dde22;
s33tr = s33old + 2.*mu0*dde33;
s12tr = s12old + 2.*mu0*dde12;
Qh sds = s11tr*s11tr + s22tr*s22tr + s33tr*s33tr + 2.*s12tr*s12tr;
Qh q = sqrt(max(0., 1.5*sds));
sy = 0.0; // initialise
H = 0.0;
// ternary
// Yield stress sy (element-wise)
sy = syv[5]; // default to highest
sy = (epbar < epv[5]) ? (syv[4] + (syv[5] - syv[4]) * (epbar - epv[4]) / (epv[5] - epv[4])) : sy;
sy = (epbar < epv[4]) ? (syv[3] + (syv[4] - syv[3]) * (epbar - epv[3]) / (epv[4] - epv[3])) : sy;
sy = (epbar < epv[3]) ? (syv[2] + (syv[3] - syv[2]) * (epbar - epv[2]) / (epv[3] - epv[2])) : sy;
sy = (epbar < epv[2]) ? (syv[1] + (syv[2] - syv[1]) * (epbar - epv[1]) / (epv[2] - epv[1])) : sy;
sy = (epbar < epv[1]) ? (syv[0] + (syv[1] - syv[0]) * (epbar - epv[0]) / (epv[1] - epv[0])) : sy;
sy = (epbar <= epv[0]) ? syv[0] : sy;
// Hardening modulus H (element-wise)
H = 0.0; // default
H = (epbar < epv[5]) ? ((syv[5] - syv[4]) / (epv[5] - epv[4])) : H;
H = (epbar < epv[4]) ? ((syv[4] - syv[3]) / (epv[4] - epv[3])) : H;
H = (epbar < epv[3]) ? ((syv[3] - syv[2]) / (epv[3] - epv[2])) : H;
H = (epbar < epv[2]) ? ((syv[2] - syv[1]) / (epv[2] - epv[1])) : H;
H = (epbar < epv[1]) ? ((syv[1] - syv[0]) / (epv[1] - epv[0])) : H;
// checking yield function
fY = q - sy;
// NR loop for dgam.
dgam = 0; //init
for (int k=0; k<10; k++){
epbartrial = epbar + dgam; //epbarold here is updated with the previous dgam alr
sytrial = 0.0;
Htrial = 0.0;
// ternary operator test..
// Trial yield stress sytrial (element-wise)
sytrial = syv[5];
sytrial = (epbartrial < epv[5]) ? (syv[4] + (syv[5] - syv[4]) * (epbartrial - epv[4]) / (epv[5] - epv[4])) : sytrial;
sytrial = (epbartrial < epv[4]) ? (syv[3] + (syv[4] - syv[3]) * (epbartrial - epv[3]) / (epv[4] - epv[3])) : sytrial;
sytrial = (epbartrial < epv[3]) ? (syv[2] + (syv[3] - syv[2]) * (epbartrial - epv[2]) / (epv[3] - epv[2])) : sytrial;
sytrial = (epbartrial < epv[2]) ? (syv[1] + (syv[2] - syv[1]) * (epbartrial - epv[1]) / (epv[2] - epv[1])) : sytrial;
sytrial = (epbartrial < epv[1]) ? (syv[0] + (syv[1] - syv[0]) * (epbartrial - epv[0]) / (epv[1] - epv[0])) : sytrial;
sytrial = (epbartrial <= epv[0]) ? syv[0] : sytrial;
// Trial hardening modulus Htrial (element-wise)
Htrial = 0.0;
Htrial = (epbartrial < epv[5]) ? ((syv[5] - syv[4]) / (epv[5] - epv[4])) : Htrial;
Htrial = (epbartrial < epv[4]) ? ((syv[4] - syv[3]) / (epv[4] - epv[3])) : Htrial;
Htrial = (epbartrial < epv[3]) ? ((syv[3] - syv[2]) / (epv[3] - epv[2])) : Htrial;
Htrial = (epbartrial < epv[2]) ? ((syv[2] - syv[1]) / (epv[2] - epv[1])) : Htrial;
Htrial = (epbartrial < epv[1]) ? ((syv[1] - syv[0]) / (epv[1] - epv[0])) : Htrial;
dgamold = dgam;
Wh phi = q - 3.*mu0*dgam - sytrial;
Wh dphi = -3.*mu0 - Htrial;
dgam = dgamold - phi/dphi;
dgam = max(0.0, dgam);
Wh diff = dgam-dgamold;
dgamtest = ((q-sytrial)>0.) ? ( (q - sytrial) / (3.*mu0 + H) ) : 0.0; // for verifying purposes
real maxDdgam = diff[].max;
if (abs(maxDdgam) < 1e-10) {
syconv = sytrial;
cout << "MUST CHECK 2 sytrial= " << sytrial[].max << endl ;
cout << " NR converged in " << k+1 << " iterations" << endl ;
cout << "dgam NR = " << dgam[].max << endl ;
cout << "dgamtest = " << dgamtest[].max << endl ;
break; // dgam saved
}
} // epbar shld be saved here
dgamfin = dgam;
epbar = epbar + dgam;
//radial return
Wh rfac = 1.0 - (3.0 * mu0 * dgam) / (q + epsDiv);// is this a correct way
rfac = max(0.0, rfac);
sysurface = syconv; // updated with new epbar – should match the returned value
s11 = s11tr * rfac;
s22 = s22tr * rfac;
s33 = s33tr * rfac;
s12 = s12tr * rfac;
Qh sr = s11*s11 + s22*s22 + s33*s33 + 2.*s12*s12;
sreturn = sqrt(max(0., 1.5*sr)); // return should map the updated yield stress sysurface
Wh isPlastic = (dgam > 1e-12) ? 1.0 : 0.0;
Wh plmultiplier = 1.5*dgam / (q + epsDiv);
ep11 = ep11 + isPlastic * s11tr * plmultiplier;
ep22 = ep22 + isPlastic * s22tr * plmultiplier;
ep12 = ep12 + isPlastic * s12tr * plmultiplier;
ep33 = ep33 + isPlastic * s33tr * plmultiplier;
Wh GEff = rfac * mu0;
Wh LEff = (3.*K0 - 2.*GEff)/3.;
Wh muEff = min(GEff, mu0); // min. value
Wh lamEff = min(LEff, lambda0);
problem stepSolve([u1, u2], [du1, du2]) =
int2d(Th)(
lamEff * div(u1, u2) * div(du1, du2) +
2. * muEff * (epsilon(u1, u2)’ * epsilon(du1, du2))
)
- int1d(Th, 3)(f1 * du1 + (f2 * alpha) * du2) // External traction forces on boundary 3
+ on(1, u1 = 0, u2 = 0); // Boundary conditions (fixed displacement at boundary 1)
stepSolve;
here’s is the main code for loading loop and iteration loop as well as the newton raphson loop for getting plastic strain…… I have implemented almost 1:1 of what is written in the Abaqus manual except the consistent Jacobian matrix for solving…… I can’t figure out where when wrong, the results just don’t match… would appreciate it if there’s any suggestion..
So the main idea of what I am doing is elastoplasticity computation for linear isotropic hardening in 2D plane strain case. I am computing Von Mises Stress using the deviatoric stresses derived from deviatoric elastic strain and using return mapping method for plastic region…..
You should look at the examples of elasticity, either in the pdf documentation
or in the examples files at
(type “elast” in the “Type to start searching” input, with “Search in examples” on).