Using pdepe with different material properties

I want to solve a system of 7 PDE (second order mass balances) and 1 ODE (second order electric potential balance) within different materials .I'm writing a vector of ones for c with one value equals to zero (ODE), if statements for f and s and creating mesh points at interfaces, but I have a few questions on how to specify subdomains conditions:

1. How do I specify equal flux and a jump in concentrations at internal interfaces? (eg. D1*(du(1)/dx)=D2(du(1)/dx) & u(1)=K*u(1) where K is constant).

2. How should I specify that the potential balance (ODE) applies only at 1 subdomain? (eg. u(8)=0 everywhere except in subdomain 3 with boundary conditions du(8)/dx=0 and (u(8)+E(u(2),u(5))) -B du(8)/dx =0 with B constant).

Equations

c*du(i)/dt = d/dx(f(x)*du(i)/dx) + s(x,u)  ;  i=1,2..8
c=[1;1;1;1;1;0;0;0]
r1=k1*u(1)*(1/(1+exp(-F*u(8)/R/T)))/(Ks1+u(1))
r2=k2*u(1)/(Ks2+u(1))
r3=k3*u(5)*exp(beta*F*na/R/T)/(Ko+u(5))

On subdomain 1 and 4 c(6)=1; f=[Ds1;Dh1;Dco1;Dch1;Do1;0;0;0] s=[-u(6)*r2;0;u(6)*r2;u(6)*r2;0;Y2*u(6)*r2 - Kinam*u(6);u(7);u(8)] u0=[Ss0;Hs0;COs0;CHs0;Os0;Xm0;0;0]

On subdomain 2

f=[Ds2;Dh2;Dco2;Dch2;Do2;0?;0?;0?]
s=[0;0;0;0;0;u(6);u(7);u(8)]
u0=[Se0;He0;COe0;CHe0;Oe0;0?;0?;0?]

On subdomain 6

f=[Ds2;Dh2;Dco2;Dch2;Do2;0;0;0]
s=[0;-4*r3;0;0;-r3;u(6);u(7);u(8)]
u0=[Se0;He0;COe0;CHe0;Oe0;0?;0?;0?]

On subdomain 3 c(7)=1; f=[Ds3;Dh3;Dco3;Dch3;Do3;0?;0?;kbio] s=[-u(7)*r1;8*u(7)*r1;u(7)*r1;0;0;u(6);Y1*u(7)*r1-Kinae*u(7);-F*gamma*fe*r1] u0=[Sb0;Hb0;COb0;CHb0;Ob0;0;Xe0;0]

On subdomain 5 and 7

f=[Ds4;Dh4;Dco4;Dch4;Do4;0?;0?;0?]
s=[0;0;0;0;0;u(6);u(7);u(8)]
u0=[Sw0;Hw0;COw0;CHw0;Ow0;0;0;0]

The boundary conditions between subdomains are:

for i=1:5
Continuous flux on internal interfaces
u(i)_subdomain1,2,3,4,5,6 = K(x)*u(i)_subdomain2,3,4,5,6,7
for i=6,7,8
du(6)/dx = 0 on interfaces between subdomains 1-2 , 3-4, 4-5 
du(7)/dx = 0 on interfaces between subdomains 2-3 , 3-4
du(8)/dx = 0 on interface between subdomains  3-4
(u(8)+E{u(2),u(5)}) - B*kbio*du(8)/dx = 0 on interface between subdomains 2-3

So variable u(6) only takes values at subdomains 1 and 4 and is equal to 0 elsewhere. Variable u(7) and u(8) only take values at subdomain 3 and are equal to 0 elsewhere.

The missing boundary conditions are external (left and right ones), for left side every flux is zero and for right size every flux is zero, excepting u(3)=0, u(4)=0, u(5)-G=0.

Equation 8 is only function of x and applies on all points of subdomain 3. Because u(6), u(7) and u(8) are non-zero on specific subdomains, i'm not sure how to write c,f,s,u0 on the subdomains they are zero.

[Bill Greene]

Handling different materials with pdepe is straightforward. The documentation discusses this specifically: "Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface." You simply check the value of x in your pdefun and return the appropriate coefficients for that x.

However, if you really have an ODE, i.e. an equation that applies only at a single point in the domain, there is no straightforward way to handle that with pdepe. I suggest you modify your post to include the equations.

[Juan Diego Carvajalino Olave]

Ok I edited my post. Could be a better option to solve it with the pde toolbox?.

[Bill Greene]

When I suggested adding the equations to your post, I really meant in mathematical form rather than trying to write pdepe code. It is hard to understand this system in the form you have written it.

Apparently, you misunderstood my comment about discontinuous material properties. pdepe generally makes the flux continuous between domains but does not permit any "internal" boundary conditions.

You did not specify equation 8 but did say u(8) is a function of x so must be satisfied everywhere in the domain. It is probably possible for pdepe to handle this.

Because your PDE system is 1D, pde toolbox offers no advantages over pdepe.

Do you have a document that describes this system of equations in more detail?

[Juan Diego Carvajalino Olave]

I attached the mathematical form and a better description. The model is mainly based on the work found here: link.springer.com/content/pdf/10.1007%2Fs40095-014-0093-5.pdf

[Bill Greene]

I don't know how to handle a jump condition like u(1)=K*u(1) in pdepe. For this problem, what does it mean to have a sharp discontinuity of the solution?

[Juan Diego Carvajalino Olave]

Means that the mass concentration of a substance in one phase (solid side interface) is different in the adjacent phase (liquid side interface), and the ratio between the two concentrations is constant. It is due to different solubility of a substance in each phase (see https://en.wikipedia.org/wiki/Distribution_constant). For a start, i could assume that those discontinuities do not exist and concentrations are equal across internal interfaces.

[Juan Diego Carvajalino Olave]

About the solution u(6), u(7) and u(8) having no flux across internal boundaries (and taking non-zero values in specific subdomains); For the subdomains where u(6) u(7) or u(8) are zero, should i specify the coefficients c f s to make the equation u(6)=0 in those subdomains for example? or Should I represent the no flux condition to those subdomains as very little diffusion coefficients?

[Bill Greene]

I don't really understand the difference in the two options you propose for those subdomains. It will probably be useful to specify a small amount of diffusion in those subdomains to smooth out the numerical effects of the sharp change in coefficients.

I suggest doing some numerical experiments on a simpler problem, particularly if this is your first experience with pdepe.

[Juan Diego Carvajalino Olave]

Thank you I will.

Answer 1

[Tao Chen]

We are solving similar diffusion problems, trying Matlab and COMSOL.

In COMSOL one could specify internal interfaces at boundaries and specify the flux; e.g. J = a * (C_1 - K C_2), where 'a' is a big number to ensure equilibrium at the interface. But it seems Matlab pdepe or PDE toolbox does not allow boundary conditions between subdomains.

You may try to reformat the diffusion equation, so that the variable being simulated is not concentration, but chemical potential. After all chemical potential is the actual driving force for diffusion, and it should be continuous across the interface. Then postprocessing can convert the chemical potential back into concentration. It works for us.

It's worth checking the conservation of mass in the results though for the finite element based solvers.