solving a system of COUPLED partial differential equations

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: Alina Stoian on 8 Jun 2010 06:59

---

Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1639689742.292475.1275920746356.JavaMail.root(a)gallium.mathforum.org>...
> > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> > in message
> > <735098493.260351.1275405399232.JavaMail.root(a)gallium.
> > mathforum.org>...
> > > > "Hello,
> > > >
> > > > I'm trying to solve a system of 3 partial
> > > > differential equations (3 diffusion equations,
> > 2nd
> > > > law of Fick) using ode15s! First of all, I
> > > > transformed my system in an ODE system using the
> > > > finite difference method to discretize each
> > equation.
> > > > After, I attached the boundary conditions , that
> > > > after the discretization are become algebraic
> > > > equations , and so in the end I' m obliged to
> > solve a
> > > > DAE system. From my readings I understood that I
> > have
> > > > to use ode15s to solve M*y'=f(t,y)! I defined the
> > > > matrix M, by putting 0 where ever I have an
> > algebraic
> > > > equation, and I 've created a separate function
> > file
> > > > for f(t,y), in which I defined the coefficients
> > of
> > > > the variables. The function f(t,y) has the form
> > > > B(n,n)*C(n), where B is the matrix with the
> > > > variable's coefficients and C is the variable
> > vector!
> > > > After the first run of the code I obtained
> > negative
> > > > results so I thought that it's best to
> > > > adimensionalize the variables, and I did so, and
> > now
> > > > I'm confronting with a warning message :'Warning:
> > > > Failure at t=3.124127e+00. Unable to meet
> > integration
> > > > tolerances without reducing the step
> > > > size below the smallest value allowed
> > (7.105427e-15)
> > > > at time t. '' . I know that the coefficient 's
> > matrix
> > > > is ill conditionned but what I don't know is how
> > to
> > > > handle this problem!""
> > > >
> > > >
> > > > I forgot to mention that the equations are
> > coupled
> > > > and that I'm trying to simulate the transfer of a
> > > > specie through three different phases ( that's
> > why I
> > > > have 3 basic diffusion equations).
> > >
> > >
> > > Did you try MATLAB's pdepe ?
> > >
> > > Best wishes
> > > Torsten.
> >
> >
> > Hello,
> > So I did what you told me and the answer form matlab
> > is :
> > ''aliq =
> >
> > 1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00
> > +00 1.0000e+00
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9999e-01
> > -01 9.9357e-01
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9996e-01
> > -01 9.8743e-01
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9990e-01
> > -01 9.8149e-01
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9983e-01
> > -01 9.7574e-01
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9974e-01
> > -01 9.7013e-01
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9963e-01
> > -01 9.6466e-01
> > 1.0000e+00 1.0000e+00 1.0000e+00 9.9950e-01
> > -01 9.5932e-01
> > 1.0000e+00 1.0000e+00 9.9999e-01 9.9936e-01
> > -01 9.5414e-01
> > 1.0000e+00 1.0000e+00 9.9999e-01 9.9920e-01
> > -01 9.4951e-01
> > == 1.0000e+00 1.0000e+00 9.9999e-01
> > 9.9903e-01 9.5065e-01=== 11th time step
> > 1.0000e+00 1.0000e+00 9.9999e-01 9.9894e-01
> > -01 1.0235e+00
> > 1.0000e+00 1.0000e+00 1.0000e+00 1.0166e+00
> > 00 1.4504e+01
> > 1.0000e+00 1.0000e+00 1.0003e+00 1.2295e+00
> > +00 1.7303e+02
> > 1.0000e+00 1.0000e+00 1.0039e+00 3.9296e+00
> > +00 2.1823e+03
> > 1.0000e+00 1.0001e+00 1.0499e+00 3.8166e+01
> > +01 2.7658e+04''
> > it seems that is working well until it reach the
> > 11th time step! Anyway I know that my J matrix is not
> > stable because the rcond is 1.475e-10, but from my
> > point of view my equations are written correctly.
> > I tried to use' pdepe ' but I have problems to impose
> > two conditions at the same boundary ( that are the
> > flow continuity and an equilibrium condition; the
> > flow continuity is like :
> > D_liq*dcliq/dx=D_gaz*dcgaz/dx (x= liquid height )
> > and cgaz=P*cliq, where P is a partition
> > n coefficient)!
> > Thank you very much for answering me and if you have
> > any other ideas please tell me!
> > Best regards,
> > A. Stoian
>
> If I got it right, you have interface conditions
> at the two phase boundaries. How did you discretize
> them ?
>
> Best wishes
> Torsten.
Hello,
First I have to tell you that I 'm trying to use the finite volumes method. That means that I divide each phase in a number of volumes , and I suppose that in each volume the concentration is homogeneous.

Yes , I have in fact 3 interfaces :, for the first one I discretize the equations like :
-liquid-gas interface :
Dliq*(Cliq,n - Cliq,i) / hliq,n / 2 = Dgas*(Cgas,i - Cgas,1) / hgas,1 / 2
Cgas,i = K * Cliq, i
where : K =partition coefficient
Dliq, Dgas = diffusion coefficient in of my compound in the liquid and in the gas
hliq,n/2 = the height from the center of the Nth liquid volume to the interface
hgas,1/2 = the height from the interface to the center of the first gas volume
Cliq,n = concentration in the Nth liquid volume
Cgas,1 = concentration in the first gas volume
So I obtain two equations !
....and so on for the others 2!
In fact now I am confronting with another issue, I get negative concentration when I run the non adimensionalized version and even for adimensionalized one I get that!
If it is necessary I can paste my code here !
Thank you !
Best regards,
Alina Stoian

From: Torsten Hennig on 8 Jun 2010 04:36

---

> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <1639689742.292475.1275920746356.JavaMail.root(a)gallium
> .mathforum.org>...
> > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de>
> wrote
> > > in message
> > >
> <735098493.260351.1275405399232.JavaMail.root(a)gallium.
> > > mathforum.org>...
> > > > > "Hello,
> > > > >
> > > > > I'm trying to solve a system of 3 partial
> > > > > differential equations (3 diffusion
> equations,
> > > 2nd
> > > > > law of Fick) using ode15s! First of all, I
> > > > > transformed my system in an ODE system using
> the
> > > > > finite difference method to discretize each
> > > equation.
> > > > > After, I attached the boundary conditions ,
> that
> > > > > after the discretization are become algebraic
> > > > > equations , and so in the end I' m obliged to
> > > solve a
> > > > > DAE system. From my readings I understood
> that I
> > > have
> > > > > to use ode15s to solve M*y'=f(t,y)! I defined
> the
> > > > > matrix M, by putting 0 where ever I have an
> > > algebraic
> > > > > equation, and I 've created a separate
> function
> > > file
> > > > > for f(t,y), in which I defined the
> coefficients
> > > of
> > > > > the variables. The function f(t,y) has the
> form
> > > > > B(n,n)*C(n), where B is the matrix with the
> > > > > variable's coefficients and C is the variable
> > > vector!
> > > > > After the first run of the code I obtained
> > > negative
> > > > > results so I thought that it's best to
> > > > > adimensionalize the variables, and I did so,
> and
> > > now
> > > > > I'm confronting with a warning message
> :'Warning:
> > > > > Failure at t=3.124127e+00. Unable to meet
> > > integration
> > > > > tolerances without reducing the step
> > > > > size below the smallest value allowed
> > > (7.105427e-15)
> > > > > at time t. '' . I know that the coefficient
> 's
> > > matrix
> > > > > is ill conditionned but what I don't know is
> how
> > > to
> > > > > handle this problem!""
> > > > >
> > > > >
> > > > > I forgot to mention that the equations are
> > > coupled
> > > > > and that I'm trying to simulate the transfer
> of a
> > > > > specie through three different phases (
> that's
> > > why I
> > > > > have 3 basic diffusion equations).
> > > >
> > > >
> > > > Did you try MATLAB's pdepe ?
> > > >
> > > > Best wishes
> > > > Torsten.
> > >
> > >
> > > Hello,
> > > So I did what you told me and the answer form
> matlab
> > > is :
> > > ''aliq =
> > >
> > > 1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00
> > > +00 1.0000e+00
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9999e-01
> > > -01 9.9357e-01
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9996e-01
> > > -01 9.8743e-01
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9990e-01
> > > -01 9.8149e-01
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9983e-01
> > > -01 9.7574e-01
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9974e-01
> > > -01 9.7013e-01
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9963e-01
> > > -01 9.6466e-01
> > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9950e-01
> > > -01 9.5932e-01
> > > 1.0000e+00 1.0000e+00 9.9999e-01 9.9936e-01
> > > -01 9.5414e-01
> > > 1.0000e+00 1.0000e+00 9.9999e-01 9.9920e-01
> > > -01 9.4951e-01
> > > == 1.0000e+00 1.0000e+00 9.9999e-01
> > > 9.9903e-01 9.5065e-01=== 11th time step
> > > 1.0000e+00 1.0000e+00 9.9999e-01 9.9894e-01
> > > -01 1.0235e+00
> > > 1.0000e+00 1.0000e+00 1.0000e+00 1.0166e+00
> > > 00 1.4504e+01
> > > 1.0000e+00 1.0000e+00 1.0003e+00 1.2295e+00
> > > +00 1.7303e+02
> > > 1.0000e+00 1.0000e+00 1.0039e+00 3.9296e+00
> > > +00 2.1823e+03
> > > 1.0000e+00 1.0001e+00 1.0499e+00 3.8166e+01
> > > +01 2.7658e+04''
> > > it seems that is working well until it reach
> the
> > > 11th time step! Anyway I know that my J matrix is
> not
> > > stable because the rcond is 1.475e-10, but from
> my
> > > point of view my equations are written correctly.
>
> > > I tried to use' pdepe ' but I have problems to
> impose
> > > two conditions at the same boundary ( that are
> the
> > > flow continuity and an equilibrium condition; the
> > > flow continuity is like :
> > > D_liq*dcliq/dx=D_gaz*dcgaz/dx (x= liquid height )
> > > and cgaz=P*cliq, where P is a partition
> > > n coefficient)!
> > > Thank you very much for answering me and if you
> have
> > > any other ideas please tell me!
> > > Best regards,
> > > A. Stoian
> >
> > If I got it right, you have interface conditions
> > at the two phase boundaries. How did you discretize
> > them ?
> >
> > Best wishes
> > Torsten.
> Hello,
> First I have to tell you that I 'm trying to use the
> finite volumes method. That means that I divide each
> phase in a number of volumes , and I suppose that in
> each volume the concentration is homogeneous.
>
> Yes , I have in fact 3 interfaces :, for the first
> one I discretize the equations like :
> -liquid-gas interface :
> Dliq*(Cliq,n - Cliq,i) / hliq,n / 2 = Dgas*(Cgas,i
> ,i - Cgas,1) / hgas,1 / 2
> Cgas,i = K * Cliq, i
> where : K =partition coefficient
> Dliq, Dgas = diffusion coefficient in
> coefficient in of my compound in the liquid and in
> the gas
> hliq,n/2 = the height from the center of
> center of the Nth liquid volume to the interface
> hgas,1/2 = the height from the interface to
> erface to the center of the first gas volume
> Cliq,n = concentration in the Nth liquid
> liquid volume
> Cgas,1 = concentration in the first gas volume
> So I obtain two equations !
> ...and so on for the others 2!
> In fact now I am confronting with another issue, I
> get negative concentration when I run the non
> adimensionalized version and even for
> adimensionalized one I get that!
> If it is necessary I can paste my code here !
> Thank you !
> Best regards,
> Alina Stoian
>

I think it will be the best if you write out the
mathematical description of your problem
(equations, boundary, interface and initial conditions)
and the corresponding MATLAB code.

Best wishes
Torsten.

From: Alina Stoian on 9 Jun 2010 05:54

---

Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1628642300.298169.1276000626806.JavaMail.root(a)gallium.mathforum.org>...
> > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> > in message
> > <1639689742.292475.1275920746356.JavaMail.root(a)gallium
> > .mathforum.org>...
> > > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de>
> > wrote
> > > > in message
> > > >
> > <735098493.260351.1275405399232.JavaMail.root(a)gallium.
> > > > mathforum.org>...
> > > > > > "Hello,
> > > > > >
> > > > > > I'm trying to solve a system of 3 partial
> > > > > > differential equations (3 diffusion
> > equations,
> > > > 2nd
> > > > > > law of Fick) using ode15s! First of all, I
> > > > > > transformed my system in an ODE system using
> > the
> > > > > > finite difference method to discretize each
> > > > equation.
> > > > > > After, I attached the boundary conditions ,
> > that
> > > > > > after the discretization are become algebraic
> > > > > > equations , and so in the end I' m obliged to
> > > > solve a
> > > > > > DAE system. From my readings I understood
> > that I
> > > > have
> > > > > > to use ode15s to solve M*y'=f(t,y)! I defined
> > the
> > > > > > matrix M, by putting 0 where ever I have an
> > > > algebraic
> > > > > > equation, and I 've created a separate
> > function
> > > > file
> > > > > > for f(t,y), in which I defined the
> > coefficients
> > > > of
> > > > > > the variables. The function f(t,y) has the
> > form
> > > > > > B(n,n)*C(n), where B is the matrix with the
> > > > > > variable's coefficients and C is the variable
> > > > vector!
> > > > > > After the first run of the code I obtained
> > > > negative
> > > > > > results so I thought that it's best to
> > > > > > adimensionalize the variables, and I did so,
> > and
> > > > now
> > > > > > I'm confronting with a warning message
> > :'Warning:
> > > > > > Failure at t=3.124127e+00. Unable to meet
> > > > integration
> > > > > > tolerances without reducing the step
> > > > > > size below the smallest value allowed
> > > > (7.105427e-15)
> > > > > > at time t. '' . I know that the coefficient
> > 's
> > > > matrix
> > > > > > is ill conditionned but what I don't know is
> > how
> > > > to
> > > > > > handle this problem!""
> > > > > >
> > > > > >
> > > > > > I forgot to mention that the equations are
> > > > coupled
> > > > > > and that I'm trying to simulate the transfer
> > of a
> > > > > > specie through three different phases (
> > that's
> > > > why I
> > > > > > have 3 basic diffusion equations).
> > > > >
> > > > >
> > > > > Did you try MATLAB's pdepe ?
> > > > >
> > > > > Best wishes
> > > > > Torsten.
> > > >
> > > >
> > > > Hello,
> > > > So I did what you told me and the answer form
> > matlab
> > > > is :
> > > > ''aliq =
> > > >
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00
> > > > +00 1.0000e+00
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9999e-01
> > > > -01 9.9357e-01
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9996e-01
> > > > -01 9.8743e-01
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9990e-01
> > > > -01 9.8149e-01
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9983e-01
> > > > -01 9.7574e-01
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9974e-01
> > > > -01 9.7013e-01
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9963e-01
> > > > -01 9.6466e-01
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 9.9950e-01
> > > > -01 9.5932e-01
> > > > 1.0000e+00 1.0000e+00 9.9999e-01 9.9936e-01
> > > > -01 9.5414e-01
> > > > 1.0000e+00 1.0000e+00 9.9999e-01 9.9920e-01
> > > > -01 9.4951e-01
> > > > == 1.0000e+00 1.0000e+00 9.9999e-01
> > > > 9.9903e-01 9.5065e-01=== 11th time step
> > > > 1.0000e+00 1.0000e+00 9.9999e-01 9.9894e-01
> > > > -01 1.0235e+00
> > > > 1.0000e+00 1.0000e+00 1.0000e+00 1.0166e+00
> > > > 00 1.4504e+01
> > > > 1.0000e+00 1.0000e+00 1.0003e+00 1.2295e+00
> > > > +00 1.7303e+02
> > > > 1.0000e+00 1.0000e+00 1.0039e+00 3.9296e+00
> > > > +00 2.1823e+03
> > > > 1.0000e+00 1.0001e+00 1.0499e+00 3.8166e+01
> > > > +01 2.7658e+04''
> > > > it seems that is working well until it reach
> > the
> > > > 11th time step! Anyway I know that my J matrix is
> > not
> > > > stable because the rcond is 1.475e-10, but from
> > my
> > > > point of view my equations are written correctly.
> >
> > > > I tried to use' pdepe ' but I have problems to
> > impose
> > > > two conditions at the same boundary ( that are
> > the
> > > > flow continuity and an equilibrium condition; the
> > > > flow continuity is like :
> > > > D_liq*dcliq/dx=D_gaz*dcgaz/dx (x= liquid height )
> > > > and cgaz=P*cliq, where P is a partition
> > > > n coefficient)!
> > > > Thank you very much for answering me and if you
> > have
> > > > any other ideas please tell me!
> > > > Best regards,
> > > > A. Stoian
> > >
> > > If I got it right, you have interface conditions
> > > at the two phase boundaries. How did you discretize
> > > them ?
> > >
> > > Best wishes
> > > Torsten.
> > Hello,
> > First I have to tell you that I 'm trying to use the
> > finite volumes method. That means that I divide each
> > phase in a number of volumes , and I suppose that in
> > each volume the concentration is homogeneous.
> >
> > Yes , I have in fact 3 interfaces :, for the first
> > one I discretize the equations like :
> > -liquid-gas interface :
> > Dliq*(Cliq,n - Cliq,i) / hliq,n / 2 = Dgas*(Cgas,i
> > ,i - Cgas,1) / hgas,1 / 2
> > Cgas,i = K * Cliq, i
> > where : K =partition coefficient
> > Dliq, Dgas = diffusion coefficient in
> > coefficient in of my compound in the liquid and in
> > the gas
> > hliq,n/2 = the height from the center of
> > center of the Nth liquid volume to the interface
> > hgas,1/2 = the height from the interface to
> > erface to the center of the first gas volume
> > Cliq,n = concentration in the Nth liquid
> > liquid volume
> > Cgas,1 = concentration in the first gas volume
> > So I obtain two equations !
> > ...and so on for the others 2!
> > In fact now I am confronting with another issue, I
> > get negative concentration when I run the non
> > adimensionalized version and even for
> > adimensionalized one I get that!
> > If it is necessary I can paste my code here !
> > Thank you !
> > Best regards,
> > Alina Stoian
> >
>
> I think it will be the best if you write out the
> mathematical description of your problem
> (equations, boundary, interface and initial conditions)
> and the corresponding MATLAB code.
>
> 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 =partition coefficient

- Gas phase - dcg/dt = Dg * d2cg/dx2
- initial condition : cg = 0
- boundary condition :
- left side - x = hliq :
Dl * dcl/dx =Dg* dcg/dx
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* dcm/dx
cg = cm
- right side - x = ht :
Smem*Dm * dcm/dx = S*Kext *(cmext &#8211; cair)
Cm = cmext ; cair = 0 ; Kext &#8211; mass transfer coefficient superior compartment side

The finite volumes method &#8211; for the equations:
-liquid
VL,1*dcL,1/dt = 0 &#8211; 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) &#8211; 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) &#8211; S*Dl*(cL,N-cL*)/(hL,i/2) ; i=1..N
but VL,1 = S*hL,1 and with this my system becomes :
dcL,1/dt = 0 &#8211; Dl*(cL,1-cL,2)/(hL,1*(hL,1+hL,2)/2)
dcL,i/dt = Dl*(cL,i-1-cL,i)/[hL,i *((hL,i-1+hL,i)/2) &#8211; 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) &#8211; Dl*(cL,N-cL*)/(hL,N/2))] ; i=1..N

gas
dcG,1/dt = Dg*(CG*-cG,1)/(hG,1*hG,1/2) &#8211; 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-1+hG,i)/2) &#8211; 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) &#8211; Dg*(cG,P-cgm*)/(hG,P/2))] ; i=1..P
membrane
dcM,1/dt = Dm*(Cgm*-cM,1)/(hM,1*hM,1/2) &#8211; 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-1+hM,i)/2) &#8211; 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) &#8211; 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 =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* 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 = S*Kext *(cmext &#8211; cair)
Cm = cmext ; cair = 0 ; Kext &#8211; 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 (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_air2= 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
% %-------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
% 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(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','JConstant','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 (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));
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 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 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)
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- 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: Torsten Hennig on 9 Jun 2010 22:27

---

Hi,

it will take a while to look through your equations.
I think I'll be back on Monday.

Best wishes
Torsten.

From: Torsten Hennig on 11 Jun 2010 08:05

---

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

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.