1@DSL RungeKutta; 2@Algorithm rk54; 3@Behaviour PolyCrystals_DD_CC; 4 5@Author Jean-Michel Proix; 6@Date 08/10/2014; 7@Description{ 8 Polycristal BZ du modele DD_CC 9} 10@Epsilon 1.e-7 ; 11 12// Number of phases 13@IntegerConstant Np = 30 ; 14// Number of sliding systems 15@IntegerConstant Nss = 12; 16 17@ModellingHypothesis Tridimensional; 18@OrthotropicBehaviour; 19@RequireStiffnessTensor; 20 21@MaterialProperty real b; 22@MaterialProperty real H ; 23@MaterialProperty real DeltaG_0 ; 24@MaterialProperty real tau_0; 25@MaterialProperty real tau_f; 26@MaterialProperty real gamma0; 27@MaterialProperty real pn; 28@MaterialProperty real omega_mob; 29@MaterialProperty real d; 30@MaterialProperty real d_lath; 31@MaterialProperty real y_at; 32@MaterialProperty real K_f; 33@MaterialProperty real K_self ; 34@MaterialProperty real k_b; 35@MaterialProperty real epsi_1; 36@MaterialProperty real Mu; 37 38// pour la matrice d'interaction 39@MaterialProperty real h0; 40@MaterialProperty real h1; 41@MaterialProperty real h2; 42@MaterialProperty real h3; 43@MaterialProperty real h4; 44@MaterialProperty real h5; 45 46@StateVariable real pg; 47@StateVariable StrainStensor epsg[Np]; 48@StateVariable strain omega[Nss*Np]; 49@StateVariable StrainStensor epsp; 50omega.setErrorNormalisationFactor(1.e-7); 51 52// pour utiliser les sytemes de glissement CC 53// et aussi d'autres systemes de glissement (24) 54@Import "SlidingSystemsCC.mfront"; 55@Import "MonoCrystal_DD_CC_InteractionMatrix.mfront"; 56 57@Includes{ 58#include"TFEL/Material/PolyCrystalsSlidingSystems.hxx" 59} 60 61@Members{ 62 // Mono cristal sliding system 63 typedef SlidingSystems<real> MCSlidingSystems; 64 // Poly cristals sliding system 65 typedef PolyCrystalsSlidingSystems<Np,MCSlidingSystems,real> PCSlidingSystems; 66} 67 68//! fraction volumique 69@LocalVariable real fv[Np]; 70@LocalVariable real Lc ; 71 72@InitLocalVariables<Append>{ 73 // fractions volumiques 74 for(unsigned short i=0;i!=Np;++i){ 75 fv[i]=1.0/static_cast<real>(Np) ; 76 } 77 Lc = 500.*b*(T/300.)*(T/300.); 78} 79 80@ComputeStress{ 81 sig = D*eel; 82} 83 84@Derivative{ 85 const PCSlidingSystems& gs = 86 PCSlidingSystems::getPolyCrystalsSlidingSystems("PolyCrystalsAngles-30.txt"); 87 88 StressStensor sigg(real(0)) ; 89 // constantes monocristal 90 const real deuxpi = 6.2831853071795862; 91 real omegap[Nss]; 92 real small = 1.e-20 ; 93 94 depsp=Stensor(real(0)) ; 95 96 // boucle sur le nombre de phases (ou grains) 97 for(unsigned short k=0;k!=Np;++k){ 98 99 // localisation BZ 100 real seq=sigmaeq(sig); 101 if (seq>0){ 102 real alph=1./(1.0+1.5*Mu*pg/seq) ; 103 sigg = sig+ Mu*alph*(epsp-epsg[k]); 104 } 105 106 depsg[k]=Stensor(real(0)) ; 107 108 for(unsigned short s=0;s!=Nss;++s){ 109 omegap[s]=omega[Nss*k+s]; 110 if(omegap[s] < 0.) { 111 throw(runtime_error("omega négatif")); 112 } 113 } 114 for(unsigned short s=0;s!=Nss;++s){ 115 116 // tenseurs mus de chaque grain / systeme dans le repere global 117 const StrainStensor& mu_ki = gs.mus[k][s]; 118 const real tau = mu_ki | sigg ; 119 120 real omega_tot=0.; 121 for (unsigned short j=0;j!=Nss;++j){ 122 if (j!=s){ 123 omega_tot += omegap[j] ; 124 } 125 } 126 const real DG_app = k_b*T*log(omega_mob*H/sqrt(omega_tot)/epsi_1) ; 127 128 real t1 = 1.0-DG_app/DeltaG_0 ; 129 real Rs = 0. ; 130 if ( t1 <= small ) { 131 Rs = 1./small ; 132 } else { 133 Rs = Mu*b/(2.0*tau_0*t1*t1) ; 134 } 135 136 const real lambda = 1.0/ ( min ( sqrt(omega_tot)/b , omega_tot*(d+2.0*Rs)/b/b) ) -d ; 137 138 real alphat2 = 0.; 139 for (unsigned short j=0;j!=Nss;++j){ 140 if (j!=s){ 141 alphat2 += mh(s,j)*omegap[j] ; 142 } 143 } 144 145 if(alphat2 < 0.) { 146 throw(runtime_error("alphat négatif")); 147 } 148 const real alphat = sqrt(alphat2/omega_tot) ; 149 150 const real ls = max ( (lambda - 2.0 * alphat * Rs) , Lc ); 151 152 const real tau_LT = max ( 0.0, alphat*Mu*b* (1.0/lambda - 1.0/(2.0 * alphat * Rs + Lc)) ); 153 154 const real Homegas = sqrt (mh(s,s) * omegap[s] ); 155 156 const real tau_LR = Mu * Homegas ; 157 158 const real tau_c = tau_f + sqrt ( tau_LT*tau_LT + tau_LR*tau_LR ) ; 159 160 const real tau_eff = abs ( tau ) - tau_c ; 161 162 real sgn_tau = 0.0 ; 163 if (abs(tau) > 0. ) { 164 sgn_tau =tau/abs(tau) ; 165 } 166 167 real DG_eff = DeltaG_0 ; 168 169 if ((tau_eff > 0.0) && ( tau_eff < tau_0)) { 170 DG_eff = DG_eff * (1.0 - sqrt ( tau_eff / tau_0 ) ) ; 171 } 172 173 const real gamman = omega_mob / b* H * ls * exp ( - DG_eff / k_b / T ) * sgn_tau ; 174 175 if(abs(tau_c) < small ) { 176 ostringstream msg; 177 msg << "invalid value for tau_c (" << tau_c << ")."; 178 throw(runtime_error(msg.str())); 179 } 180 181 const real gammap = gamma0 * pow( ( abs(tau) / tau_c ), pn) * sgn_tau ; 182 183 real inv_gamman = 0. ; 184 if (abs(gamman) > 0. ) { 185 inv_gamman = 1.0 / gamman ; 186 } 187 188 real inv_gammap = 0. ; 189 if (abs(gammap) > 0. ) { 190 inv_gammap = 1.0 / gammap ; 191 } 192 real Dgamma =0. ; 193 if (abs( inv_gammap + inv_gamman) > 0. ) { 194 Dgamma = 1. / ( inv_gammap + inv_gamman ) ; 195 } 196 197 depsg[k] += Dgamma * mu_ki ; 198 199 const real ys = 1.0 / (1.0 / y_at + deuxpi * tau_eff / Mu / b) ; 200 201 real c_eff = 1.0 ; 202 203 if ((tau_eff > 0.0)&& (tau_0 > 0)){ 204 c_eff = 1.0 - tau_eff / tau_0 ; 205 } 206 207 const real Hs = b/d_lath + c_eff*( Homegas/K_self + alphat*lambda*omega_tot/K_f/b ) - ys*omegap[s]/b; 208 209 domega[Nss*k+s]= abs( Dgamma ) * Hs ; 210 211 } 212 depsp+=depsg[k]*fv[k] ; 213 } 214 dpg = sqrt(depsp|depsp)/sqrt(1.5); 215 deel = -depsp+deto; 216} 217 218 219 220 221