code 1705

constants_lx_rk1.py

+# -*- coding: utf-8 -*-
+#________________________________Constant values used by Ling Xu_________________________________
+
+#Algorithm constants
+delta_t = 3e-2 #s. It is the time step of the algorithm.
+tf = 500 #s. It is the simulation time.
+
+#Specific constants of the problem
+Tfc = 343 #K. It is the temperature of the fuel cell.
+Pfc = 101325 #Pa. It is the pressure of the fuel cell.
+k = 3 #. It is the ratio of the inlet molar flow of O2 over the reaction molar flow of O2.
+Re = 1e-5 #ohm.m². It is the electron conduction resistance of the circuit.
+
+#The physical characteristics of the PEMFC
+    #Gas channel
+Lch = 1.298 #m. It is the length of the channel.
+Hgc= 1e-3 #m. It is the thickness of the gas channel.
+    #Gas diffusion layer
+Hgdl= 2.1e-4 #m. It is the thickness of the gas diffusion layer.
+epsilon_gdl  = 0.6 #It is the anode/cathode GDL porosity.
+    #Catalyst layer
+Hcl = 1e-5 #m. It is the thickness of the anode or cathode catalyst layer.
+epsilon_cl = 0.6 #It is the porosity of the catalyst layer, without units.
+i0_c = 0.1 #A.m-2.It is the exchange current density at the cathode.
+s_lim = 0.3 #It is the limit liquid saturation. 
+    #Membrane
+Hmem= 2.5e-5 #m. It is the thickness of the membrane.
+rho_mem = 1980 #kg.m-3. It is the density of the dry membrane.
+M_eq = 1.1e0 #kg.mol-1. It is the equivalent molar mass of ionomer.
+
+#The constants based on the interaction between the structure and the flow
+    #Vapor
+gamma_v = 1.3 #s-1. It is the sorption/desorption rate of vapor at the ionomer/CL interface.
+Dvc = 2.236e-5 #m².s-1. It is the vapor diffusivity in cathode.
+Dva = 5.457e-5 #m².s-1. It is the vapor diffusivity in anode.
+h_mv = 0.08 #m.s-1. It is the convective mass transfer coefficient of vapor.
+    #Liquid water
+gamma_l = 100 #s-1. It is the sorption/desorption rate of liquid at the ionomer/CL interface.
+Kc = 6.875e-13 #m². It is the anode/cathode GDL permeability for liquid water.
+theta_c = 110 #°. It is the contact angle of GDL for liquid water. 
+    #Dioxygen
+D_O2 = 2.9e-5 #m².s-1. It is the O2 diffusivity in cathode.
+h_O2 = 0.078 #m.s-1. It is the convective mass transfer coefficient of O2. 
+
+#The physical constants
+F = 96485 #C.mol-1. It is the Faraday constant.
+R = 8.314 #J. mol-1.K-1. It is the universal gas constant.
+M_H2O = 0.018 #kg.mol-1. It is the water molecular mass.
+rho_H2O = 1000 #kg.m-3. It is the water density. 
+nu_l = 3.7e-7 #m².s-1. It is the liquid water kinematic viscosity. 
+mu_l = 3.56e-4 #Pa.s. It is the liquid water dynamic viscosity. 
+sigma = 0.0625 #N.m-1. It is the water surface tension. 
+Csat = 10.9 #mol.m-3. It is the saturated vapor concentration at 70°C for a perfect gas.
+Psat = 31200 #Pa. It is the saturated partial pressure of vapor at 70°C.
+C_O2ref = 40.88 #mol.m-3. It is the reference concentration of oxygen.
+y_O2 = 0.2095 #. It is the molar fraction of O2 in dry air.
+
+#Mathematical factors
+Kshape = 5 #It is a shape mathematical factor, without physical sense.
+
+#Not understood yet
+mu_cg = 1.881e-5 #Pa.s. It is the cathode gas dynamic viscosity. 
+alpha_c = 0.5 #It is the transfer coefficient of the cathode.
\ No newline at end of file

current_density.py

+# -*- coding: utf-8 -*-
+import numpy as np
+#_______________________________________Current density function______________________________________
+
+#The current density raises in ~5s from 0.0 to 0.4 A.m-2 at tf/2.
+
+def current_density(n,delta_t):
+    
+    i_fc = np.zeros(n)
+    c = round(n/2) #The current density value changes at half-time.
+    tlim=c*delta_t #time at which the current density raises.
+    tau_c=1 #In 3*tau_c, 95% of the final value is reached.
+
+    for j in range(c):
+        i_fc[j]=0.0 #A.m-2. Start at 0.0 A.m-2
+        
+    for j in range(c,n):
+        t=j*delta_t #current time.
+        i_fc[j]=0.0+0.4*(1-np.exp( -(t-tlim)/tau_c )) #A.m-2. 
+        
+    return i_fc
\ No newline at end of file

main_rk1.py

+# -*- coding: utf-8 -*-
+from IPython import get_ipython
+import sys
+import time
+import numpy as np
+import numpy.polynomial.polynomial as nppol
+import matplotlib.pyplot as plt
+
+from current_density import current_density
+
+#Clean the console and reset all the variables
+get_ipython().magic('clear')
+get_ipython().magic('reset -f')
+
+#Start time
+start_time = time.time()
+
+#Select the value of the constants needed for the simulation.
+#"lx" : data used by Ling Xu in his articles + other data to complete them.
+version = "lx"
+
+#_______________________________Initialisation and constants value______________________________
+
+if version == "lx":
+    from constants_lx_rk1 import Re,delta_t,tf,Tfc,Pfc,k,Hgc,Hgdl,epsilon_gdl,Hcl,epsilon_cl,i0_c,s_lim,\
+         Hmem,rho_mem,M_eq,gamma_v,Dvc,Dva,h_mv,Kc,theta_c,D_O2,h_O2,F,R,M_H2O,rho_H2O,nu_l,mu_l,\
+         sigma,Csat,Psat,C_O2ref,y_O2,mu_cg,alpha_c
+    from transitory_functions_lx import lambda_eq,a_w,Cw,gamma,D,Q
+
+#Initialisation. Is has to start in the vapor regime ! (for now ?)
+Cini = 10.6
+C_O2ini = 1.55
+s5 = 0
+
+Cagc,C1,C5,Ccgc,C_O2cgc = Cini,Cini,Cini,Cini,C_O2ini
+lambda2,lambda3,lambda4 = lambda_eq(a_w(Cini)),lambda_eq(a_w(Cini)),lambda_eq(a_w(Cini))
+
+#Total number of step    
+n=round(tf/delta_t)
+
+#Target variables
+i_fc_t=current_density(n,delta_t) #table with all the values of i_fc
+Ja_t,Jc_t,J_O2_t,Jmem_acl_t,Jmem_ccl_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+Cagc_t,C0_t,C1_t,lambda2_t,lambda3_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+lambda4_t,s5_t,C5_t,C6_t,Ccgc_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+C_O2cgc_t,C_O2inter_t,C_O2ccl_t,Vcell_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+regime_t,a_w_c_t,lambda_eq_t = [],np.zeros(n),np.zeros(n)
+
+#_______________________________________Loop implementation_______________________________________
+for j in range (0,n):
+    
+    #____________________________________Current density value____________________________________
+    i_fc=i_fc_t[j]
+    
+    #_______________________Concentration flows functions of i_fc_______________________
+    
+    if j==0 or i_fc != i_fc_t[j-1]: #otherwise it is useless to recalculate them
+        #Intermediates 
+        Ccgc_eq = Csat - Psat/(Pfc-Psat)*i_fc/(2*F)*( Hgdl/(Dvc*epsilon_gdl**1.5) + 1/h_mv )
+        Phi_eq = Ccgc_eq*R*Tfc/Psat
+        
+        #Concentration flows at the GC
+        Wagc_in = Psat/(Pfc-Psat)*i_fc/(2*F*Hgc)
+        W_O2_in = k*i_fc/(4*F*Hgc)
+        W_O2_out = (k-1)*i_fc/(4*F*Hgc)
+        Wcgc_in = Phi_eq*Psat/(Pfc-Phi_eq*Psat)*k*i_fc/(4*F*y_O2*Hgc)
+        Wcgc_out = i_fc/(2*F*Hgc)*(1+ Psat/(Pfc-Psat) + Phi_eq*Psat/(Pfc-Phi_eq*Psat)*k/(2*y_O2))
+        
+        #Oxygen transport flow in the CGDL:
+        J_O2 = -i_fc/(4*F)
+    
+    #__________________________________Calculation of the flows__________________________________
+    
+    #Water transport flow in the membrane
+    Jmem_acl = 2.5/22*i_fc/F*lambda2 - 2*rho_mem/M_eq*D(lambda2)*(lambda3-lambda2)/Hmem
+    Jmem_ccl = 2.5/22*i_fc/F*lambda4 - 2*rho_mem/M_eq*D(lambda4)*(lambda4-lambda3)/Hmem
+    
+    #Intermediate values
+    Cw_c = Cw(s5,C5)
+    a_w_a, a_w_c= a_w(C1), a_w(Cw_c)
+    gamma_c = gamma(a_w_c)
+    
+    #Water transport flow in the two GDLs:
+    Ja = gamma_v*epsilon_cl*Hcl*rho_mem/M_eq*(lambda_eq(a_w_a)-lambda2)
+    
+    if lambda4 >= 22:
+        Jc = max( Jmem_ccl + i_fc/(2*F) ,\
+                  gamma_c*epsilon_cl*Hcl*rho_mem/M_eq*(lambda4-lambda_eq(a_w_c)) )
+    else:
+        Jc = gamma_c*epsilon_cl*Hcl*rho_mem/M_eq*(lambda4-lambda_eq(a_w_c))
+
+    #________________________Determination of the water transport regime_______________________
+    
+    if C5 < Csat :
+        regime = "V"
+    elif C5 >= Csat and Ccgc < Csat :
+        regime = "M"
+    elif C5 >= Csat and Ccgc >= Csat :
+        regime = "L"
+    else:
+        print(C5,Ccgc,Csat)
+        sys.exit("The concentrations are not relevant")
+
+    #________________________Calculation of the instantaneous quantities________________________
+    
+    #At the anode :
+    C0 = Cagc - Ja/h_mv
+    if C0>Csat:
+        C0=Csat
+    if C0<0:
+        print (C0,j)
+        sys.exit("C0 concentration cannot be negative.")
+        
+    C1 = C0 - Ja*Hgdl/(Dva*epsilon_gdl**1.5)
+    if C1>Csat:
+        C1=Csat
+    if C1<0:
+        print (C1,j)
+        sys.exit("C1 concentration cannot be negative.")
+        
+    #At the cathode :
+    if regime == "V":
+        s5,s6 = 0,0
+        
+        C6 = Ccgc + Jc/h_mv
+        if C6>Csat:
+            C6=Csat
+        if C6<0:
+            print (C6,j)
+            sys.exit("C6 concentration cannot be negative.")
+            
+        C5 = C6 + Jc*Hgdl/(Dvc*epsilon_gdl**1.5)
+        if C5>Csat:
+            C5=Csat
+        if C5<0:
+            print (C5,j)
+            sys.exit("5 concentration cannot be negative.")
+        
+        C_O2inter = C_O2cgc + J_O2/h_O2
+        if C_O2inter<0:
+            print (C_O2inter,j)
+            sys.exit("C_O2inter concentration cannot be negative.")
+            
+        C_O2ccl = C_O2inter + J_O2*Hgdl/(D_O2*epsilon_gdl**1.5)
+        if C_O2ccl<0:
+            print (C_O2ccl,j)
+            sys.exit("C_O2ccl concentration cannot be negative.")
+        
+        
+    elif regime == "M":
+        C5,s6 = Csat,0
+        
+        C6 = Ccgc + Jc/h_mv
+        if C6>Csat:
+            C6=Csat
+        if C6<0:
+            print (C6,j)
+            sys.exit("C6 concentration cannot be negative.")
+
+        xs = Hgdl + Dvc*epsilon_gdl**1.5*(C6-Csat)/Jc
+        
+        a0 = M_H2O*Jc*xs / (sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5)
+        solv = nppol.polyroots([a0, 0, 0, 0, 0.35425, -0.848, 0.6315])
+        c=0; #number of time it enters in the following if :
+        for i in solv:
+            if abs(i.imag) < 1e-10 and i.real>=0.0 and i.real<=1.0001: 
+                if s5 >= 1.0:           
+                    s5=0.999999
+                else :
+                    s5 = i.real
+                c += 1
+        if c != 1 :
+            print (solv,j)
+            sys.exit("The polynomial solver gives more or less than a unique useful value.")
+            
+        C_O2inter = C_O2cgc + J_O2/h_O2
+        if C_O2inter<0:
+            print (C_O2inter,j)
+            sys.exit("C_O2inter concentration cannot be negative.")
+            
+        C_O2_xs = C_O2inter + J_O2*(Hgdl-xs)/(D_O2*epsilon_gdl**1.5)
+        if C_O2_xs<0:
+            print (C_O2_xs,j)
+            sys.exit("C_O2_xs concentration cannot be negative.")
+            
+        C_O2ccl = C_O2_xs - (J_O2*sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5) / \
+                            (M_H2O*Jc*D_O2*epsilon_gdl**1.5) * (Q(s5)-Q(0))
+        if C_O2ccl<0:
+            print (C_O2ccl,j)
+            sys.exit("C_O2ccl concentration cannot be negative.")
+                
+    elif regime == "L":
+        C5,C6 = Csat,Csat
+        
+        s6 = 1/(1+ (rho_H2O*Wcgc_out*Hgc*mu_cg / (M_H2O*abs(Jc)*Ccgc_eq*mu_l))**(1/3) )
+        
+        a0 = M_H2O*Jc*Hgdl / (sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5) \
+            -0.35425*s6**4 + 0.848*s6**5 -0.6315*s6**6
+        solv = nppol.polyroots([a0, 0, 0, 0, 0.35425, -0.848, 0.6315])
+        c=0; #number of time it enters in the following if :
+        for i in solv:
+            if abs(i.imag) < 1e-10 and i.real>=0.0 and i.real<=1.0001:
+                if s5 >= 1.0:                              
+                    s5=0.999999
+                else :
+                    s5 = i.real
+                c += 1
+        if c != 1 :
+            print (solv,j)
+            sys.exit("The polynomial solver gives more or less than a unique useful value.")
+            
+        C_O2inter = C_O2cgc + J_O2/h_O2
+        if C_O2inter<0:
+            print (C_O2inter,j)
+            sys.exit("C_O2inter concentration cannot be negative.")
+            
+        C_O2ccl = C_O2inter - (J_O2*sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5) / \
+                            (M_H2O*Jc*D_O2*epsilon_gdl**1.5) * (Q(s5)-Q(s6))
+        if C_O2ccl<0:
+            print (C_O2ccl,j)
+            sys.exit("C_O2ccl concentration cannot be negative.")
+    
+    #___________________________Control of the water transport regime___________________________
+    
+    #I am going to look at this later.
+    # if regime == "V":
+    #     if C5 > Csat :
+    #         regime == "M"
+    # elif regime == "M":
+    #     if C5 < Csat:
+    #         regime == "V"
+    #     if C5 >= Csat and Ccgc >= Csat :
+    #         regime == "L"
+    # elif regime == "L":    
+    #     if  Ccgc < Csat:
+    #         regime == "M"
+    
+    #_____________________________Calculation of the slow quantities_____________________________
+    
+    Cagc += delta_t*(Wagc_in - Ja/Hgc)
+    if Cagc>Csat:
+        Cagc=Csat
+    if Cagc<0:
+        print (Cagc,j)
+        sys.exit("Cagc concentration cannot be negative.")
+        
+    Ccgc += delta_t*(Wcgc_in + Jc/Hgc - Wcgc_out)
+    if Ccgc>Csat:
+        Ccgc=Csat
+    if Ccgc<0:
+        print (Ccgc,j)
+        sys.exit("Ccgc concentration cannot be negative.")
+        
+    C_O2cgc += delta_t*(W_O2_in + J_O2/Hgc - W_O2_out)
+    if C_O2cgc<0:
+            print (C_O2cgc,j)
+            sys.exit("C_O2cgc concentration cannot be negative.")
+    
+    lambda2 += delta_t*M_eq/(Hcl*rho_mem)*(Ja-Jmem_acl)
+    if lambda2<0:
+        print (lambda2,j)
+        sys.exit("lambda2 concentration cannot be negative.")
+        
+    lambda3 += delta_t*M_eq/(Hmem*rho_mem)*(Jmem_acl-Jmem_ccl)
+    if lambda3<0:
+        print (lambda3,j)
+        sys.exit("lambda3 concentration cannot be negative.")
+        
+    lambda4 += delta_t*M_eq/(Hcl*rho_mem)*(Jmem_ccl-Jc+i_fc/(2*F))
+    if lambda4<0:
+        print (lambda4,j)
+        sys.exit("lambda4 concentration cannot be negative.")
+    
+    #______________________________Calculation of the cell voltage______________________________
+    
+    Ueq = 1.229 - 8.5e-4*(Tfc-298.15)
+    
+    if i_fc == 0.0:
+        eta_c = 0.0
+        Rccl = 0.0
+        Rmem = 0.0
+    else:
+        eta_c = -R*Tfc/(alpha_c*F)*np.log(s_lim/(s_lim-s5) * C_O2ref/C_O2ccl * i_fc/i0_c )
+        Rccl = Hcl/( epsilon_cl**1.5 * (0.5139*lambda4-0.326) * np.exp(1268*(1/303.15 - 1/Tfc)) )
+        Rmem = Hmem/( 2*np.exp(1268*(1/303.15 - 1/Tfc)) ) *\
+        ( np.log(0.5139*lambda4-0.326)-np.log(0.5139*lambda3-0.326) ) / (0.5139*(lambda4-lambda3))
+        
+    Rp = Rmem + Rccl/3
+    
+    Vcell = Ueq + eta_c - i_fc*(Rp+Re)
+    
+    #_______________________________Record of the target variables_______________________________
+    
+    Ja_t[j],Jc_t[j],J_O2_t[j],Jmem_acl_t[j],Jmem_ccl_t[j] = Ja,Jc,J_O2,Jmem_acl,Jmem_ccl
+    Cagc_t[j],C0_t[j],C1_t[j],lambda2_t[j],lambda3_t[j] = Cagc,C0,C1,lambda2,lambda3
+    lambda4_t[j],s5_t[j],C5_t[j],C6_t[j],Ccgc_t[j] = lambda4,s5,C5,C6,Ccgc
+    C_O2cgc_t[j],C_O2inter_t[j],C_O2ccl_t[j],Vcell_t[j] = C_O2cgc,C_O2inter,C_O2ccl,Vcell
+    a_w_c_t[j],lambda_eq_t[j] = a_w_c,lambda_eq(a_w_c)
+    regime_t += [regime]
+
+#_____________________________________________Display_____________________________________________
+
+t = np.linspace(0,tf,n)
+
+#Zoom on the timeline [0.49*tf,0.52*tf]s
+n1_z1=round(0.49*tf/delta_t) #number of step to reach 245s.
+n2_z1=round(0.52*tf/delta_t) #number of step to reach 260s.
+delta_n = n2_z1-n1_z1 #number of step from 245s to reach 260s.
+t_z1 = np.linspace(0.49*tf,0.52*tf,delta_n)
+
+#Zoom on the timeline [0.4*tf,0.7*tf]s
+n1_z2=round(0.4*tf/delta_t) #number of step to reach 245s.
+n2_z2=round(0.7*tf/delta_t) #number of step to reach 260s.
+delta_n = n2_z2-n1_z2 #number of step from 245s to reach 260s.
+t_z2 = np.linspace(0.4*tf,0.7*tf,delta_n)
+
+#Zoom on the timeline [0.0*tf,0.01*tf]s
+n1_z3=round(0.0*tf/delta_t) #number of step to reach 245s.
+n2_z3=round(0.01*tf/delta_t) #number of step to reach 260s.
+delta_n = n2_z3-n1_z3 #number of step from 245s to reach 260s.
+t_z3 = np.linspace(0.0*tf,0.01*tf,delta_n)
+
+# #Test
+# n1_z4=n//2
+# n2_z4=25038
+# delta_n = n2_z4-n1_z4 #number of step from 245s to reach 260s.
+# t_z4 = np.linspace(n1_z4*delta_t,n2_z4*delta_t,delta_n)
+
+# plt.figure(80)
+# plt.plot(t_z4,a_w_c_t[n1_z4:n2_z4],label='a_w_c')
+# plt.legend()
+
+# plt.figure(81)
+# plt.plot(t_z4,lambda_eq_t[n1_z4:n2_z4],label='lambda_eq')
+# plt.legend()
+
+# #Figure 1 : zoom1 on i_fc
+# plt.figure(1)
+# i_fc_z=i_fc_t[n1_z1:n2_z1] #It is a zoom on i_fc.
+# plt.plot(t_z1,i_fc_z,label='i_fc')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Current density (A/m²)')
+# plt.title('The current density\nbehaviour over time')
+# plt.legend()
+
+#Figure 2 : the different flows J
+plt.figure(2)
+i_fc_flow = i_fc_t/(2*F)
+plt.plot(t_z1,Ja_t[n1_z1:n2_z1],label='Ja')
+plt.plot(t_z1,Jc_t[n1_z1:n2_z1],label='Jc')
+plt.plot(t_z1,i_fc_flow[n1_z1:n2_z1],label='i_fc/2F')
+plt.plot(t_z1,Jmem_acl_t[n1_z1:n2_z1],label='Jmem_acl')
+plt.plot(t_z1,Jmem_ccl_t[n1_z1:n2_z1],label='Jmem_ccl')
+# plt.plot(t,J_O2_t,label='J_O2')
+plt.xlabel('Time (s)')
+plt.ylabel('Flow (mol/m²/s)')
+plt.title('The flow\nbehaviour over time')
+plt.legend()
+
+# #Figure 3 : the different vapor concentration C
+# plt.figure(3)
+# plt.plot(t,Cagc_t,label='Cagc')
+# plt.plot(t,C0_t,label='C0')
+# plt.plot(t,C1_t,label='C1')
+# plt.plot(t,C5_t,label='C5')
+# plt.plot(t,C6_t,label='C6')
+# plt.plot(t,Ccgc_t,label='Ccgc')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Vapor concentration (mol/m^3)')
+# plt.title('The vapor concentration\nbehaviour over time')
+# plt.legend()
+
+# #Figure 4 : the water content
+# plt.figure(4)
+# plt.plot(t,lambda2_t,label='lambda2')
+# plt.plot(t,lambda3_t,label='lambda3')
+# plt.plot(t,lambda4_t,label='lambda4')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Water content')
+# plt.title('The three water content\nbehaviour over time')
+# plt.legend()
+
+# # Figure 5 : the liquid water saturation s5
+# plt.figure(5)
+# plt.plot(t,s5_t,label='s5')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Liquid water saturation at CCL')
+# plt.title('The CCL liquid water saturation\nbehaviour over time')
+# plt.legend()
+
+# #Figure 6 : the O2 concentrations C_O2
+# plt.figure(6)
+# plt.plot(t,C_O2cgc_t,label='C_O2cgc')
+# plt.plot(t,C_O2inter_t,label='C_O2inter')
+# plt.plot(t,C_O2ccl_t,label='C_O2ccl')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Oxygen concentration (mol/m3)')
+# plt.title('The oxygen concentrations\nbehaviour over time')
+# plt.legend()
+
+# # Figure 7 : zoom1 on Vcell
+# plt.figure(7)
+# Vcell_z=Vcell_t[n1_z1:n2_z1] #It is a zoom on Vcell.
+# plt.plot(t_z1,Vcell_z,label='Vcell')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Cell voltage (V)')
+# plt.title('The cell voltage\nbehaviour over time')
+# plt.legend()
+
+# #Figure 8 : the different vapor concentration C zoom1
+# plt.figure(8)
+# plt.plot(t_z1,Cagc_t[n1_z1:n2_z1],label='Cagc')
+# plt.plot(t_z1,C0_t[n1_z1:n2_z1],label='C0')
+# plt.plot(t_z1,C1_t[n1_z1:n2_z1],label='C1')
+# plt.plot(t_z1,C5_t[n1_z1:n2_z1],label='C5')
+# plt.plot(t_z1,C6_t[n1_z1:n2_z1],label='C6')
+# plt.plot(t_z1,Ccgc_t[n1_z1:n2_z1],label='Ccgc')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Vapor concentration (mol/m^3)')
+# plt.title('The vapor concentration\nbehaviour over time')
+# plt.legend()
+
+# #Figure 9 : the water content zoom 1
+# plt.figure(9)
+# plt.plot(t_z1,lambda2_t[n1_z1:n2_z1],label='lambda2')
+# plt.plot(t_z1,lambda3_t[n1_z1:n2_z1],label='lambda3')
+# plt.plot(t_z1,lambda4_t[n1_z1:n2_z1],label='lambda4')
+# # plt.plot(t_z1,lambda_eq(a_w_c_t)[n1_z1:n2_z1],label='lambda_eq(a_w_c)')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Water content')
+# plt.title('The three water content\nbehaviour over time')
+# plt.legend()
+
+# #Figure 10 : test
+# plt.figure(10)
+# test=C5_t
+# plt.plot(t_z1,test[n1_z1:n2_z1],label='C5')
+# plt.xlabel('Time (s)')
+# plt.legend()
+
+#Algorithm time
+algo_time = time.time() - start_time
+print('Time of the algorithm in second :',algo_time)
\ No newline at end of file

transitory_functions_lx.py

+# -*- coding: utf-8 -*-
+import numpy as np
+#_________________________________Initialisation and constants value________________________________
+from constants_lx_rk1 import Tfc,gamma_v,gamma_l,M_H2O,rho_H2O,Csat,Kshape
+
+#_______________________________________Transitory functions_______________________________________
+
+def lambda_eq(a_w):
+    return 0.5*( 0.043 +17.81*a_w -39.85*a_w**2 +36.0*a_w**3 )*( 1 -np.tanh(100*(a_w -1)) ) \
+         + 0.5*( 14 +8*( 1-np.exp(-Kshape*(a_w -1)) ) )*( 1+np.tanh(100*(a_w -1)) )
+        
+def a_w(Cw):
+    return Cw/Csat
+
+def Cw(s,C):
+    return C+s*rho_H2O/M_H2O
+
+def gamma(a_w):
+    ome = omega(a_w)
+    if a_w <=1 :
+        return gamma_v
+    else:
+        return (1-ome)*gamma_v + ome*gamma_l
+
+def omega(a_w):
+    lambdaEq = lambda_eq(a_w)
+    return (lambdaEq-14)/(22-14)*0.5 + (a_w -1)/(3 -1)*0.5
+
+def D(lambdaa): #the diffusion coefficient of water in the membrane
+    if lambdaa <= 3:
+        return 3.1e-7*lambdaa*(np.exp(0.28*lambdaa)-1)*np.exp(-2346/Tfc)
+    else:
+        return 4.17e-8*lambdaa*(161*np.exp(-lambdaa)+1)*np.exp(-2346/Tfc)
+    
+def Q(s):
+    return -(14735*s**5-300*s**4+9439*s**3+18878*s**2+75512*s-151024) / (17500*np.sqrt(1-s))
\ No newline at end of file