function [out]=modele_emballage(); %% ***************************************** % GCH 2535 - Devoir No2 % Equipe: Minh Dung Nguyen % % % %% ***************************************** close all clear all clc % close all global R; %% ***************************************** % Prétraitement % ***************************************** %Données d'entrée du problème %Constantes physiques Meau=(1.01*2+16.00)*10^-3; %Masse moléculaire de l'eau (kg/mol) R=8.314; %Constante des gaz parfaits (m^3.Pa/K/mol) isotherm='GAB'; %Type d'isotherme de sorption (GAB, BET, lin) if lower(isotherm)=='gab' %isotherme GAB et ces constantes kGAB=0.897; %Constante pour l'isotherme GAB isotherm CGAB=2.31; %Constante de Guggenheim pour l'isotherme GAB M0GAB=0.066; %Moisture content of monolayer for GAB isotherm (kg water/kg solids) end %Conditions initiales RHi=90; %Taux d'humidité ambient(%) Ti=273.15+38; %Température (K) p0i=exp(23.4795-3990.56/(Ti-273.15+233.833)); %Pression de vapeur saturée de l'eau pure à T (Pa) pai=RHi/100*p0i; %Pression de vapeur ambiante (Pa) Mi=0.0357; %Contenu initial en humidité des croustilles (kg eau/kg solides) if lower(isotherm)=='gab' %Calcul de la pression de vapeur initial à l'équilibre avec les croustilles à l'intérieur de l'emballage (Pa) pfi=p0i*(-(kGAB*(CGAB-CGAB*M0GAB/Mi-2))-sqrt((kGAB*(CGAB-CGAB*M0GAB/Mi-2))^2-4*kGAB^2*(1-CGAB)))/(2*kGAB^2*(1-CGAB)); end %Données sur l'emballage Rint= .0326 ; %Rayon intérieur de l'emballage (m) A= 0.2025*2*pi*Rint; %Surface intérieure de l'emballage (m^2) ms=0.160/(1+Mi); %Masse de solide sec dans les croustilles (kg) %Couche 1 de PE X(1)=45e-6; %Épaisseur (m) [D(1),S(1),P(1)]=DSPFun(1,Ti); %Calcul de la diffusion (m^2/s), de la solubilité (mol/(m^3·Pa)) et de la perméabilité (mol·m/(m^2·s·Pa)) %Couche 2 d'Al X(2)=5e-6; %Épaisseur (m) [D(2),S(2),P(2)]=DSPFun(2,Ti); %Calcul de la diffusion (m^2/s), de la solubilité (mol/(m^3·Pa)) et de la perméabilité (mol·m/(m^2·s·Pa)) %Couche 3 de PE X(3)=40e-6; %Épaisseur (m) [D(3),S(3),P(3)]=DSPFun(3,Ti); %Calcul de la diffusion (m^2/s), de la solubilité (mol/(m^3·Pa)) et de la perméabilité (mol·m/(m^2·s·Pa)) %Couche 4 de carton X(4)=300e-6; %Épaisseur (m) [D(4),S(4),P(4)]=DSPFun(4,Ti); %Calcul de la diffusion (m^2/s), de la solubilité (mol/(m^3·Pa)) et de la perméabilité (mol·m/(m^2·s·Pa)) %Couche 5 de couchage d'impression X(5)=10e-6; %Épaisseur (m) [D(5),S(5),P(5)]=DSPFun(5,Ti); %Calcul de la diffusion (m^2/s), de la solubilité (mol/(m^3·Pa)) et de la perméabilité (mol·m/(m^2·s·Pa)) XT=sum(X); %Total thickness of packaging (m) %Discrétisation du domaine J=7; %nombre de noeuds à l'intérieur de chaque couche excluant les noeuds aux interfaces et aux frontières(-) dr=X./(J+1); %pas en espace pour chacune des couches (m) % position des differents noeuds par rapport au centre de l'emballage (m) nbcouches=5; nn=nbcouches*(J+1)+1; r=repmat(dr,[J+1 1]).*repmat([1:J+1]',[1 nbcouches]); r=r+repmat([0 r(end,1:end-1)*triu(ones(nbcouches-1))]+Rint,[J+1 1]);r=[Rint; r(:)]'; % Initilisation des vecteurs/matrices p=(ones(1,nbcouches*(J+1)+1)*pai)'; p(1)=pfi; %% ***************************************** % Solveur % ***************************************** % Boucle en temps stemps=365; %Temps de simulation (jours) tfin=stemps*24*3600; % Temps final (s) dt=1*24*3600; % pas de temps (s) %======= construction de la matrice ================= drn=repmat(dr,[J+1 1]); drn=[drn(:); 1]'; Dn=repmat(D,[J+1 1]); Dn=[Dn(:); 1]'; % Noeuds entre les interfaces Noeuds aux interfaces cst1=(dt.*Dn./(2.*r.*drn)-dt.*Dn./(drn.^2)); cst1(J+2:J+1:end-1)=-D(1:end-1)./dr(1:end-1); cst2=2.*Dn.*dt./drn.^2+1; cst2(J+2:J+1:end-1)= D(1:end-1)./dr(1:end-1)+D(2:end)./dr(2:end); cst3=(-dt.*Dn./(2.*r.*drn)-dt.*Dn./(drn.^2)); cst3(J+2:J+1:end-1)=-D(2:end)./dr(2:end); matA=full(gallery('tridiag',nn,cst1(2:end),cst2,cst3(1:end-1))); matA(1,:)=0; matA(end,:)=0; matA(1,1)=1; matA(end,end)=1; %======= construction matrice fin ==================== %vecteur de Masse M=zeros(round(tfin/dt),1); M(1)=Mi; %======= résolution ================================== for t=1:1:round(tfin/dt) %résolution prression vapeur vectB=p; pf=p0i*(-(kGAB*(CGAB-CGAB*M0GAB/M(t)-2))-sqrt((kGAB*(CGAB-CGAB*M0GAB/M(t)-2))^2-4*kGAB^2*(1-CGAB)))/(2*kGAB^2*(1-CGAB)); vectB(1)=pf; vectB(2+J:J+1:end-1)=0; p=linsolve(matA,vectB); % résolution masse eau croustille M(t+1)=dt*Meau/ms*(P(1)*A*(-p(3)+4*p(2)-3*p(1))/(2*dr(1)))+M(t); if (M(t+1)>0.10&&M(t)<=0.10), temps=(t*dt)/3600/24+(0.10-M(t))/(M(t+1)-M(t)); fprintf('Temps pour atteindre la limite de 0.10 kg eau/kg solide: %8.2f jours\n',temps) end end fprintf('Contenu en humidité final: %9.4f kg eau/kg solide\n',M(end)) fprintf('Poids final des croustilles: %9.4f g\n',ms*(1+M(end))*1000) %======= résolution fin ============================= %% ***************************************** % Post-traitement des résultats % ***************************************** % % NB: Vous pouvez mettre ces graphiques à votre main et modifier cette % partie comme il vous plaira % % Graphique du contenu en humidité dans les croustilles en fonction du temps % hFig = figure(2); set(hFig, 'Position', [200 50 1100 750]) subplot(2,1,1); hold on plot([0:1:round(tfin/dt)+1]*dt/3600/24,[Mi; M],'b-','LineWidth',2); title('Contenu d''humidité dans les croustilles','fontsize',18); xlabel('Temps (jours)','fontsize',18); ylabel('M (kg eau/kg solide)','fontsize',18); set(gca,'fontsize',16,'LineWidth',2) box on %Graphique de la pression de vapeur d''eau dans l'emballage subplot(2,1,2); plot((r-Rint)*1e6,p,'b-','LineWidth',2); title(['Pression de vapeur d''eau dans l''emballage après ',num2str(stemps),' jours'],'fontsize',18); xlabel('Épaisseur de l''emballage (µm)','fontsize',18); ylabel('p (Pa)','fontsize',18); set(gca,'fontsize',16,'LineWidth',2) box on end %%Diffusivity, Solubility and permeability of Water Vapour in Packaging - Function File function [D,S,P]=DSPFun(n,T) % DESCRIPTION: cette routine calcule les valeurs de permeabilité, % solubilité et diffusivité à une température T donnée pour une couche donnée n à % partir de valeurs de ces mêmes propriétés qui y sont données à température de référence de 38oC. % % n : numéro de la couche % T : température (K) % D : diffusivité (m^2/s) % S : solubilité (mol/(m^3·Pa)) % P : perméabilité (mol·m/(m^2·s·Pa) % % ATTENTION: si vous désirez changer les propriétés de vos matériaux, je % vous suggère de seulement changer les valeurs de référence à 38oC pour la % perméabilité (c-a-d Pref(n)) seulement % global R; switch n case 1 Pref(n)=9.4e-14; % Reference permeability of water vapour in layer 1 (PE) of packaging (mol·m/(m^2·s·Pa)) Sref(n)=0.001; %Reference solubility of water vapour in layer 1 (PE) of packaging (mol/m^3/Pa) Dref(n)=Pref(n)/Sref(n); %Reference diffusivity of water vapour in layer 1 of packaging (m^2/s) TrefS(n)=273.15+38; %Temperature of reference solubility of water vapour in layer 1 of packaging (K) dHS(n)=-3.5e4; %Partial molar enthalpy of sorption of water vapour in layer 1 of packaging (J/mol) TrefD(n)=273.15+38; %Temperature of reference diffusivity of water vapour in layer 1 of packaging (K) ED(n)=2.5e4-dHS(n); %Activation energy of diffusion of water vapour in layer 1 of packaging (J/mol) case 2 Pref(n)=2e-15; % Reference permeability of water vapour in layer 2 (Al) of packaging (mol·m/(m^2·s·Pa)) Sref(n)=0.001; %Reference solubility of water vapour in layer 2 (Al) of packaging (mol/m^3/Pa) Dref(n)=Pref(n)/Sref(n); %Reference diffusivity of water vapour in layer 2 of packaging (m^2/s) TrefS(n)=273.15+38; %Temperature of reference solubility of water vapour in layer 2 of packaging (K) dHS(n)=-3.5e4; %Partial molar enthalpy of sorption of water vapour in layer 2 of packaging (J/mol) TrefD(n)=273.15+38; %Temperature of reference diffusivity of water vapour in layer 2 of packaging (K) ED(n)=2.5e4-dHS(n); %Activation energy of diffusion of water vapour in layer 2 of packaging (J/mol) case 3 Pref(n)=9.4e-14; % Reference permeability of water vapour in layer 3 (PE) of packaging (mol·m/(m^2·s·Pa)) Sref(n)=0.001; %Reference solubility of water vapour in layer 3 (PE) of packaging (mol/m^3/Pa) Dref(n)=Pref(n)/Sref(n); %Reference diffusivity of water vapour in layer 3 of packaging (m^2/s) TrefS(n)=273.15+38; %Temperature of reference solubility of water vapour in layer 3 of packaging (K) dHS(n)=-3.5e4; %Partial molar enthalpy of sorption of water vapour in layer 3 of packaging (J/mol) TrefD(n)=273.15+38; %Temperature of reference diffusivity of water vapour in layer 3 of packaging (K) ED(n)=2.5e4-dHS(n); %Activation energy of diffusion of water vapour in layer 3 of packaging (J/mol) case 4 Pref(n)=1.6e-11; % Reference permeability of water vapour in layer 4 (Carton) of packaging (mol·m/(m^2·s·Pa)) Sref(n)=0.001; %Reference solubility of water vapour in layer 4 (Carton) of packaging (mol/m^3/Pa) Dref(n)=Pref(n)/Sref(n); %Reference diffusivity of water vapour in layer 4 of packaging (m^2/s) TrefS(n)=273.15+38; %Temperature of reference solubility of water vapour in layer 4 of packaging (K) dHS(n)=-3.5e4; %Partial molar enthalpy of sorption of water vapour in layer 4 of packaging (J/mol) TrefD(n)=273.15+38; %Temperature of reference diffusivity of water vapour in layer 4 of packaging (K) ED(n)= 2.5e4-dHS(n); %Activation energy of diffusion of water vapour in layer 4 of packaging (J/mol) case 5 Pref(n)=8.5e-13; % Reference permeability of water vapour in layer 5 (Couchage) of packaging (mol·m/(m^2·s·Pa)) Sref(n)=0.001; %Reference solubility of water vapour in layer 5 (Couchage) of packaging (mol/m^3/Pa) Dref(n)=Pref(n)/Sref(n); %Reference diffusivity of water vapour in layer 5 of packaging (m^2/s) TrefS(n)=273.15+38; %Temperature of reference solubility of water vapour in layer 5 of packaging (K) dHS(n)=-3.5e4; %Partial molar enthalpy of sorption of water vapour in layer 5 of packaging (J/mol) TrefD(n)=273.15+38; %Temperature of reference diffusivity of water vapour in layer 5 of packaging (K) ED(n)= 2.5e4-dHS(n); %Activation energy of diffusion of water vapour in layer 5 of packaging (J/mol) end D=Dref(n)*exp(ED(n)/R*(1/TrefD(n)-1/T)); S=Sref(n)*exp(dHS(n)/R*(1/TrefS(n)-1/T)); P=D*S; end