1@Parser    Implicit;
2@Behaviour metalemani;
3@Includes{ #include<TFEL/Material/Hill.hxx>
4           #include<TFEL/Material/Lame.hxx>
5}
6@OrthotropicBehaviour;
7@Algorithm NewtonRaphson_NumericalJacobian;
8@Theta 1.; @Epsilon 1.e-10;
9@MaterialProperty real young;
10young.setGlossaryName("YoungModulus");
11@MaterialProperty real nu;
12nu.setGlossaryName("PoissonRatio");
13@MaterialProperty real a[3];
14@MaterialProperty real m[3];
15@MaterialProperty real pn[3];
16@MaterialProperty real Q[3];
17@MaterialProperty real M1[6];
18@MaterialProperty real M3[6];
19@StateVariable real p;
20@AuxiliaryStateVariable real seq;
21@AuxiliaryStateVariable real svi[3];
22
23@LocalVariable stress lambda;
24@LocalVariable stress mu;
25@LocalVariable tfel::math::st2tost2<N,real> H;
26@LocalVariable real T_ ;
27@LocalVariable real invn[3], f[3], gamma[3], sv[3]  ;
28// variables de commande aster
29@ExternalStateVariable real SECH,HYDR,IRRA,NEUT1,NEUT2,CORR,ALPHPUR,ALPHBETA;
30@IsTangentOperatorSymmetric true;
31@TangentOperator{using namespace tfel::material::lame;
32    StiffnessTensor Hooke;     Stensor4 Je;
33    computeElasticStiffness<N,Type>::exe(Hooke,lambda,mu);
34    getPartialJacobianInvert(Je);
35    Dt = Hooke*Je;  }
36@InitLocalVariables{
37  using namespace tfel::material::lame;
38  lambda = computeLambda(young,nu);
39  mu = computeMu(young,nu);
40  // proportion en phase alpha en milieu de pas de temps
41  const real Z = min(max(ALPHPUR  + theta*dALPHPUR+ ALPHBETA + theta*dALPHBETA,0.),1.);
42
43  // fonctions f
44  if        (Z >= 0.99) { f[0]=1. ;
45  } else if (Z >= 0.9)  { f[0] = (Z-0.9)/0.09 ;
46  } else                { f[0] = 0. ;  }
47
48  if        (Z >= 0.1)  { f[2]=0. ;
49  } else if (Z >= 0.01) { f[2] = (0.1-Z)/0.09 ;
50  } else                { f[2] = 1. ;  }
51
52  if        (Z >= 0.99) { f[1]=0. ;
53  } else if (Z >= 0.9)  { f[1] = 1.0-(Z-0.9)/0.09 ;
54  } else if (Z >= 0.1)  { f[1] = 1.0 ;
55  } else if (Z >= 0.01) { f[1] = 1.0-(0.1-Z)/0.09 ;
56  } else                { f[1] = 0. ;  }
57  // Temperature Aster en Celsius
58  T_ = 273.0 + T + theta * dT ;
59  for(unsigned short i=0;i!=3;++i){
60    invn[i] = 1.0 / pn[i] ;
61    gamma[i] = a[i]* exp(Q[i]/T_*invn[i]) ;  }
62
63// correspondance M aster (repere x,y,z) et H
64  real M[6];
65  if        (Z >= 0.99) {for(unsigned short i=0;i!=6;++i){M[i]=M1[i];}
66  } else if (Z >= 0.01) {for(unsigned short i=0;i!=6;++i){M[i]=Z*M1[i]+(1.-Z)*M3[i];}
67  } else                {for(unsigned short i=0;i!=6;++i){M[i]=M3[i];}}
68  const real H_F = 0.5*( M[0]+M[1]-M[2]);
69  const real H_G = 0.5*(-M[0]+M[1]+M[2]);
70  const real H_H = 0.5*( M[0]-M[1]+M[2]);
71  const real H_L = 2.0*M[3];
72  const real H_M = 2.0*M[4];
73  const real H_N = 2.0*M[5];
74  H = hillTensor<N,real>(H_F,H_G,H_H,H_L,H_M,H_N);
75}
76@ComputeStress{
77  sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;}
78@Integrator{
79  const real sigeq = sqrt(sig|H*sig);
80  real p_=p+theta*dp ;
81  real sigv =  0. ;  real pm[3];  real dpn[3] ;
82
83
84  for(unsigned short i=0;i!=3;++i){
85    pm[i]  = (p_ > 0.) ? pow(p_,m[i])         : 0.;
86    dpn[i] = (dp > 0.) ? pow((dp/dt),invn[i]) : 0. ;
87    sv[i]=gamma[i]*pm[i]*dpn[i] ;
88    sigv += f[i]*sv[i] ; }
89  Stensor  n(0.);
90  if(sigeq > 1.e-10*young){ n= (H*sig)/sigeq; }
91  feel += dp*n-deto;
92  fp    = (sigeq-sigv)/young;
93}
94@UpdateAuxiliaryStateVars{
95  for(unsigned short i=0;i!=3;++i){ svi[i]=sv[i] ; }
96}
97
98
99
100