Prev: date on x axis feather plot
Next: Need Help Attempted to access hx(0,0); index must be a positive integer or logical.
From: Torsten Hennig on 12 Jun 2010 04:45
Hello Alina,

two more things I would like to mention:

1.
I think the starting values for the concentrations at
the liquid-air interface are not consistent with the
algebraic equations you define in function dcdt:

In my opinion, they should read:
c0(N+1,1) = (2*Dl*c0(N,1)/h(N)+2*Dg*c0(N+3,1)/h(N+3))/
(2*Dl/h(N)+2*Dg*Kgl/h(N+3))
c0(N+2,1) = Kgl*c0(N+1,1)

2.
Because of the storage requirements,
I would not set up the Jacobian matrix in function dcdt.
Instead, I would provide J*c directly.
E.g. instead of
J(1,1) = -(2*Dl/h(1))/(h(1)+h(2))
J(1,2) = (2*Dl/h(1))/(h(1)+h(2)),
I would prescribe
dcdt(1) = -(2*Dl/h(1))/(h(1)+h(2))*c(1) +
(2*Dl/h(1))/(h(1)+h(2))*c(2)
and so on.

Best wishes
Torsten.


> Hello Alina,
> I found one error in your discretization (see below)
> After the change, all went well with my solver.
> Hope MATLAB's ODE solver now also gives satisfactory
> results.
>
> Best wishes
> Torsten.
>
>
> > Hello ,
> >
> > I am trying to simulate the transfer of a compound
> > between two compartments separated by a membrane!
> I'm
> > using constant temperature and pressure. In the
> > inferior compartment I have 2 phases, an aqueous
> > solution of my compound and the air above it
> (there
> > is no agitation in this compartment, that is why I
> > consider the molecular diffusion as transport
> > mechanism in the liquid and the gas) and in the
> > superior compartment I have an air flow such that
> the
> > volume of air is renewed each minute (the air flow
> is
> > parallel to the membrane and the transport
> mechanism
> > is the forced convection), that is why I consider
> the
> > concentration of my compound 0 in the superior
> > compartment. So my compound is vaporizing into the
> > air above the liquid, it's passing through the
> > membrane and is swept by the air flow to the
> > atmosphere. The equations for the transport in the
> > liquid, gas and membrane have the form of
> > the Fick's second law (the concentration is a
> > function of time and space), with constant
> diffusion
> > coefficients for each phase!
> > My equations are :
> > - Liquid phase - dcl/dt = Dl * d2cl/dx2
> > -initial condition : cl = c0
> > - boundary condition :
> > - left side - x = 0 : dcl/dx = 0 ;
> > - right side- x = hliq :
> > Dl * dcl/dx =Dg* dcg/dx
> > cl * Kgl = cg , Kgl
> > cl * Kgl = cg , Kgl =partition coefficient
> >
> > - Gas phase - dcg/dt = Dg *
> d2cg/dx2
> > 2
> > - initial condition : cg = 0
> > - boundary condition :
> > - left side - x = hliq :
> > Dl * dcl/dx =Dg* dcg/dx
> > cl * Kgl = cg , Kgl
> > cl * Kgl = cg , Kgl =partition coefficient
> >
> > right side - x = hcomp1 :
> > Dg * dcg/dx =Dm* dcm/dx
> > cg = cm
> > - Membrane - dcm/dt = Dm * d2cm/dx2
> > - initial condition : cm = 0
> > - boundary condition :
> > - left side - x = hcomp1 :
> > S*Dg * dcg/dx =Smem*Dm*
> > S*Dg * dcg/dx =Smem*Dm* dcm/dx
> > cg = cm
> > - right side - x = ht :
> > Smem*Dm * dcm/dx = S*Kext
> > Smem*Dm * dcm/dx = S*Kext *(cmext – cair)
> > Cm = cmext ; cair = 0 ; Kext
> > = cmext ; cair = 0 ; Kext – mass transfer
> > coefficient superior compartment side
> >
> > The finite volumes method – for the equations:
> > -liquid
> > VL,1*dcL,1/dt = 0 –
> > 1; S*Dl*(cL,1-cL,2)/((hL,1+hL,2)/2)
> > VL,i*dcL,i/dt =
> > = S*Dl*(cL,i-1-cL,i)/((hL,i-1+hL,i)/2) –
> > S*Dl*(cL,i-cL,i+1)/((hL,i+hL,i+1)/2)
> > VL,N*dcL,N/dt =
> > = S*Dl*(cL,N-1-cL,N)/((hL,N-1+hL,N)/2) –
> > S*Dl*(cL,N-cL*)/(hL,i/2) ; i=1..N
> > but VL,1 = S*hL,1 and with this my system
> becomes
> > s :
> > dcL,1/dt = 0 –
> > 1; Dl*(cL,1-cL,2)/(hL,1*(hL,1+hL,2)/2)
> > dcL,i/dt = Dl*(cL,i-1-cL,i)/[hL,i
> > ,i *((hL,i-1+hL,i)/2) –
> > Dl*(cL,i-cL,i+1)/((hL,i+hL,i+1)/2)]
> > dcL,N/dt =
> Dl*(cL,N-1-cL,N)/[hL,N*((hL,N-1+hL,N)/2)
> > 2) – Dl*(cL,N-cL*)/(hL,N/2))] ; i=1..N
> >
> > gas
> > dcG,1/dt = Dg*(CG*-cG,1)/(hG,1*hG,1/2) –
> > Dg*(CG,1-cG,2)/(hG,1*(hG,1+hG,2)/2)
> > dcG,i/dt = Dg*(cG,i-1-cG,i)/[hG,i
> > /[hG,i *((hG,i-1+hG,i)/2) –
> > Dg*(cG,i-cG,i+1)/((hG,i+hG,i+1)/2)]
> > dcG,P/dt =
> Dg*(cG,P-1-cG,P)/[hG,P*((hG,P-1+hG,P)/2)
> > ) – Dg*(cG,P-cgm*)/(hG,P/2))] ; i=1..P
> > membrane
> > dcM,1/dt = Dm*(Cgm*-cM,1)/(hM,1*hM,1/2) –
> > Dm*(CM,1-cM,2)/(hM,1*(hM,1+hM,2)/2)
> > dcM,i/dt = Dm*(cM,i-1-cM,i)/[hM,i
> > /[hM,i *((hM,i-1+hM,i)/2) –
> > Dm*(cM,i-cM,i+1)/((hM,i+hM,i+1)/2)]
> > dcM,Q/dt =
> Dm*(cM,Q-1-cM,Q)/[hM,Q*((hM,Q-1+hM,Q)/2)
> > ) – Dm*(cM,Q-cext*)/(hM,Q/2))] ; i=1..Q
> >
> > discretise the equations of boundary layer
> > liquid-gas interface: Dl * dcl/dx =Dg* dcg/dx
>
> > cl * Kgl = cg , Kgl
> > cl * Kgl = cg , Kgl =partition coefficient
> > Dl*cL,N/(hL,N/2) +
> > +
> >
> Dg*cG,1/(hG,1/2)-(Dl*cL**/(hL,N/2)+Dg*cG**/(hG,1/2))=0
>
> > cL* * Kgl=cG*
> > gas -membrane interface :
> > S*Dg * dcg/dx =Smem*Dm*
> > S*Dg * dcg/dx =Smem*Dm* dcm/dx
> > cg = cm
> > S*Dg*cG,P/(hG,P/2) +
> > +
> >
> Smem*Dm*cM,1/(hM,1/2)-(S*Dg/(hG,P/2)+Smem*Dm/(hM,1/2))
>
> > *cgm*=0
> > cgm* =cgm*
> > membrane -air compartment superior
> > Smem*Dm * dcm/dx =
> > Smem*Dm * dcm/dx = S*Kext *(cmext –
> > mext – cair)
> > Cm = cmext ; cair = 0 ; Kext
> > = cmext ; cair = 0 ; Kext – mass transfer
> > coefficient superior compartment side
> > Smem*Dm*cM,Q/(hM,Q/2) +
> > S*Kext*cair-(Smem*Dm/(hM,Q/2)+S*Kext)*cext*=0
> >
> > <<< My code
> > ''==========================main
> > program========================
> > %% Assumption
> > % The equilibrium is reached at the interface at
> any
> > time
> > format short e
> > global N P Q Kgl Dl Dg Dm S Smem h poro Kext
> > %% Definition of parameters and properties of the
> > compounds
> > % General Parameters
> > %------gas
> constant-----------------------------
> > R=8.3145; % J/(mol*K)
> >
> >
> >
> >
> %------temperature-----------------------------------
>
> > --
> > Tc=37; %° C
> > T=Tc+273.15; % K
> > %------total interior and exterior pressure
> > e (Pa)----------------
> > Pt=1.01325e5 ;
> > % Geometric parameters
> > % notation : A - compound, B- water, C1-air
> > -air compart 1, C2-air compart 2
> > %------volume compartment 1----------------
> > L=14e-2; % length of the
> > ngth of the cell, m;
> > H1=3e-2; % total height
> > otal height compartment 1, m;
> > S=L^2; % cell surface,
> > ll surface, m2;
> >
> > %------volume compartment 2----------------
> > H2=1.5e-2; % total height
> > otal height compartment 2, m;
> > %-------liquid volume liquide compartment
> > ent 1--------------------
> > hl=3e-3; % liquid height, m
> > hg=3e-2-hl; % gas height in the
> > ight in the compartment 1, m;
> > Vl=60e-6; % liquid volume, m3
> > Vg=(H1*S)-Vl; % gas volume in the
> > lume in the compart 1; m3
> > % Physical properties of pur compounds
> > %-------molecular
> > lar weight------------------------------
> > MA=102.13e-3; % compound,
> > % compound, kg/mol
> > Mw=18e-3; % water,
> > % water, kg/mol
> > Mo=32e-3; % O2,
> > % O2, kg/mol
> > Mn=28e-3; % N2,
> > % N2, kg/mol
> > Mco2=44e-3; %C02,
> > %C02, kg/mol
> > Mair1=0.21*Mo+0.79*Mn ; % molecular
> > % molecular weight of the air compart1, kg/mol
> > Mair2=0.21*Mo+0.74+Mn+0.05*Mco2;%
> > 0.05*Mco2;% molecular weight of the air compart2,
> > kg/mol
> > %------density at
> 37°C----------------------------
> > rho_Aliq=799.936 ; %
> > % compound, kg/m3
> > rho_wliq=993.342 ; %
> > % water, kg/m3
> > rho_Avap=MA*1000*273.15/(22.4*T); %
> > 22.4*T); % compound vapors, kg/m3
> > rho_wvap=Mw*1000*273.15/(22.4*T); % water
> > ); % water vapor, kg/m3
> > rho_air1=Pt*Mair1/(R*T) ; % air
> > % air compart 1, kg/m3
> > rho_air2= Pt*Mair2/(R*T); % air
> > % air compart 2, kg/m3
> > %-------viscosity at
> > at 37°C------------------------------------------
> > u_Aliq=0.46339e-3; %
> Pa*s
> > u_wliq=0.692e-3; %
> Pa*s
> > TcA= 554.5 ; %
> > % critical temperature compound, K
> > PcA=3.473; %
> > % critical pressure compound, MPa
> > Tr=T/TcA ; %
> > % reduced temperature
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> u_Avap=(46.1*Tr^0.618-20.4*exp(-0.449*Tr)+19.4*exp(-4
>
> >
> .058*Tr)+1)/((2.173424*10^11)*TcA^(1/6)*(MA*10^3)^(-1/
>
> > 2)*(PcA*10^6)^(-2/3));
> > % gas viscosity by the Yoon Thodos method,
> > os method, see Perry page 508,
> > % Pa*s
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> u_air1=1.71*(10^(-6))*((T/273.15)^(3/2))*((273.15+110
>
> > .4)/(T+110.4)); % Pa*s
> > u_air2=u_air1;
> > % Membrane characteristics
> > tort=2.71; % tortuosity
> > poro=0.7; % porosity
> > eff=poro/tort;
> > dp=0.36e-6; % pore diameter, m
> > Dimem=4e-2; % membrane actif
> > brane actif diameter, m
> > zmem=110e-6; % membrane thickness,
> > thickness, m
> > Smem=pi*(Dimem/2)^2; % memmbrane surface,
> > ne surface, m2
> > % Mass transfer parameters
> > % Diffusion coefficient
> > % in water
> > phi=2.6; % constant for the
> solvent
> > Vm=128.78; % liquid molar volume of the
> > ume of the compound, cm3/mol
> > Msolv=Mw*1000 ;
> > usolv=u_wliq*1000;
> >
> >
> >
> > Dl=coeffdiffliquid(T,phi,Msolv,usolv,Vm)*1e-4;
> > )*1e-4; % m2/s , Wilke and Chang
> > % in air
> > Tb=102+273.15; % boiling
> > % boiling temperature, K
> > Vc= 341.5 ; % critical volume,
> > al volume, cm3/mol
> > Dg=coeffdiffair(T,Pt*1e-5,MA*1000,Tb,Vc);
> > 00,Tb,Vc); % m2/s
> > Dm=eff*Dg; % effective
> > effective diffusion coefficient, m2/s
> >
> > % %-------coefficient de partage pour
> > l'arome--------------------------
> > Kgl=0.0211 ;
> > % Air flow in the superior compartment
> > Gair=(H2*L^2)/60; % m3/s
> > Sflow=L*H2; % m2
> > w_air=Gair/Sflow; % air velocity, m/s
> > ocity, m/s
> > % MAss transfer coefficient calculation
> > yair2=1; % air molar fraction
> > r fraction in the superior compart
> > Mext=yair2*Mair2; % mix molecular
> > molecular weight, kg/mol
> > rho_ext=Mext*Pt/(R*T);% mix density,
> kg/m3
> > miu_ext=yair2*u_air2; % mix viscosity, Pa*s
> > sity, Pa*s
> > Re_ext=w_air*rho_ext*L/miu_ext; % reynolds
> > % reynolds number
> > Sc_ext=miu_ext/(rho_ext*Dg); %schmidt
> > %schmidt number
> > Sh_ext=0.664*(Re_ext)^(1/2)*(Sc_ext)^(1/3);
> > xt)^(1/3); % sherwood number
> > Kext=Sh_ext*Dg/L;
> > %% Mesh
> > % space/volume discretisation
> > hl=3e-3;
> > N=10;
> > dz1=hl/N ;
> > z1=0:dz1:hl;
> > hL=zeros(1,N);
> > for i=1:N
> > hL(i)=z1(i+1)-z1(i);
> > end
> > H1=3e-2;
> > P=100;
> > dz2=(H1-hl)/P ;
> > z2=hl:dz2:H1;
> > hG=zeros(1,P);
> > for k=1:P
> > hG(k)=z2(k+1)-z2(k);
> > end
> > Dm=(poro/tort)*Dg;
> > Q=5;
> > dz3=zmem/Q;
> > z3=H1:dz3:(H1+zmem);
> > hM=zeros(1,Q);
> > for p=1:Q
> > hM(p)=z3(p+1)-z3(p);
> > end
> > % heights array
> > h=vertcat(hL',0,0,hG',0,hM',0);
> > %% Mass matrix
> > M=eye(N+P+Q+4);
> > M(N+1,N+1)=0;
> > M(N+2,N+2)=0;
> > M(N+P+3,N+P+3)=0;
> > M(N+P+Q+4,N+P+Q+4)=0;
> > % Initial condition
> > cA0=13.3 ; % kg/m3
> > c0=zeros(N+P+Q+4,1);
> > for i=1:N
> > c0(i,1)=cA0;
> > end
> >
> c0(N+1,1)=(2*Dl*c0(N,1)/h(N)+2*Dg*c0(N+3,1)/h(N+3))/(2
>
> > *Dg/h(N+3)+2*Dl*Kgl/h(N)); % conc
> > equilibre liquide
> > c0(N+2,1)=Kgl*c0(N+1,1);
> > tspan=0::3600*48;
> >
> options=odeset('Mass',M,'MassSingular','yes','JConstan
>
> > t','on'); %,'JPattern','on');
> > %
> > %
> >
> 'Refine',10,'RelTol',1e-6,'AbsTol',1e-9,'NormControl',
>
> > 'on', ...
> > % 'Vectorized','on');
> > [t,c] = ode15s(@diff1,tspan,c0,options);
> > cliq=c(:,1:N);
> > cgaz=c(:,N+3:N+P+2);
> > cmem=c(:,N+P+4:N+P+Q+3);
> > cliq_eq=c(:,N+1);
> > cgaz_eq=c(:,N+2);
> > cmem_ext=c(:,N+P+Q+4);
> > ======function with the matrix coefficient,
> > J===========
> > function dcdt = diff1(t,c)
> > % Derivative function.
> > % I am trying to solve
> > % dcdt = J * c
> > global Dl Dg Dm S Smem Kgl N P Q h poro Kext
> > % Matrix that contains the system coefficients
> > nts (Jacobian)
> > J=zeros(N+P+Q+4,N+P+Q+4);
> > % in the liquid
> > J(1,1)=-(2*Dl/h(1))/(h(1)+h(2));
> > J(1,2)=(2*Dl/h(1))/(h(1)+h(2));
> > for i=2:N-1
> > J(i,i-1)=(2*Dl/h(i))/(h(i-1)+h(i));
> >
> >
> >
> >
> >
> >
> >
> J(i,i)=(-2)*(Dl/h(i))*(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)
>
> > ));
> > J(i,i+1)=(2*Dl/h(i))/(h(i)+h(i+1));
> > end
> > J(N,N-1)=(2*Dl/h(N))/(h(N-1)+h(N));
> > J(N,N)=(-2)*(Dl/h(N))*(1/(h(N-1)+h(N))+1/h(N));
> > J(N,N+1)=(2*Dl/h(N))/h(N);
> > % first interace AE
> > % flux continuity
> > J(N+1,N)=2*Dl/h(N);
> > J(N+1,N+1)=-2*Dl/h(N) ;
> > J(N+1,N+2)=(-2)*Dg/h(N+3);
> > J(N+1,N+3)=2*Dg/h(N+3);
> > % partition coefficient
> > J(N+2,N+1)=Kgl;
> > J(N+2,N+2)=-1;
> > % in the gas
> > J(N+3,N+2)=2*Dg/(h(N+3)*h(N+3));
> >
> >
> >
> J(N+3,N+3)=(-2)*(Dg/h(N+3))*(1/h(N+3)+1/(h(N+3)+h(N+4
>
> > )));
> > J(N+3,N+4)=(2*Dg/h(N+3))/(h(N+3)+h(N+4));
> > for i=N+4:N+P+2-1
> > J(i,i-1)=(2*Dg/h(i))/(h(i-1)+h(i));
> >
> >
> >
> >
> >
> >
> >
> J(i,i)=(-2)*(Dg/h(i))*(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)
>
> > ));
> > J(i,i+1)=(2*Dg/h(i))/(h(i)+h(i+1));
> > end
> >
> >
>
> J(N+P+2,N+P+1)=(2*Dg/h(N+P+2))/(h(N+P+1)+h(N+P+2));
> >
> >
>
>
> >
> J(N+P+2,N+P+2)=(-2)*(Dg/h(N+P+2))/(1/(h(N+P+1)+h(N+P+
>
> > 2))+1/h(N+P+2));
>
> Here is an error in your discretization.
> It must read:
> J(N+P+2,N+P+2)=(-2)*(Dg/h(N+P+2))*(1/(h(N+P+1)+h(N+P+
>
> 2))+1/h(N+P+2));
>
>
> > J(N+P+2,N+P+3)=2*Dg/(h(N+P+2)*h(N+P+2));
> > % gas membrane interface
> > % flux continuity and I assumed that the partition
> > ion coefficient is 1, so
> > % the concentrations are equals
> > J(N+P+3,N+P+2)=2*S*(Dg/h(N+P+2));
> >
> >
> >
> J(N+P+3,N+P+3)=-(2*S*Dg/h(N+P+2)+2*poro*Smem*Dm/h(N+P
>
> > +4));
> > J(N+P+3,N+P+4)=2*poro*Smem*Dm/h(N+P+4);
> > % in the membrane
> > J(N+P+4,N+P+3)=2*poro*Dm/(h(N+P+4)*h(N+P+4));
> >
> >
> >
> J(N+P+4,N+P+4)=-2*poro*(Dm/h(N+P+4))*(1/h(N+P+4)+1/(h
>
> > (N+P+4)+h(N+P+5)));
> >
> >
> >
> J(N+P+4,N+P+5)=2*poro*(Dm/h(N+P+4))/(h(N+P+4)+h(N+P+5
>
> > ));
> > for i=N+P+5:N+P+Q+3-1
> > J(i,i-1)=(2*poro*Dm/h(i))/(h(i-1)+h(i));
> >
> >
> >
> >
> >
> >
> >
> J(i,i)=(-2)*poro*(Dm/h(i))*(1/(h(i-1)+h(i))+1/(h(i)+h
>
> > (i+1)));
> > J(i,i+1)=(2*poro*Dm/h(i))/(h(i)+h(i+1));
> > end
> >
> >
> >
> J(N+P+Q+3,N+P+Q+3-1)=2*poro*(Dm/h(N+P+Q+3))/(h(N+P+Q+
>
> > 3-1)+h(N+P+Q+3));
> >
> >
> >
> J(N+P+Q+3,N+P+Q+3)=-2*poro*(Dm/h(N+P+Q+3))*(1/(h(N+P+
>
> > Q+3-1)+h(N+P+Q+3))+1/h(N+P+Q+3));
> >
> >
> >
> J(N+P+Q+3,N+P+Q+4)=2*poro*Dm/(h(N+P+Q+3)*h(N+P+Q+3));
>
> >
> > % interface membrane - air
> > % continuity flux and K = 1, also the
> > the concentration in the air is
> > % assumed 0 at all times (because we have chosen
> > sen the air velocity in
> > % such an way that the holl volume is renewed each
> > ach minute)
> > J(N+P+Q+4,N+P+Q+3)=2*poro*Smem*Dm/h(N+P+Q+3);
> >
> >
> >
> J(N+P+Q+4,N+P+Q+4)=-(2*poro*Smem*Dm/h(N+P+Q+3)+Kext*S
>
> > );
> > dcdt = J * c;
> > %cond(J)
> > end
> > ======diffusion coefficient functions======
> > ======function Dl = coeffdiffliquid(
> > T,phi,Msolv,usolv,Vm)
> > % Diffusion coefficient in a liquid , cm2/s
> > % Wilke and chang method
> > % T, K; phi,-; Msolv,g/mol; usolv, cP, Vm, cm3*/mol
>
> >
> >
> >
> Dl=(7.4*10^(-8)*((phi*Msolv)^(1/2))*T/(usolv*Vm^0.6))
>
> > ;
> > end
> > =======function Dg = coeffdiffair( T,Pt,MA,Tb,Vc)
> > % Function to determine the diff coeff in airr
> > % T,K ; Pt bar; M g/mol ; Tb k; Vc cm3/mol
> > % proprietes air
> > epsilon_kair=97; % k
> > sigma_air=3.62; % Angstrom
> > Mair=0.21*32+0.79*28 ;
> > % proprietes compound
> > Vb=0.285*(Vc^1.048);
> > epsilon_k=1.15*Tb;
> > sigma=1.18*(Vb)^(1/3);
> > % proprietes mix
> > epsilon_kmix=sqrt(epsilon_k*epsilon_kair);
> > sigma_mix=(sigma_air+sigma)/2;
> > Mmix=2/(1/MA+1/Mair);
> > Tx=T*epsilon_kmix^(-1); % potentiel Lennard-
> > d- Jones
> > % calculation collision integrale ;
> > A=1.06036;
> > B=0.15610;
> > C=0.19300;
> > D=0.47635;
> > F=1.52996;
> > E=1.03587;
> > G=1.76474;
> > H=3.89411;
> >
> >
> gammaD=A/(Tx*B)+C/exp(D*Tx)+E/exp(F*Tx)+G/exp(H*Tx);
> > % calculation binary diff coeff in a gas, using
> > % Wilke and Lee method
> >
> >
> Dg=10^(-4)*(3.03-(0.98/Mmix^(1/2)))*1e-3*T^(3/2)/(Pt*
>
> > Mmix^(1/2)*sigma_mix^2*gammaD); % en m2/s
> > end''
> > So I can't make it work properly, moreover I can
> see
> > that the problem with the negative concentrations
> it
> > is still there!
> > I read that a row equilibration may help
> stabilize
> > my system , but I don't know how to do it! I
> allready
> > did an adimensionalization but not in the version
> I
> > just post and even doing this I still have
> negative
> > values!
> > Best regards ,
> > A. STOIAN
From: Alina Stoian on 14 Jun 2010 08:52

Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <909946030.322787.1276346766130.JavaMail.root(a)gallium.mathforum.org>...
> Hello Alina,
>
> two more things I would like to mention:
>
> 1.
> I think the starting values for the concentrations at
> the liquid-air interface are not consistent with the
> algebraic equations you define in function dcdt:
>
> In my opinion, they should read:
> c0(N+1,1) = (2*Dl*c0(N,1)/h(N)+2*Dg*c0(N+3,1)/h(N+3))/
> (2*Dl/h(N)+2*Dg*Kgl/h(N+3))
> c0(N+2,1) = Kgl*c0(N+1,1)
>
> 2.
> Because of the storage requirements,
> I would not set up the Jacobian matrix in function dcdt.
> Instead, I would provide J*c directly.
> E.g. instead of
> J(1,1) = -(2*Dl/h(1))/(h(1)+h(2))
> J(1,2) = (2*Dl/h(1))/(h(1)+h(2)),
> I would prescribe
> dcdt(1) = -(2*Dl/h(1))/(h(1)+h(2))*c(1) +
> (2*Dl/h(1))/(h(1)+h(2))*c(2)
> and so on.
>
> Best wishes
> Torsten.
>
>
> > Hello Alina,
> > I found one error in your discretization (see below)
> > After the change, all went well with my solver.
> > Hope MATLAB's ODE solver now also gives satisfactory
> > results.
> >
> > Best wishes
> > Torsten.


Hello,

I don't know why it is not working for me too, but I have corrected like you told me to do and it continues to give negatives values for the concentration!
I've tried the two versions, with the defined matrix coefficients and with the defined equations and the results are identical but negatives!

My code it is :

===================main program===============================
clear all
clc

%% Assumption
% The equilibrium is reached at the interface at any time
format short e
global N P Q Kgl Dl Dg Dm S Smem h poro Kext
%% Definition of parameters and properties of the compounds
% General Parameters

%------gas constant-----------------------------
R=8.3145; % J/(mol*K)
%------temperature-------------------------------------
Tc=37; %° C
T=Tc+273.15; % K
%------total interior and exterior pressure (Pa)----------------
Pt=1.01325e5 ;
% Geometric parameters
% notation : A - compound, B- water, C1-air compart 1, C2-air compart 2
%------volume compartment 1----------------
L=14e-2; % length of the cell, m;
H1=3e-2; % total height compartment 1, m;
S=L^2; % cell surface, m2;
%------volume compartment 2----------------
H2=1.5e-2; % total height compartment 2, m;
%-------liquid volume liquide compartment 1--------------------
hl=3e-3; % liquid height, m
hg=3e-2-hl; % gas height in the compartment 1, m;
Vl=60e-6; % liquid volume, m3
Vg=(H1*S)-Vl; % gas volume in the compart 1; m3
% Physical properties of pur compounds
%-------molecular weight------------------------------
MA=102.13e-3; % compound, kg/mol
Mw=18e-3; % water, kg/mol
Mo=32e-3; % O2, kg/mol
Mn=28e-3; % N2, kg/mol
Mco2=44e-3; %C02, kg/mol
Mair1=0.21*Mo+0.79*Mn ; % molecular weight of the air compart1, kg/mol
Mair2=0.21*Mo+0.74+Mn+0.05*Mco2;% molecular weight of the air compart2, kg/mol
%------density at 37°C----------------------------
rho_Aliq=799.936 ; % compound, kg/m3
rho_wliq=993.342 ; % water, kg/m3
rho_Avap=MA*1000*273.15/(22.4*T); % compound vapors, kg/m3
rho_wvap=Mw*1000*273.15/(22.4*T); % water vapor, kg/m3
rho_air1=Pt*Mair1/(R*T) ; % air compart 1, kg/m3
rho_aanir2= Pt*Mair2/(R*T); % air compart 2, kg/m3
%-------viscosity at 37°C------------------------------------------
u_Aliq=0.46339e-3; % Pa*s
u_wliq=0.692e-3; % Pa*s
TcA= 554.5 ; % critical temperature compound, K
PcA=3.473; % critical pressure compound, MPa
Tr=T/TcA ; % reduced temperature
u_Avap=(46.1*Tr^0.618-20.4*exp(-0.449*Tr)+19.4*exp(-4.058*Tr)+1)/((2.173424*10^11)*TcA^(1/6)*(MA*10^3)^(-1/2)*(PcA*10^6)^(-2/3));
% gas viscosity by the Yoon Thodos method, see Perry page 508,
% Pa*s
u_air1=1.71*(10^(-6))*((T/273.15)^(3/2))*((273.15+110.4)/(T+110.4)); % Pa*s
u_air2=u_air1;
% Membrane characteristics
tort=2.71; % tortuosity
poro=0.7; % porosity
eff=poro/tort;
dp=0.36e-6; % pore diameter, m
Dimem=4e-2; % membrane actif diameter, m
zmem=110e-6; % membrane thickness, m
Smem=pi*(Dimem/2)^2; % memmbrane surface, m2
% Mass transfer parameters
% Diffusion coefficient
% in water
phi=2.6; % constant for the solvent
Vm=128.78; % liquid molar volume of the compound, cm3/mol
Msolv=Mw*1000 ;
usolv=u_wliq*1000;
Dl=coeffdiffliquid(T,phi,Msolv,usolv,Vm)*1e-4; % m2/s , Wilke and Chang
% in air
Tb=102+273.15; % boiling temperature, K
Vc= 341.5 ; % critical volume, cm3/mol
Dg=coeffdiffair(T,Pt*1e-5,MA*1000,Tb,Vc); % m2/s
Dm=eff*Dg; % effective diffusion coefficient, m2/s
%-------partition coefficient for the compound--------------------------
Kgl=0.0211 ;

% Air flow in the superior compartment
Gair=(H2*L^2)/60; % m3/s
Sflow=L*H2; % m2
w_air=Gair/Sflow; % air velocity, m/s

% MAss transfer coefficient calculation
yair2=1; % air molar fraction in the superior compart
Mext=yair2*Mair2; % mix molecular weight, kg/mol
rho_ext=Mext*Pt/(R*T);% mix density, kg/m3
miu_ext=yair2*u_air2; % mix viscosity, Pa*s

Re_ext=w_air*rho_ext*L/miu_ext; % reynolds number
Sc_ext=miu_ext/(rho_ext*Dg); %schmidt number
Sh_ext=0.664*(Re_ext)^(1/2)*(Sc_ext)^(1/3); % sherwood number
Kext=Sh_ext*Dg/L;

%% Mesh

% space/volume discretisation
hl=3e-3;
N=10;
dz1=hl/N ;
z1=0:dz1:hl;
hL=zeros(1,N);
for i=1:N
hL(i)=z1(i+1)-z1(i);
end

H1=3e-2;
P=100;
dz2=(H1-hl)/P ;
z2=hl:dz2:H1;
hG=zeros(1,P);
for k=1:P
hG(k)=z2(k+1)-z2(k);
end

Dm=(poro/tort)*Dg;
Q=5;
dz3=zmem/Q;
z3=H1:dz3:(H1+zmem);
hM=zeros(1,Q);
for p=1:Q
hM(p)=z3(p+1)-z3(p);
end

% heights array

h=vertcat(hL',0,0,hG',0,hM',0);

%% Mass matrix
M=eye(N+P+Q+4);
M(N+1,N+1)=0;
M(N+2,N+2)=0;
M(N+P+3,N+P+3)=0;
M(N+P+Q+4,N+P+Q+4)=0;

% Initial condition
cA0=13.3 ; % kg/m3

c0=zeros(1,N+P+Q+4);
for i=1:N
c0(1,i)=cA0;
end
c0(1,N+1)=(Dl*c0(1,N)/h(N)+Dg*c0(1,N+3)/h(N+3))/(Dg*Kgl/h(N+3)+Dl/h(N)); % conc equilibre liquide
c0(1,N+2)=Kgl*c0(1,N+1);

tspan=0:3600:3600*48;
rtol=1e-6;
atol=zeros(N+P+Q+4,1);
atol(1:N)=rtol;
atol(N+1:N+P+Q+4)=rtol*1e-3;

options=odeset('Mass',M,'Jconstant','on');

[t,c] = ode15s(@diff11t,tspan,c0,options);

cliq=c(:,1:N);
cgaz=c(:,N+3:N+P+2);
cmem=c(:,N+P+4:N+P+Q+3);
cliq_eq=c(:,N+1);
cgaz_eq=c(:,N+2);
cmem_ext=c(:,N+P+Q+4);
==================================================================
===========function, version with the defined coefficients===================
function dcdt = diff1(t,c)

% Derivative function.
% I am trying to solve
% dcdt = A * c
global Dl Dg Dm S Smem Kgl N P Q h poro Kext
% Matrix that contains the system coefficients

A=zeros(N+P+Q+4,N+P+Q+4);

% in the liquid

A(1,1)=-(2*Dl/h(1))/(h(1)+h(2));
A(1,2)=(2*Dl/h(1))/(h(1)+h(2));
for i=2:N-1
A(i,i-1)=(2*Dl/h(i))/(h(i-1)+h(i));
A(i,i)=(-2)*(Dl/h(i))*(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)));
A(i,i+1)=(2*Dl/h(i))/(h(i)+h(i+1));
end
A(N,N-1)=(2*Dl/h(N))/(h(N-1)+h(N));
A(N,N)=(-2)*(Dl/h(N))*(1/(h(N-1)+h(N))+1/h(N));
A(N,N+1)=(2*Dl/h(N))/h(N);

% first interace AE
% flux continuity
A(N+1,N)=2*Dl/h(N);
A(N+1,N+1)=-2*(Dl/h(N)+Dg*Kgl/h(N+3)) ;
% change
A(N+1,N+3)=2*Dg/h(N+3);
% partition coefficient
A(N+2,N+1)=Kgl;
A(N+2,N+2)=-1;

% in the gas

A(N+3,N+2)=2*Dg/(h(N+3)*h(N+3));
A(N+3,N+3)=(-2)*(Dg/h(N+3))*(1/h(N+3)+1/(h(N+3)+h(N+4)));
A(N+3,N+4)=(2*Dg/h(N+3))/(h(N+3)+h(N+4));
for i=N+4:N+P+2-1
A(i,i-1)=(2*Dg/h(i))/(h(i-1)+h(i));
A(i,i)=(-2)*(Dg/h(i))*(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)));
A(i,i+1)=(2*Dg/h(i))/(h(i)+h(i+1));
end
A(N+P+2,N+P+1)=(2*Dg/h(N+P+2))/(h(N+P+1)+h(N+P+2));
A(N+P+2,N+P+2)=(-2)*(Dg/h(N+P+2))*(1/(h(N+P+1)+h(N+P+2))+1/h(N+P+2));
A(N+P+2,N+P+3)=2*Dg/(h(N+P+2)*h(N+P+2));
% gas membrane interface
% flux continuity and I assumed that the partition coefficient is 1, so
% the concentrations are equals

A(N+P+3,N+P+2)=2*S*(Dg/h(N+P+2));
A(N+P+3,N+P+3)=-(2*S*Dg/h(N+P+2)+2*poro*Smem*Dm/h(N+P+4));
A(N+P+3,N+P+4)=2*poro*Smem*Dm/h(N+P+4);

% in the membrane

A(N+P+4,N+P+3)=2*poro*Dm/(h(N+P+4)*h(N+P+4));
A(N+P+4,N+P+4)=-2*poro*(Dm/h(N+P+4))*(1/h(N+P+4)+1/(h(N+P+4)+h(N+P+5)));
A(N+P+4,N+P+5)=2*poro*(Dm/h(N+P+4))/(h(N+P+4)+h(N+P+5));
for i=N+P+5:N+P+Q+3-1
A(i,i-1)=(2*poro*Dm/h(i))/(h(i-1)+h(i));
A(i,i)=(-2)*poro*(Dm/h(i))*(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)));
A(i,i+1)=(2*poro*Dm/h(i))/(h(i)+h(i+1));
end

A(N+P+Q+3,N+P+Q+3-1)=2*poro*(Dm/h(N+P+Q+3))/(h(N+P+Q+3-1)+h(N+P+Q+3));
A(N+P+Q+3,N+P+Q+3)=-2*poro*(Dm/h(N+P+Q+3))*(1/(h(N+P+Q+3-1)+h(N+P+Q+3))+1/h(N+P+Q+3));
A(N+P+Q+3,N+P+Q+4)=2*poro*Dm/(h(N+P+Q+3)*h(N+P+Q+3));

% interface membrane - air
% continuity flux and K = 1, also the concentration in the air is
% assumed 0 at all times (because we have chosen the air velocity in
% such an way that the holl volume is renewed each minute)
A(N+P+Q+4,N+P+Q+3)=2*poro*Smem*Dm/h(N+P+Q+3);
A(N+P+Q+4,N+P+Q+4)=-(2*poro*Smem*Dm/h(N+P+Q+3)+Kext*S);

dcdt = A * c;
end

=====================function with the equations defined==================

function dcdt = diff11t(t,c)
% Derivative function.
global Dl Dg Dm S Smem Kgl N P Q h poro Kext
% System of Differential equations
dcdt=zeros(N+P+Q+4,1);
% in the liquid
dcdt(1)=-(2*Dl/h(1))*(c(1)-c(2))/(h(1)+h(2));
for i=2:N-1
dcdt(i)=(2*Dl/h(i))*(c(i-1)/(h(i-1)+h(i))-...
(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)))*c(i)+...
c(i+1)/(h(i)+h(i+1)));
end
dcdt(N)=(2*Dl/h(N))*(c(N-1)/(h(N-1)+h(N))-...
(1/(h(N-1)+h(N))+1/h(N))*c(N)+...
c(N+1)/h(N));
% first interace AE
% flux continuity
dcdt(N+1)=(Dl/h(N))*c(N)+(Dg/h(N+3))*c(N+3)-...
(Dl/h(N)+Dg*Kgl/h(N+3))*c(N+1);
% partition coefficient
dcdt(N+2)=Kgl*c(N+1)-c(N+2);
% in the gas
dcdt(N+3)=(2*Dg/h(N+3))*(c(N+2)/h(N+3)-...
(1/h(N+3)+1/(h(N+3)+h(N+4)))*c(N+3)+...
c(N+4)/(h(N+3)+h(N+4)));
for i=N+4:N+P+2-1
dcdt(i)=(2*Dg/h(i))*(c(i-1)/(h(i-1)+h(i))-...
(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)))*c(i)+...
c(i+1)/(h(i)+h(i+1)));
end
dcdt(N+P+2)=(2*Dg/h(N+P+2))*(c(N+P+1)/(h(N+P+1)+h(N+P+2))-...
(1/(h(N+P+1)+h(N+P+2))+1/h(N+P+2))*c(N+P+2)+...
c(N+P+3)/h(N+P+2));
% gas membrane interface
% flux continuity and I assumed that the partition coefficient is 1, so
% the concentrations are equals
dcdt(N+P+3)=2*S*(Dg/h(N+P+2))*c(N+P+2)-...
(2*S*Dg/h(N+P+2)+2*poro*Smem*Dm/h(N+P+4))*c(N+P+3)+...
2*poro*Smem*(Dm/h(N+P+4))*c(N+P+4);
% in the membrane
dcdt(N+P+4)=(2*poro*Dm/(h(N+P+4)))*(c(N+P+3)/h(N+P+4)-...
(1/h(N+P+4)+1/(h(N+P+4)+h(N+P+5)))*c(N+P+4)+...
c(N+P+5)/(h(N+P+4)+h(N+P+5)));
for i=N+P+5:N+P+Q+3-1
dcdt(i)=(2*poro*Dm/h(i))*(c(i-1)/(h(i-1)+h(i))-...
(1/(h(i-1)+h(i))+1/(h(i)+h(i+1)))*c(i)+...
c(i+1)/(h(i)+h(i+1)));
end
dcdt(N+P+Q+3)=2*poro*(Dm/h(N+P+Q+3))*(c(N+P+Q+2)/(h(N+P+Q+3-1)+h(N+P+Q+3))-...
(1/(h(N+P+Q+2)+h(N+P+Q+3))+1/h(N+P+Q+3))*c(N+P+Q+3)+...
c(N+P+Q+4)/h(N+P+Q+3));
% interface membrane - air
% continuity flux and K = 1, also the concentration in the air is
% assumed 0 at all times (because we have chosen the air velocity in
% such an way that the holl volume is renewed each minute)
dcdt(N+P+Q+4)=2*poro*Smem*Dm*(1/h(N+P+Q+3))*c(N+P+Q+3)-...
(2*poro*Smem*Dm/h(N+P+Q+3)+Kext*S)*c(N+P+Q+4);
end

========functions for the coeff diff calculation==========
function Dl = coeffdiffliquid( T,phi,Msolv,usolv,Vm)
% Diffusion coefficient in a liquid , cm2/s
% Wilke and chang method

% T, K; phi,-; Msolv,g/mol; usolv, cP, Vm, cm3*/mol

Dl=(7.4*10^(-8)*((phi*Msolv)^(1/2))*T/(usolv*Vm^0.6));

end
==========diff coeff air==================
function Dg = coeffdiffair( T,Pt,MA,Tb,Vc)
% Fonction pour calculer le coefficient de diffusion dans l'air
% T,K ; Pt bar; M g/mol ; Tb k; Vc cm3/mol

% proprietes air
epsilon_kair=97; % k
sigma_air=3.62; % Angstrom
Mair=0.21*32+0.79*28 ;

% proprietes composé
Vb=0.285*(Vc^1.048);
epsilon_k=1.15*Tb;
sigma=1.18*(Vb)^(1/3);
% proprietes melange
epsilon_kmix=sqrt(epsilon_k*epsilon_kair);
sigma_mix=(sigma_air+sigma)/2;
Mmix=2/(1/MA+1/Mair);
Tx=T*epsilon_kmix^(-1); % potentiel Lennard- Jones

% calcul integrale de collision;
A=1.06036;
B=0.15610;
C=0.19300;
D=0.47635;
F=1.52996;
E=1.03587;
G=1.76474;
H=3.89411;

gammaD=A/(Tx*B)+C/exp(D*Tx)+E/exp(F*Tx)+G/exp(H*Tx);

% calcul coeff diffusion binaire dans un gas, utilisant la methode de
% Wilke and Lee

Dg=10^(-4)*(3.03-(0.98/Mmix^(1/2)))*1e-3*T^(3/2)/(Pt*Mmix^(1/2)*sigma_mix^2*gammaD); % en m2/s
end
===============================================================
Here you have some results that I obtain:
for cliq using the function with the defined coefficients
======= =======================================================
cliq =

Columns 1 through 7

1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01
8.3106e+00 8.2076e+00 8.0028e+00 7.6987e+00 7.2989e+00 6.8082e+00 6.2325e+00
4.5059e+00 4.4497e+00 4.3382e+00 4.1726e+00 3.9550e+00 3.6882e+00 3.3754e+00
2.4425e+00 2.4121e+00 2.3516e+00 2.2618e+00 2.1438e+00 1.9992e+00 1.8296e+00
1.3247e+00 1.3082e+00 1.2754e+00 1.2267e+00 1.1627e+00 1.0843e+00 9.9229e-01
7.1826e-01 7.0931e-01 6.9152e-01 6.6512e-01 6.3043e-01 5.8789e-01 5.3802e-01
3.8936e-01 3.8451e-01 3.7487e-01 3.6056e-01 3.4175e-01 3.1869e-01 2.9165e-01
2.1114e-01 2.0851e-01 2.0328e-01 1.9552e-01 1.8532e-01 1.7281e-01 1.5815e-01
1.1437e-01 1.1295e-01 1.1012e-01 1.0591e-01 1.0039e-01 9.3612e-02 8.5671e-02
6.1928e-02 6.1156e-02 5.9622e-02 5.7346e-02 5.4355e-02 5.0687e-02 4.6387e-02
3.3531e-02 3.3113e-02 3.2283e-02 3.1050e-02 2.9431e-02 2.7444e-02 2.5116e-02
1.8139e-02 1.7913e-02 1.7464e-02 1.6797e-02 1.5921e-02 1.4847e-02 1.3587e-02
9.8194e-03 9.6970e-03 9.4539e-03 9.0929e-03 8.6186e-03 8.0370e-03 7.3552e-03
5.3145e-03 5.2483e-03 5.1167e-03 4.9213e-03 4.6647e-03 4.3499e-03 3.9809e-03
2.8782e-03 2.8424e-03 2.7711e-03 2.6653e-03 2.5263e-03 2.3558e-03 2.1560e-03
1.5610e-03 1.5415e-03 1.5029e-03 1.4455e-03 1.3701e-03 1.2776e-03 1.1693e-03
8.4564e-04 8.3511e-04 8.1416e-04 7.8308e-04 7.4224e-04 6.9215e-04 6.3343e-04
4.5858e-04 4.5287e-04 4.4151e-04 4.2466e-04 4.0251e-04 3.7534e-04 3.4350e-04
2.4848e-04 2.4539e-04 2.3923e-04 2.3010e-04 2.1810e-04 2.0338e-04 1.8613e-04
1.3437e-04 1.3270e-04 1.2937e-04 1.2443e-04 1.1794e-04 1.0998e-04 1.0065e-04
7.2750e-05 7.1843e-05 7.0042e-05 6.7367e-05 6.3854e-05 5.9544e-05 5.4493e-05
3.9243e-05 3.8754e-05 3.7782e-05 3.6339e-05 3.4444e-05 3.2120e-05 2.9395e-05
2.0902e-05 2.0641e-05 2.0124e-05 1.9355e-05 1.8346e-05 1.7108e-05 1.5657e-05
1.1110e-05 1.0972e-05 1.0696e-05 1.0288e-05 9.7512e-06 9.0930e-06 8.3214e-06
5.9981e-06 5.9233e-06 5.7747e-06 5.5542e-06 5.2644e-06 4.9090e-06 4.4924e-06
3.2636e-06 3.2230e-06 3.1421e-06 3.0221e-06 2.8645e-06 2.6712e-06 2.4446e-06
1.7270e-06 1.7054e-06 1.6627e-06 1.5992e-06 1.5158e-06 1.4135e-06 1.2936e-06
7.1937e-07 7.1041e-07 6.9259e-07 6.6614e-07 6.3138e-07 5.8874e-07 5.3876e-07
1.8669e-07 1.8436e-07 1.7971e-07 1.7282e-07 1.6376e-07 1.5266e-07 1.3965e-07
-1.9395e-08 -1.9177e-08 -1.8743e-08 -1.8101e-08 -1.7259e-08 -1.6228e-08 -1.5024e-08
-9.3559e-08 -9.2426e-08 -9.0173e-08 -8.6830e-08 -8.2439e-08 -7.7056e-08 -7.0752e-08
-2.0588e-08 -2.0360e-08 -1.9904e-08 -1.9228e-08 -1.8337e-08 -1.7240e-08 -1.5949e-08
3.2893e-08 3.2470e-08 3.1631e-08 3.0386e-08 2.8753e-08 2.6754e-08 2.4415e-08
5.2068e-08 5.1422e-08 5.0138e-08 4.8234e-08 4.5733e-08 4.2666e-08 3.9071e-08
2.8747e-08 2.8395e-08 2.7694e-08 2.6654e-08 2.5286e-08 2.3605e-08 2.1632e-08
-2.3290e-09 -2.2970e-09 -2.2337e-09 -2.1404e-09 -2.0189e-09 -1.8716e-09 -1.7013e-09
-1.1730e-07 -1.1584e-07 -1.1295e-07 -1.0866e-07 -1.0303e-07 -9.6120e-08 -8.8032e-08
-2.6730e-07 -2.6399e-07 -2.5740e-07 -2.4763e-07 -2.3480e-07 -2.1907e-07 -2.0064e-07
-3.4790e-07 -3.4359e-07 -3.3502e-07 -3.2229e-07 -3.0558e-07 -2.8510e-07 -2.6109e-07
-2.6957e-07 -2.6622e-07 -2.5957e-07 -2.4969e-07 -2.3670e-07 -2.2077e-07 -2.0210e-07
-3.3115e-07 -3.2704e-07 -3.1888e-07 -3.0676e-07 -2.9084e-07 -2.7131e-07 -2.4841e-07
-3.4598e-07 -3.4170e-07 -3.3319e-07 -3.2057e-07 -3.0397e-07 -2.8361e-07 -2.5974e-07
-2.8122e-07 -2.7776e-07 -2.7087e-07 -2.6065e-07 -2.4722e-07 -2.3075e-07 -2.1145e-07
-1.5148e-07 -1.4964e-07 -1.4599e-07 -1.4056e-07 -1.3343e-07 -1.2468e-07 -1.1444e-07
-7.3131e-08 -7.2271e-08 -7.0562e-08 -6.8026e-08 -6.4698e-08 -6.0621e-08 -5.5850e-08
4.8581e-10 4.3055e-10 3.2024e-10 1.5531e-10 -6.3779e-11 -3.3656e-10 -6.6288e-10
4.9693e-08 4.9035e-08 4.7728e-08 4.5787e-08 4.3237e-08 4.0110e-08 3.6445e-08
6.3442e-08 6.2631e-08 6.1019e-08 5.8628e-08 5.5490e-08 5.1647e-08 4.7151e-08
4.6635e-08 4.6052e-08 4.4893e-08 4.3173e-08 4.0915e-08 3.8147e-08 3.4906e-08

Columns 8 through 10

1.3300e+01 1.3300e+01 1.3300e+01
5.5787e+00 4.8550e+00 4.0702e+00
3.0205e+00 2.6280e+00 2.2027e+00
1.6372e+00 1.4244e+00 1.1939e+00
8.8795e-01 7.7255e-01 6.4752e-01
4.8145e-01 4.1888e-01 3.5109e-01
2.6099e-01 2.2707e-01 1.9032e-01
1.4153e-01 1.2313e-01 1.0320e-01
7.6663e-02 6.6700e-02 5.5905e-02
4.1510e-02 3.6115e-02 3.0270e-02
2.2475e-02 1.9554e-02 1.6390e-02
1.2159e-02 1.0578e-02 8.8664e-03
6.5818e-03 5.7264e-03 4.7996e-03
3.5623e-03 3.0993e-03 2.5977e-03
1.9293e-03 1.6785e-03 1.4069e-03
1.0463e-03 9.1032e-04 7.6300e-04
5.6683e-04 4.9316e-04 4.1335e-04
3.0739e-04 2.6744e-04 2.2416e-04
1.6656e-04 1.4491e-04 1.2146e-04
9.0069e-05 7.8363e-05 6.5681e-05
4.8763e-05 4.2426e-05 3.5559e-05
2.6304e-05 2.2885e-05 1.9182e-05
1.4010e-05 1.2189e-05 1.0216e-05
7.4461e-06 6.4780e-06 5.4291e-06
4.0198e-06 3.4972e-06 2.9309e-06
2.1875e-06 1.9031e-06 1.5951e-06
1.1575e-06 1.0071e-06 8.4404e-07
4.8206e-07 4.1934e-07 3.5139e-07
1.2490e-07 1.0859e-07 9.0927e-08
-1.3666e-08 -1.2175e-08 -1.0577e-08
-6.3607e-08 -5.5716e-08 -4.7182e-08
-1.4474e-08 -1.2826e-08 -1.1017e-08
2.1770e-08 1.8853e-08 1.5706e-08
3.4993e-08 3.0479e-08 2.5581e-08
1.9391e-08 1.6906e-08 1.4210e-08
-1.5105e-09 -1.3019e-09 -1.0774e-09
-7.8860e-08 -6.8716e-08 -5.7717e-08
-1.7974e-07 -1.5663e-07 -1.3156e-07
-2.3386e-07 -2.0373e-07 -1.7106e-07
-1.8090e-07 -1.5743e-07 -1.3200e-07
-2.2241e-07 -1.9365e-07 -1.6246e-07
-2.3265e-07 -2.0267e-07 -1.7017e-07
-1.8954e-07 -1.6531e-07 -1.3905e-07
-1.0284e-07 -9.0011e-08 -7.6133e-08
-5.0449e-08 -4.4494e-08 -3.8065e-08
-1.0429e-09 -1.4772e-09 -1.9665e-09
3.2289e-08 2.7693e-08 2.2716e-08
4.2064e-08 3.6456e-08 3.0406e-08
3.1232e-08 2.7171e-08 2.2774e-08 ====
====and for the function with the defined equations===============
cliq =

Columns 1 through 7

1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01 1.3300e+01
8.3106e+00 8.2076e+00 8.0029e+00 7.6988e+00 7.2990e+00 6.8082e+00 6.2325e+00
4.5059e+00 4.4498e+00 4.3382e+00 4.1726e+00 3.9551e+00 3.6882e+00 3.3754e+00
2.4425e+00 2.4120e+00 2.3516e+00 2.2618e+00 2.1438e+00 1.9991e+00 1.8295e+00
1.3247e+00 1.3081e+00 1.2753e+00 1.2266e+00 1.1627e+00 1.0842e+00 9.9222e-01
7.1819e-01 7.0924e-01 6.9145e-01 6.6505e-01 6.3036e-01 5.8782e-01 5.3795e-01
3.8931e-01 3.8446e-01 3.7482e-01 3.6051e-01 3.4170e-01 3.1864e-01 2.9161e-01
2.1110e-01 2.0847e-01 2.0324e-01 1.9548e-01 1.8529e-01 1.7278e-01 1.5812e-01
1.1435e-01 1.1293e-01 1.1010e-01 1.0589e-01 1.0037e-01 9.3597e-02 8.5658e-02
6.1920e-02 6.1149e-02 5.9615e-02 5.7339e-02 5.4349e-02 5.0681e-02 4.6382e-02
3.3527e-02 3.3109e-02 3.2279e-02 3.1047e-02 2.9428e-02 2.7442e-02 2.5114e-02
1.8138e-02 1.7912e-02 1.7463e-02 1.6796e-02 1.5920e-02 1.4846e-02 1.3587e-02
9.8192e-03 9.6969e-03 9.4537e-03 9.0928e-03 8.6186e-03 8.0370e-03 7.3553e-03
5.3147e-03 5.2485e-03 5.1169e-03 4.9215e-03 4.6648e-03 4.3501e-03 3.9811e-03
2.8763e-03 2.8405e-03 2.7693e-03 2.6635e-03 2.5246e-03 2.3543e-03 2.1546e-03
1.5571e-03 1.5377e-03 1.4991e-03 1.4419e-03 1.3667e-03 1.2745e-03 1.1664e-03
8.4271e-04 8.3221e-04 8.1134e-04 7.8037e-04 7.3967e-04 6.8976e-04 6.3125e-04
4.5614e-04 4.5046e-04 4.3917e-04 4.2240e-04 4.0037e-04 3.7335e-04 3.4169e-04
2.4683e-04 2.4376e-04 2.3765e-04 2.2857e-04 2.1665e-04 2.0203e-04 1.8489e-04
1.3353e-04 1.3187e-04 1.2856e-04 1.2365e-04 1.1721e-04 1.0930e-04 1.0003e-04
7.2259e-05 7.1359e-05 6.9569e-05 6.6913e-05 6.3424e-05 5.9144e-05 5.4128e-05
3.9012e-05 3.8526e-05 3.7560e-05 3.6126e-05 3.4242e-05 3.1931e-05 2.9223e-05
2.0876e-05 2.0616e-05 2.0099e-05 1.9332e-05 1.8324e-05 1.7088e-05 1.5639e-05
1.1182e-05 1.1042e-05 1.0766e-05 1.0356e-05 9.8164e-06 9.1551e-06 8.3800e-06
6.0556e-06 5.9803e-06 5.8306e-06 5.6085e-06 5.3166e-06 4.9586e-06 4.5389e-06
3.2787e-06 3.2379e-06 3.1567e-06 3.0362e-06 2.8779e-06 2.6837e-06 2.4562e-06
1.6932e-06 1.6721e-06 1.6302e-06 1.5679e-06 1.4861e-06 1.3859e-06 1.2685e-06
7.6411e-07 7.5460e-07 7.3570e-07 7.0766e-07 6.7084e-07 6.2570e-07 5.7281e-07
2.0090e-07 1.9851e-07 1.9376e-07 1.8671e-07 1.7744e-07 1.6606e-07 1.5271e-07
2.8061e-08 2.7939e-08 2.7699e-08 2.7348e-08 2.6899e-08 2.6373e-08 2.5799e-08
-2.2149e-09 -1.9222e-09 -1.3410e-09 -4.8002e-10 6.4782e-10 2.0251e-09 3.6294e-09
4.0430e-08 4.0124e-08 3.9511e-08 3.8592e-08 3.7365e-08 3.5826e-08 3.3963e-08
6.2005e-08 6.1269e-08 5.9802e-08 5.7615e-08 5.4726e-08 5.1158e-08 4.6943e-08
4.9604e-08 4.8943e-08 4.7631e-08 4.5683e-08 4.3129e-08 4.0005e-08 3.6356e-08
1.3059e-08 1.2850e-08 1.2438e-08 1.1835e-08 1.1060e-08 1.0136e-08 9.0905e-09
-6.0103e-08 -5.9339e-08 -5.7814e-08 -5.5537e-08 -5.2520e-08 -4.8782e-08 -4.4350e-08
-2.1118e-07 -2.0846e-07 -2.0305e-07 -1.9497e-07 -1.8430e-07 -1.7112e-07 -1.5554e-07
-3.1537e-07 -3.1130e-07 -3.0320e-07 -2.9114e-07 -2.7523e-07 -2.5561e-07 -2.3247e-07
-2.6701e-07 -2.6359e-07 -2.5678e-07 -2.4668e-07 -2.3341e-07 -2.1714e-07 -1.9810e-07
-3.0575e-07 -3.0179e-07 -2.9390e-07 -2.8221e-07 -2.6688e-07 -2.4811e-07 -2.2617e-07
-3.2649e-07 -3.2217e-07 -3.1360e-07 -3.0089e-07 -2.8421e-07 -2.6379e-07 -2.3992e-07
-2.6485e-07 -2.6123e-07 -2.5404e-07 -2.4338e-07 -2.2938e-07 -2.1222e-07 -1.9212e-07
-1.3340e-07 -1.3138e-07 -1.2735e-07 -1.2138e-07 -1.1351e-07 -1.0384e-07 -9.2461e-08
-5.3216e-08 -5.2158e-08 -5.0052e-08 -4.6916e-08 -4.2776e-08 -3.7663e-08 -3.1611e-08
2.1812e-08 2.1926e-08 2.2154e-08 2.2500e-08 2.2970e-08 2.3576e-08 2.4332e-08
7.0395e-08 6.9830e-08 6.8705e-08 6.7032e-08 6.4830e-08 6.2120e-08 5.8932e-08
1.0100e-07 9.9954e-08 9.7873e-08 9.4771e-08 9.0674e-08 8.5612e-08 7.9617e-08
1.0665e-07 1.0543e-07 1.0299e-07 9.9369e-08 9.4582e-08 8.8669e-08 8.1669e-08
7.4125e-08 7.3208e-08 7.1383e-08 6.8671e-08 6.5102e-08 6.0715e-08 5.5562e-08

Columns 8 through 10

1.3300e+01 1.3300e+01 1.3300e+01
5.5787e+00 4.8550e+00 4.0702e+00
3.0205e+00 2.6280e+00 2.2027e+00
1.6372e+00 1.4244e+00 1.1939e+00
8.8787e-01 7.7247e-01 6.4744e-01
4.8138e-01 4.1882e-01 3.5103e-01
2.6095e-01 2.2703e-01 1.9029e-01
1.4149e-01 1.2310e-01 1.0318e-01
7.6651e-02 6.6689e-02 5.5897e-02
4.1505e-02 3.6111e-02 3.0267e-02
2.2473e-02 1.9553e-02 1.6388e-02
1.2158e-02 1.0578e-02 8.8662e-03
6.5819e-03 5.7266e-03 4.7999e-03
3.5625e-03 3.0995e-03 2.5980e-03
1.9281e-03 1.6775e-03 1.4061e-03
1.0437e-03 9.0809e-04 7.6114e-04
5.6488e-04 4.9148e-04 4.1195e-04
3.0576e-04 2.6603e-04 2.2298e-04
1.6545e-04 1.4395e-04 1.2065e-04
8.9512e-05 7.7881e-05 6.5279e-05
4.8437e-05 4.2142e-05 3.5323e-05
2.6151e-05 2.2753e-05 1.9071e-05
1.3995e-05 1.2177e-05 1.0207e-05
7.5008e-06 6.5285e-06 5.4751e-06
4.0626e-06 3.5358e-06 2.9650e-06
2.1981e-06 1.9126e-06 1.6034e-06
1.1354e-06 9.8826e-07 8.2891e-07
5.1283e-07 4.4650e-07 3.7461e-07
1.3753e-07 1.2071e-07 1.0244e-07
2.5220e-08 2.4693e-08 2.4297e-08
5.4326e-09 7.3985e-09 9.4804e-09
3.1758e-08 2.9178e-08 2.6170e-08
4.2125e-08 3.6761e-08 3.0930e-08
3.2239e-08 2.7719e-08 2.2873e-08
7.9533e-09 6.7535e-09 5.5139e-09
-3.9266e-08 -3.3593e-08 -2.7422e-08
-1.3773e-07 -1.1790e-07 -9.6387e-08
-2.0609e-07 -1.7681e-07 -1.4509e-07
-1.7654e-07 -1.5278e-07 -1.2718e-07
-2.0138e-07 -1.7411e-07 -1.4477e-07
-2.1292e-07 -1.8320e-07 -1.5119e-07
-1.6936e-07 -1.4422e-07 -1.1704e-07
-7.9488e-08 -6.5040e-08 -4.9243e-08
-2.4655e-08 -1.6825e-08 -8.1479e-09
2.5263e-08 2.6401e-08 2.7786e-08
5.5296e-08 5.1247e-08 4.6820e-08
7.2724e-08 6.4964e-08 5.6363e-08
7.3627e-08 6.4585e-08 5.4582e-08
4.9704e-08 4.3216e-08 3.6187e-08
============================================================
It is possible that my ode15s function be changed or not be the same for you ( I have Matlab 7.7.0 (R2008b))?
Thank you for your time and for helping me with my problem!

Best regards,
A. Stoian



From: Torsten Hennig on 14 Jun 2010 07:12
Hi Alina,

Let's compare our results:

At which times did you output the concentrations listed
in you last mail, and to which h-values do they
belong ?

Best wishes
Torsten.
From: Alina Stoian on 14 Jun 2010 14:42
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <9593155.331588.1276528402857.JavaMail.root(a)gallium.mathforum.org>...
> Hi Alina,
>
> Let's compare our results:
>
> At which times did you output the concentrations listed
> in you last mail, and to which h-values do they
> belong ?
>
> Best wishes
> Torsten.

Hello,
I specified in tspan to output the results at each 3600 s ! So I displayed the values of the liquid concentration at each volume center (I remind you that I devised the liquid volume in 10 , and that the heights are all equal to 3e-4).
Best regards,
A. Stoian
From: Torsten Hennig on 14 Jun 2010 22:47
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <9593155.331588.1276528402857.JavaMail.root(a)gallium.ma
> thforum.org>...
> > Hi Alina,
> >
> > Let's compare our results:
> >
> > At which times did you output the concentrations
> listed
> > in you last mail, and to which h-values do they
> > belong ?
> >
> > Best wishes
> > Torsten.
>
> Hello,
> I specified in tspan to output the results at each
> 3600 s ! So I displayed the values of the liquid
> concentration at each volume center (I remind you
> that I devised the liquid volume in 10 , and that the
> heights are all equal to 3e-4).
> Best regards,
> A. Stoian

Hi Alina,

tightening the tolerances for rtol and atol to a very
small value (1e-11) keeps concentrations positive
with my solver.
To be sure about the results, you should make a second
calculation with a finer grid (maybe 5 times the number
of points in each compartment).

Best wishes
Torsten.
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4
Prev: date on x axis feather plot
Next: Need Help Attempted to access hx(0,0); index must be a positive integer or logical.