1@DSL Implicit; 2@Behaviour Hayhurst; 3@Author Jean-Michel Proix; 4@Date 07/02/2013; 5@Description{ 6 resolution implicite NEWTON classique 7} 8 9@Epsilon 1.e-15; 10@Theta 1; 11 12@MaterialProperty stress young; 13young.setGlossaryName("YoungModulus"); 14@MaterialProperty real nu; 15nu.setGlossaryName("PoissonRatio"); 16 17@MaterialProperty real K; 18@MaterialProperty real epsi0; 19@MaterialProperty real sigma0; 20@MaterialProperty real h1; 21@MaterialProperty real h2; 22@MaterialProperty real H1star; 23@MaterialProperty real H2star; 24@MaterialProperty real A0; 25@MaterialProperty real alphaD; 26@MaterialProperty real delta1; 27@MaterialProperty real delta2; 28 29@StateVariable real d; 30@StateVariable real p; 31@StateVariable real H1; 32@StateVariable real H2; 33 34@LocalVariable stress lambda; 35@LocalVariable stress mu; 36 37@Parameter dmax = 0.999999999999999; 38 39/* Initialize Lame coefficients */ 40@InitLocalVariables{ 41 lambda = computeLambda(young,nu); 42 mu = computeMu(young,nu); 43} 44 45@ComputeStress{ 46 if(d > dmax){ 47 sig = StressStensor(stress(0)); 48 } else { 49 sig = (1-d)*(lambda*trace(eel)*Stensor::Id()+2*mu*eel); 50 } 51} 52 53@Integrator{ 54 const real seq = sigmaeq(sig); 55 if(seq > 1e-8*young){ 56 const real d_ = min(d+theta*dd,dmax); 57 const real seq0 = seq/(1-d_); 58 const real H1_ = H1+theta*dH1; 59 const real H2_ = H2+theta*dH2; 60 const real H_ = H1_+H2_; 61 const real shp = sinh(seq*(1-H_)/K/(1-(d_))); 62 const real chp = sqrt(1+shp*shp) ; 63 const real trsig = max(trace(sig),stress(0)); 64 const real trsig0 = trsig/(1-d); 65 const real shd = sinh( (alphaD*trsig + (1-alphaD)*seq )/sigma0 ) ; 66 const real chd = sqrt(1+shd*shd) ; 67 const real iseq = 1/seq; 68 const real dtrsde=(3*lambda+2*mu)*theta*(1-d_)*trsig/trace(sig); 69 const Stensor n = 3*deviator(sig)*iseq/2; 70 Stensor dsvmde=2*mu*theta*(1-d_)*n; 71 Stensor dsvm0de=2*mu*theta*n; 72 // systeme a résoudre 73 feel += dp*n-deto; 74 fp = dp-epsi0*dt*shp; 75 fH1 = dH1-h1*dp*(H1star-delta1*H1_)*iseq ; 76 fH2 = dH2-h2*dp*(H2star-delta2*H2_)*iseq ; 77 fd = dd - A0 * shd*dt; 78 // jacobienne 79 dfeel_ddeel += 2*mu*theta*dp*iseq*(1-d_)*(Stensor4::M()-(n^n)); 80 dfeel_ddp = n ; 81 dfp_ddeel = - (dt*epsi0*chp*(1-H_)/K)*dsvm0de ; 82 dfp_ddH1 = dt*epsi0*chp*theta*seq0/K ; 83 dfp_ddH2 = dfp_ddH1 ; 84 dfH1_ddeel = h1*dp*(H1star - delta1*H1_)*iseq*iseq*dsvmde; 85 dfH2_ddeel = h2*dp*(H2star - delta2*H2_)*iseq*iseq*dsvmde; 86 dfH1_ddp = - h1*iseq*(H1star - delta1*H1_) ; 87 dfH2_ddp = - h2*iseq*(H2star - delta2*H2_) ; 88 dfH1_ddH1+= h1*delta1*dp*theta*iseq ; 89 dfH2_ddH2+= h2*delta2*dp*theta*iseq ; 90 dfH1_ddd= h1*dp*(H1star - delta1*H1_)*theta*seq0*iseq*iseq; 91 dfH2_ddd= h2*dp*(H2star - delta2*H2_)*theta*seq0*iseq*iseq; 92 dfd_ddd += A0*dt*chd/sigma0 *(1-alphaD)*theta*seq0 ; 93 if (trsig > 1e-8*young) { 94 dfd_ddeel = -A0*dt*chd/sigma0 *(alphaD*dtrsde*Stensor::Id()+ (1-alphaD)*dsvmde); 95 dfd_ddd += A0*dt*chd/sigma0 *alphaD*theta*trsig0 ; 96 } else { 97 dfd_ddeel = -A0*dt*chd/sigma0 *(1-alphaD)*dsvmde; 98 } 99 } else { 100 feel -= deto; 101 } 102} 103 104@TangentOperator{ 105 if((smt==ELASTIC)||(smt==SECANTOPERATOR)){ 106 Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id(); 107 } else if(smt==CONSISTENTTANGENTOPERATOR){ 108 const StiffnessTensor De = lambda*Stensor4::IxI()+2*mu*Stensor4::Id(); 109 Stensor4 Je; 110 Stensor Jd; 111 getPartialJacobianInvert(Je,Jd); 112 Dt = (1-d)*De*Je-((Jd)^(De*eel)); 113 } else { 114 return false; 115 } 116} 117