Synthèse des résultats ====================== Pour les modélisations A à E, on obtient un écart entre la solution SCILAB et la solution *Code_Aster* inférieur à :math:`\text{1.5 \%}` ce qui permet de valider l’implémentation des différentes lois de séchage dans le code. Notons simplement qu’on observe une violation du principe du maximum au début de la simulation avec Aster pour la loi de Mensi. Ceci peut s’expliquer (par analogie avec la thermique) par le «chochydrique» important dû à la façon d’imposer les conditions aux limites (concentration en eau imposée). Ce problème devrait pouvoir être résolu par l’utilisation des éléments lumpés de la même façon qu’en thermique. Pour la modélisation F, on peut conclure que la loi SECH_RFT est équivalente à la loi SECH_NAPPE avec la formule python fournie. **Fichier de commande Scilab** **Main.sci:** getf('/home/xxxx/librairie.sci') ; // PARAMETRES DE LA SIMULATION NUMERIQUE // // discretisation de la largeur x0 = 0.08 ; x = [:ref:`-0.080:0.001:+0.080 <-0.080:0.001:+0.080>`] ; [:ref:`n1 n2 `] = size(x) ; // teneur en eau initiale Cinit = 128.8 ; Ci = Cinit*ones(1,n2) ; // conditions aux limites à 50%HR CL = [:ref:`58.8 58.8 <58.8 58.8>`] ;. Ci(1) = CL(1) ; Ci($) = CL(2) ; Ci_bazant = Ci ; // pas de temps dt = 60 ; // [s] // coefficients de la loi de Bazant D1 = 3.0E-10 ; //[m2/s] a = 0.04 ; n = 6; TMAX = 5 ; // années //______________________________________________________________________ // // SIMULATION NUMERIQUE // j = 0 ; u=file('open','resultat_g','unknown'); for annee = 0:TMAX, annee for jour = 0:364, for heure = 0:23, minute = 0 ; for minute = 0:59, D_bazant = diffusion_bazant(D1,a,n,Cinit,58.8,Ci_bazant,293,293*ones(Ci),4700); Ci_bazant = linear_drying(D_bazant,Ci_bazant,CL,dt,x,"polar") ; if( (annee == 0 & jour == 0 & heure == 1 & minute == 0) | ... (annee == 0 & jour == 3 & heure == 0 & minute == 0) | ... (annee == 0 & jour == 28 & heure == 0 & minute == 0) | ... (annee == 1 & jour == 91 & heure == 0 & minute == 0) | ... (annee == 3 & jour == 0 & heure == 0 & minute == 0) | ... (annee == 5 & jour == 0 & heure == 0 & minute == 0) ) then, annee, jour, heure t=81:1:161; for tk=t, fprintf(u,'%6.3f %6.3f',x(tk),Ci_bazant(tk)) ; end, end, // if end, // for minute end, // for heure end, // for jour end, // for annee file('close',u) ; **Librairie.sci:** ________________________________________________________________ // // COEFFICIENT DE DIFFUSION NON LINEAIRE POUR LES SECHAGE DU BETON // LOI DE MENSI D(C) = a.exp(b.C) // ACTIVATION THERMIQUE D(C,T) = D(C,T0).(T/T0).exp[-Q/R*(1/T-1/T0)] // // a coefficient de la loi de Mensi // b coefficient de la loi de Mensi // C vecteur des teneurs en eau [-] // T0 température de référence [K] // T vecteur des températures [K] // Q_R Q/R (valant 4700 K) function D = diffusion_mensi(a,b,C,T0,T,Q_R), D = a*ones(C).*exp(b*C) ; D = D.*(T./(T0*ones(T))) ; D = D.*exp(Q_R*((ones(T)./T0)-(ones(T)./T))) ; endfunction, // // ________________________________________________________________ // // COEFFICIENT DE DIFFUSION NON LINEAIRE POUR LES SECHAGE DU BETON // LOI DE BAZANT // ACTIVATION THERMIQUE D(C,T) = D(C,T0).(T/T0).exp[-Q/R*(1/T-1/T0)] // // D1 coefficient de la loi de Bazant // a coefficient de la loi de Bazant(alpha) // n coefficient de la loi de Bazant // C0 teneur en eau à 100%HR // Cext teneur en eau du milieu environnant // C vecteur des teneurs en eau [-] // T0 température de référence [K] // T vecteur des températures [K] // Q_R Q/R (valant 4700 K-1) function D = diffusion_bazant(D1,a,n,C0,Cext,C,T0,T,Q_R), h = ones(C)-0.5*((C-C0*ones(C))/(Cext-C0))**2 ; D = (((1-a)*ones(C)./(ones(C)+(4**n)*(ones(C)-h)**n))+a*ones(C))*D1 ; D = D.*(T./(T0*ones(T))) ; D = D.*exp(Q_R*((ones(T)./T0)-(ones(T)./T))) ; endfunction, // ________________________________________________________________ // // DIFFUSION // Résolution par la méthode des différences finies // // D vecteur des coefficients de diffusion // Ci vecteur des teneurs en eau à l'instant j [-] // CL condition aux limites en xmin et xmax de type Dirichlet (C=C0) // dt pas de temps [s] // x vecteur des abscisses [m] // mode_ polaire/cartesien function Cf = linear_drying(D,Ci,CL,dt,x,mode_), [:ref:`n1,n2 `] = size(Ci) ; dx_ = zeros(1,n2-2) ; dx_(1:$) = (x(3:$)-x(1:$-2))*0.5 ; // Cf_ = (D*dt*(ones(dx_)./(dx_**2)).*(Ci(3:$)-2*Ci(2:$-1)+Ci(1:$-2)))+Ci(2:$-1); dx3 = ( ( ... (x(2:$-1)-x(1:$-2)).* ... (x(3:$ )-x(1:$-2)) ... ).* ... (x(3:$)-x(2:$-1)) ... ) ; d2C_dx2 = 2*( Ci(3:$ ).*( x(2:$-1)-x(1:$-2) ) ... -Ci(2:$-1).*( x(3:$ )-x(1:$-2) ) ... +Ci(1:$-2).*( x(3:$ )-x(2:$-1) ) ) ; d2C_dx2 = d2C_dx2./dx3 ; if (mode_ == "polar") then, dC_dx = ( Ci(3:$ ).*( x(2:$-1)-x(1:$-2) )**2 ... -Ci(1:$-2).*( x(3:$ )-x(2:$-1) )**2 ) ; // -Ci(2:$-1).*((x(2:$-1)-x(1:$-2) )**2 - (x(3:$ )-x(2:$-1))**2 )... dD_dx = ( D(3:$ ).*( x(2:$-1)-x(1:$-2) )**2 ... - D(1:$-2).*( x(3:$ )-x(2:$-1) )**2 ) ; // - D(2:$-1).*((x(2:$-1)-x(1:$-2) )**2 - (x(3:$ )-x(2:$-1))**2 )... dC_dx = dC_dx./dx3 ; dD_dx = dD_dx./dx3 ; i = find(x==0) ; [:ref:`k1 k2 `] = size(i) ; if(~(k1==0)) then, x(i) = x(i+1)/10, end, // printf("1er ordre %s ; 2ème ordre %s",string(min(dC_dx)),string(min(d2C_dx2))); d2C_dx2 = d2C_dx2 + dC_dx./x(2:$-1) ; end, Cf_ = Ci(2:$-1)+dt*(D(2:$-1).*d2C_dx2) ; if(mode_ == "polar") then, Cf_ = Cf_ +dt*(dD_dx.*dC_dx) ; end, Cf = zeros(1,n2) ; Cf(2:$-1) = Cf_ ; Cf(1) = CL(1) ; Cf($) = CL(2) ; endfunction, // // ________________________________________________________________ // **Comparaison Aster/Scilab** A2.1 SECH_MENSI/ SECH_GRANGER/SECH_NAPPE .. image:: images/10012C0E00003F9A00002AF98ADB245BC26FFA0A.svg :width: 521 :height: 351 .. _RefImage_10012C0E00003F9A00002AF98ADB245BC26FFA0A.svg: .. image:: images/1002CFB200003F9A00002AF9AC5F84FD3F0A0409.svg :width: 521 :height: 351 .. _RefImage_1002CFB200003F9A00002AF9AC5F84FD3F0A0409.svg: A2.2 SECH_BAZANT .. image:: images/1000A65800003FBA00002B18540144C59F410010.svg :width: 521 :height: 351 .. _RefImage_1000A65800003FBA00002B18540144C59F410010.svg: .. image:: images/10014EF800003FBA00002B1814619BCC3E120D1A.svg :width: 521 :height: 351 .. _RefImage_10014EF800003FBA00002B1814619BCC3E120D1A.svg: