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