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