Modeling DAE system with discontinous solution (jump condition domain)

Hello,

I am trying build steady-state 1D electrodialysis model in matlab, cromprizing three sub domains, namely the analyte channel, the catholyte channel and an ion exchange membrane inbetween. I have provided the key eqations, while there are still some other minor algebraic equations.

DAE system:

(dc/dt=0)

Boundary condition: Neuman boundary conditions (also shown in the sketch)

Internal boundary condition (jump condition)

(membran solution interface continuity condition)

The masstransport is modeled using the nerst plank equation for each species and the electricfield is calculated through ohms law. Additionally I use the electron neutrality condition to close the mass balance. The individuall sub domains are connetced through an internal flux continuity boundary condition.

So far i have a rough understanding of the various ode solvers but from what i understand, the ode solvers only work when the whole problem is fully defined. In my case however, I have these sub domains connected via flux continuity boundary conditions, which is why i can not solve the sub domains individually. I have to solve them all at once. Additionally the boundary continuity boundary condition at the membrane interface introduces a jump (discontinuity) to the conentration and the potential profile. (qualitatively sketched below)

I am now looking for an appropirate way to model this system. The tricky part is the continuity condition at the membrane interface.

I hope to solve this problem with bvp5c where one can define multiple boundary condition inside the domain as described in this documentation:

https://de.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

But I am not entirlly sure how to implement the countinuity condition.

I have figured that it is possible to have jumps in the solution but since the continuity boundary condition which i am aiming to use is more complex that what has been done in the documentation, I am not quite sure on how to implent it. I would need the derivertives at the interface, instad of simply relating the absolute values at the interface to each other. (as it is done here)

function res = bc(YL,YR) % boundary conditions

res = [YL(1,1) % v(0) = 0

YR(1,1) - YL(1,2) - 1 % Here, I introduced a jump with a constant value at x=1

YR(2,1) - YL(2,2) % Continuity of C(x) at x=1

YR(2,end) - 1]; % C(lambda) = 1

end

Does anyone have an Idea, how i can use this solver for my problem. Any alternative approaches and clues are also very welcome.

Thank you^^

Bonus:

Ideally I want to build a hirachical model where all these three domains are modeled as individual systems (ode systems) and conntect them through a continiuty condtion for the concentraiton and the donnanpotential for the electric field. Also for the purpose of reusing the individuall models for later simulations. I know that this is common practice in gPROMS but i couldnt find similar approaches to modeling using matlab. Is there a similar modeling framework in matlab? I have only seen people using the various ode solvers for modeling in matlab for fully defined ODE systems.

[Torsten]

You want to solve for the steady-state profiles ?

[Berat Cagan Türkmen]

Yes, for the concentration of the individual species aswell as the potential.

[Berat Cagan Türkmen]

I have been trying to figure it out and i cant seem to find a solution to my problem, that is i do not see a way to implement the internal boundary conditions as continuities.

I now concider coding the mathematical method my self. However i have a hard time finding recources to code it my self. I am thinking about solving it using explicit finite difference method.

My main questions are as follows:

1. The balane equation contains the variables conentration c as well as potential phi. In all the recources i could find online there was only one variable, which were put in side a vector and are mulitpied with the matix to form the system of linear eqauations. but how do i handle 2 veriable in one equation?
2. additionally i have a system of ode eqations. how to I hadle that? Can i simply solve them indiviually by forming linear quations system for each ode equation?

But of course if there is a way to solve the problem with a suffisticated ode solver from matlab i would glad.

[Torsten]

I don't understand the background of your equations which I think would be necessary to help.

What variables are known, what are unknown ?

How are ϕ and φ related ?

Is c related to the C_i ?

[Berat Cagan Türkmen]

Hey thank for feedback,

you are right my intial question was strongly missleading because i screenshoted them form different recources. i rewrote the equations with some additional expainations. I hope that helps.

to still briefly answer your questions:

1. ϕ and φ are the same
2. c and C_i are the same

Answer 1

[Berat Cagan Türkmen]

I will try to explain the phyiscal backgroung of the equations used. Aigan, i want model the masstransport of ions in a chemical cell.

I start with the general continuity equation and since i am only interested in in the steady state solution i set

Flux of of ions are describe with the so called nernst plank equation.

Here, the frist term descirbes the ion flux due to diffusion (conenctraion gradient is the driving force) and the second term is the flux due to the potential gradient (electric field is the driving force)

by combing the flux with the continuity equation we get the first equation.

these are the balance equation for the ionic species in our system. when whe have N diffierence species we need N-1 of these balance equations for each region. Remember we had three different regions anaolyte, catholyte and membrane.

The last species is calculated via the electroneutrality condition, which basically assumes that at each point in our system the potential must cancel out.

to solve the system we still need an equation to describe the potential field (electric field) inside the system. We do this using the relationship between the current density and flux of ions. (we basicly say that that the flowing current i our system is basicly the flow of the electroylte...i hope this is intuitive)

here we can insert the ion flux given by the nerst plank equation and isolate , wich gives is the following equation

(i made a small adjustment to this equation by removing the third term in the numerator)

This should be the whole set of equaitons where every variable is known exept for the species concentreation and the potential ϕ, which i want to solve for.

As for the boundary conditions:

I have neuman boundary condition. L is basicly the maximal lenghth of the system. i compute form 0 to L. is known.

As sketched above i have two membrane/solution interfaces where the concentration jumps. At these interfaces i want to implement the continuity condition. is the ion flux intering the interface form the solution side and is the ion flux leaving form the membrane site (and vice verca)

Aigan i am inserting the ion flux from the nerst plank equation which yield the final equation.

I need to someone incorporate this equation as an internal boundary condition, which i am currently failing at.

[Torsten]

By inserting

dphi/dx = (I/F + sum_j [z_j*D_j/(R*T)]/(-F/(R*T) * sum_j [z_j^2*D_j*c_j])

into the equation

0 = -grad(D_i*dc_i/dx + z_i*D_i/(R*T)*c_i*dphi/dx) (i=1,...,N-1)

you get (N-1) second-order differential equations only in the c_i - thus phi is eliminated.

Also the two necessary transmission conditions at the interfaces can be written only in terms of the c_i and the dc_i/dx.

One is the relation between the fluxes, another will concern the concentrations on both sides, I guess.

So it should be no problem to use bvp4c or bvp5c as a boundary value problem with internal boundary conditions, shouldn't it ?

But this is only an ad hoc suggestion for a solution. There should be advices in the literature on how this well-studied problem can be solved stably and efficiently.

[Berat Cagan Türkmen]

Your first point is very helpful I will deffinetly do that.

But as for the the interface condition I don't know how to properly code it.

From the documentation you can see, that you deffinetly use the specific values at the boundary in the interface condition, by using YL and YL.

But in my case u also need the derivative of c_i at the boundary as the condition contains the dc/dx. And that is where I am failing.

Do you by chance know wether that us even possible?

[Torsten]

You have a second-order ODE, so you will have to write one equation for c_i and a second one for dc_i/dx in order to use bvp4c/bvp5c.

This is done as

dy(1) = y(2)

dy(2) = your equation for c_i solved for d^2c_i/dx^2.

Thus YL(1) = c_i from left, YR(1) = c_i from right, YL(2) = dc_i/dx from left, YR(2) = dc_i/dx from right.

[Berat Cagan Türkmen]

Okey that makes perfect sense i will try that! Thanks alot!

[Berat Cagan Türkmen]

I also want to follow this apporach a little further, becasue i found the physically correct interface boundary condition.

the new system is as follows with N being the numer of species:

Boundary conditions:

internal boundary condition (donnan potential)

I wrote down the whole set of equations for sake of completeness but only the internal boundry condittion chagned.

I think also for this set of equations the bvp5c solver reaches its limits becasue i can not think of a way to give all the internal boundary condtion with only this internal bodniary condition. If there is please enlighten me.^^

thefore im looking to implement a solver my self. is there a mathematical method i can use to also include the internal condition? the recources i could find on finite difference etc. where only for basic single ode and not for a DAE system as i have it here.

Thank you!!!

[Torsten]

You try to solve N first-order equations. Thus at each internal boundary, bvp4c/bvp5c expects N transmission conditions. Further you have to set a total of N conditions at the outer boundaries. Thus you have to derive a total of 3*N transmission/boundary conditions.

This is a mathematical fact - you won't be able to circumvent this if you start writing your own code. And bvp4c/bvp5c will try to handle these conditions if they make sense.

But in order to be able to solve for phi and not only dphi/dx, you will have to fix phi somewhere. Do you intend to set phi at one or both outer boundaries ?

[Berat Cagan Türkmen]

Thanks for the insight. I though coding the solver my self is necceray because bvp5c only acepts a Set of first order odes. But I also have the algebraic equation.

I want to fix Phi at both total boundaries.

[Torsten]

Why do you have an algebraic equation ? You have a differential equation for phi:

dphi/dx = ....

[Berat Cagan Türkmen]

I still have the electroneutrality condition:

to compute c_i of the Nth species

[Torsten]

But you can solve this equation for c_N - you don't need a solver for this ?!

c_N = -1/z_N * sum_{i=1}^{i=N-1} z_i*c_i

[Berat Cagan Türkmen]

I acutally need the concentration of c_N because it is part of the dphi/dx equation right?

or do you mean i should just insert into the your equation into the ode maually?

[Torsten]

I acutally need the concentration of c_N because it is part of the dphi/dx equation right?

or do you mean i should just insert into the your equation into the ode maually?

I don't understand what you mean. If you need c_N and you solve for c_1,...,c_(N-1), you can deduce c_N using the above formula. Same for dc_N/dx.

[Berat Cagan Türkmen]

For some reason i though, bvp5c can not handle algebraic equaitons and all the functions i provide have to be frist order odes. But that is acutally not the case. thank you for the hint.

Now here is my current progress. i can provide 4 boundary condtion at the outer boundaries in total:

And on the internal boundaries i can provide 3 for each species:

You try to solve N first-order equations. Thus at each internal boundary, bvp4c/bvp5c expects N transmission conditions. Further you have to set a total of N conditions at the outer boundaries. Thus you have to derive a total of 3*N transmission/boundary conditions.

I tried to reduce the system into first order ode's but i get 4 equations. That means i get N+1 instead of N. Did i do something wrong?

This way i have to provied 12 conditions in total. 4 at the outer boudnaries and 4+4 at the internal boundaries.

is there a way of formulating the set of equations into 3 odes?

I tired for example to derive dphi/dx manually to get dphi^2/dx^2 and insert it into the second equation and solve it for d^2c/dx^2. but it gets impossibe.

function dydx = fun(x,y,region,p) % equations being solved

DH = p(1);

DHM = p(2);

F = p(3);

I = p(4);

R = p(5);

T = p(6);

f = p(7);

zH = p(8);

RH = p(9);

PHI = p(10);

PHIM = p(11);

dydx = zeros(4,1);

switch region

case 1 % x in [0 1]

dydt(1)=y(2);

dydt(2)=-zH*f*y(2)*y(3)-zH*f*y(1)*y(4);

dydt(3)=((I/F)+zH*DH*y(2))/(-f*zH^2*DH*y(1));

dydt(4)=y(3);

case 2 % x in [1 lambda]

dydt(1)=y(2);

dydt(2)=-zH*f*y(2)*y(3)-zH*f*y(1)*y(4);

dydt(3)=((I/F)+zH*DH*y(2))/(-f*zH^2*DH*y(1));

dydt(4)=y(3);

case 3 % x in [1 lambda]

dydt(1)=y(2);

dydt(2)=-zH*f*y(2)*y(3)-zH*f*y(1)*y(4);

dydt(3)=((I/F)+zH*DH*y(2))/(-f*zH^2*DH*y(1));

dydt(4)=y(3);

end

end

[Torsten]

For some reason i though, bvp5c can not handle algebraic equaitons and all the functions i provide have to be frist order odes. But that is acutally not the case. thank you for the hint.

bvp5c cannot solve algebraic equations, but as long as you can derive their solutions "manually", your system can also contain algebraic equations.

I tried to reduce the system into first order ode's but i get 4 equations. That means i get N+1 instead of N. Did i do something wrong?

I don't know why you need 4 equations. You need one for c_i, one for dc_i/dx and one for phi. To get d^2phi/dx^2 which has to be inserted in the first equation, differentiate the expression for dphi/dx for x.

As discussed before, in principle you only need 2 equations because dphi/dx and d^2phi/dx^2 can be expressed in terms of c_i, dc_i/dx and d^2c_i/dx^2 alone. But as you noted, you want to set conditions for phi at the outer boundaries.

[Berat Cagan Türkmen]

I don't know why you need 4 equations. You need one for c_i, one for dc_i/dx and one for phi. To get d^2phi/dx^2 which has to be inserted in the first equation, differentiate the expression for dphi/dx for x.

So that was what i tried, but i get this expression for dphi^2/dx^2:

after inserting it into my my first equaiton it looks like this:

please correct me if i am worng but this seems to be impossible to solve for dc_i^2/dx^2, since it is in the sum. I can solve for a specific species N but the remaining dc_i^2/dx^2's will still stay in the equation. that way i would not have first oder ode right?

[Torsten]

The index in the sums must be j, not i.

As far as I can see, you get (N-1) such equations which are linear in the unknows d^2c_i/dx^2 (i=1,...,N-1).

So a linear system of equations in the unknowns d^2c_1/dx^2, d^2c_2/dx^2,...,d^2c_(N-1)/dx^2. Easily solved using backslash.

But as I said previously: a literature research about the usual solution methods for this kind of physical system should be made in order to avoid meanders. I can only tell you what is possible - if it is efficient and constructive: I don't know.

What is the COMSOL system being solved ?

[Berat Cagan Türkmen]

lets us form some concrete equations to not get confused: lets assume i have 2 components c_1 and c_2. for these two componentes my quations system should look like this.

theoretically i would need to solve the first equation for right?

but even if i did so. i still have which bvp5c wont accept because it would not be an first order ode.

The index in the sums must be j, not i.

I can understand why this should be the case the dc^2/dx i am inserting should be the same as the dc^2/dx which is allready in the equation.

So a linear system of equations in the unknowns d^2c_1/dx^2, d^2c_2/dx^2,...,d^2c_(N-1)/dx^2. Easily solved using backslash.

i can not really follow.

EDIT: and as you can see form the equation above. the dphi/dx equation also has to include the Nth species even if we dont spcificlally formulate a differential eqation for it.

But as I said previously: a literature research about the usual solution methods for this kind of physical system should be made in order to avoid meanders. I can only tell you what is possible - if it is efficient and constructive: I don't know.

in litterature i could find people solving almost the same system with pdepe but they did not have internal boudnaries as i have them. they only calulated the flow over the membrane. but i also want to include the solution phase. but yeah, other than that the differential eqation system is the same.

What is the COMSOL system being solved ?

i dont know if i understand you right. but the comsol system is more or less the exact same as i am trying to do. for the outer boudnary conditions i want to use neumann boudaries, while they used dirichlet. but the set of eqations we are using are the same.

[Torsten]

but even if i did so. i still have which bvp5c wont accept because it would not be an first order ode.

d^2c_2/dx^2 = -z1/z2 * d^2c_1/dx^2

by the relation

z1*c_1 + z2*c_2 = 0

And then you solve the equation for d^2c_1/dx^2, substitute y(1) = c_1, y(2) = dc_1/dx as usual. What's the problem ?

The index in the sums must be j, not i.

I can understand why this should be the case the dc^2/dx i am inserting should be the same as the dc^2/dx which is allready in the equation.

i is the index of the equation being solved (namely for c_i), j is the summation index over all c_j coming from substituting the expression for dphi/dx.

in litterature i could find people solving almost the same system with pdepe but they did not have internal boudnaries as i have them. they only calulated the flow over the membrane. but i also want to include the solution phase. but yeah, other than that the differential eqation system is the same.

At the moment, we do not discuss about the internal boundaries, but about how to set up the system to be solved (which solution variables are to be used; should there be an extra equation for phi or should phi be completely substituted by the c_i). I think if you find the way of solving the "usual" equations in the literature (without extra difficulties) would already be of great help.

[Berat Cagan Türkmen]

d^2c_2/dx^2 = -z1/z2 * d^2c_1/dx^2

by the relation

z1*c_1 + z2*c_2 = 0

And then you solve the equation for d^2c_1/dx^2, substitute y(1) = c_1, y(2) = dc_1/dx as usual. What's the problem ?

this works for 2 species, you are right. but i have more then 2 species i cant do that.

i is the index of the equation being solved (namely for c_i), j is the summation index over all c_j coming from substituting the expression for dphi/dx.

okey, i believe i understand what you mean but in the end i still have all the dc_1^2/dx^2,dc_2^2/dx^2,dc_3^2/dx^2,dc_4^2/dx^2 etc. which prevent me from solving the eqation for each individual dc_i^2/dx^2. the only difference is that the i is for N-1 species and and j is for all the N species

At the moment, we do not discuss about the internal boundaries, but about how to set up the system to be solved (which solution variables are to be used; should there be an extra equation for phi or should phi be completely substituted by the c_i). I think if you find the way of solving the "usual" equations in the literature (without extra difficulties) would already be of great help.

okey, I will look for information on that regard.

The problem i see however, is that the bvp5c model forces me to form first order odes. but if i could code a finite difference method based solver, shouldnt i be able to manually discritize me second order and the dphi/dx equation?

[Torsten]

d^2c_2/dx^2 = -z1/z2 * d^2c_1/dx^2

by the relation

z1*c_1 + z2*c_2 = 0

And then you solve the equation for d^2c_1/dx^2, substitute y(1) = c_1, y(2) = dc_1/dx as usual. What's the problem ?

this works for 2 species, you are right. but i have more then 2 species i cant do that.

Write down the system for N species. If you arrange the equations with regard to d^2c_i/dx^2 (i=1,...,N-1) and substitute d^2c_N/dx^2 as you substituted d^2c_2/dx^2 above, you get a linear system of equations of the form

a11*d^2c_1/dx^2 + a12*d^2c_2/dx^2 + ... + a1,N-1*d^2c_(N-1)/dx^2 = b_1
a21*d^2c_1/dx^2 + a22*d^2c_2/dx^2 + ... + a2,N-1*d^2c_(N-1)/dx^2 = b_2
    ...
aN-1,1*d^2c_1/dx^2 + aN-1,2*d^2c_2/dx^2 + ... + aN-1,N-1*d^2c_(N-1)/dx^2 = b_N-1

Here, a_ij and bi depend on c_i and dc_i/dx.

This system can be solved for the d^2c_i/dx^2 using backslash, as I already told you in a previous response.

The problem i see however, is that the bvp5c model forces me to form first order odes. but if i could code a finite difference method based solver, shouldnt i be able to manually discritize me second order and the dphi/dx equation?

Where do you see a problem in writing second-order derivatives as two first-order derivatives ?

In my experience, self-written numerical code (e.g. to solve a multipoint boundary value problem) should only come into play if the available solvers are not applicable. I don't see that this is the case here.