How to solve a system of pdes coupled in their boundary condition

Initially i had the following system of pdes to solve:

with the following initial and boundary conditions

at

at and

Both domains are coupled by a first order reversible reaction such as

at

where are constants.

I have successfuly solved the following system using an implicit euler scheme for the x-space derivative and second order central difference for the y-space second derivative for the internal nodes. Points at and were discretized using second order backward discritization.

Then to link between both domains the interface conditions was discritized also using backward discritization.

Now i want to developp more my model by considering the following system of 4 pdes:

the boundary conditions are the same as the previous system except that this time the interface condition is a non linear condition of the following form

at

and now i am confused on how to solve numerically such a system i have read about the segregated approach to solve coupled pdes. My understanding is that we fix all dependant variable but one where we solve for it and do the same for each variable. Then we iterate until convergence. However, i am not really sure how to implement it.

If anyone can point me in the right direction or useful literature, i would be grateful. Note that if it is feasable, i want to stick with my previously developed scheme and then elaborate it to take into account the coupling.

Thanks in advance for any help

[Torsten]

Any condition like

A = A1
B = B1 

at the coupling point ?

Or is there a jump in the profile ?

[hamza karim]

Hello Torsten,

Thanks for your reply, there is no such condition at the coupling plane. The pdes are only related by the flux conditions as mentioned above.

however to satisfy the system of pdes there is also no flux boundary conditions () at y=0 for A and A1 and at y=1 for B and B1

Answer 1

[Torsten]

Ok, then the discretization at y = alpha should be for the nonlinear case

(A(n+1)-A(n-1))/(2*dy)=-Bi_A * A(n)*A1(n)^2*B(1)
(A1(n+1)-A1(n-1))/(2*dy) = -2*Bi_A1 * A(n)*A1(n)^2*B(1)
(B(2)-B(0))/(2*dy) = 2*Bi_B * A(n)*A1(n)^2*B(1)
(B1(2)-B1(0))/(2*dy) = -Bi_B1 * A(n)*A1(n)^2*B(1)

These are four equations for A and A1 in the artificial point yA(n+1) and for B and B1 in the artificial point yB(0).

Assuming that your PDE equations also hold at the interface point, you can use the usual central difference discretization for the 2nd derivatives through the points yA(n-1), yA(n), yA(n+1) and yB(0),yB(1) and yB(2) to get A(n), A1(n), B(1) and B1(1).

[hamza karim]

sorry but i don't think i am following your approach.

in my discritization the domain is divided as follows. Supose we have Ny nodal points in the y direction.

at j=0 and j=Ny i discritize the no flux condition using second order backward discritization.

now for the interface, i have already divided it into two nodal points, at j=Ny/2 i express the variables in the first phase (0<y<h) and at j=Ny/2+1 i express those in the second phase (h<y<H).

but i must note that both j=Ny/2 and Ny/2+1 represents the same point h.

as far as the discritization. when it was a linear condition (the first system)

i obtain the following for the first condition

and

and i solve for both A and B by putting the in one matrix C

so in this case instead of your central difference i make my backward difference ? and how about the non linear constraint should there not be an iteration process until convergence ?

thanks again!

[Torsten]

so in this case instead of your central difference i make my backward difference ?

Yes, it seems so.

and how about the non linear constraint should there not be an iteration process until convergence ?

Sure. I always discretize in space and leave the time derivative in its continuous form. The resulting system of ordinary differential equations is then solved by ODE15S, e.g. . (Method-of-lines approach).

For this approach it is irrelevant if the Jacobian is constant or depends on the solution - so it doesn't matter whether the problem is linear or nonlinear. The ODE integrator does the job for me.

How about the values you get for A(ny/2) and B(ny/2) ? Is there really a jump in temperatures ? If yes: Isn't that unphysical ?

[hamza karim]

Hello torsten,

Thanks a lot ! i will try and implement the method of lines approach, especially that as you said the ode integrators take into account the non linearity in the system.

and yes there is a jump in conditions for A and B but they are not temperatures but concentrations and in mass transfer we are often exposed to such phenomena.

[Torsten]

I thought Biot-number refers to heat transfer, but maybe I'm mistaken.

[hamza karim]

no you haven't ! indeed Biot number refers to heat transfer but there is also an equivalent for the mass transfer