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