Solving 2 nonlinear elliptic pdes in 2 separate domains (2D), coupled by a function of solution u

I was trying to model voltage distribution in 2 parallel planar electrode by conductive media DC equations. The details of the system is described in the attached pdf file. After going through the pde toolbox user guide, I formulated the problem as following

*for domain 1. –div(sigma.grad(v1))=(OCV-v1-v2)R

for domain 2. 1. –div(sigma.grad(v1))=-(OCV-v1-v2)R*

where OCV and R are constants.

the boundary conditions are applied as described in the document

As the source q is function of solution v1 and v2, I used non linear solver. And used a coefficient function for q by calling a matlab function. But as it 2 different pde in 2 region I could not use the solution from 2 pde in one function. Tried different syntax by using domain variable, for example u(sd==1) for solution v1.

As the q coefficient functions are evaluated one by one for each equations, after solving 1 equation for 1 domain, all the variables are cleared. So I could not use the solution for coupled equations (value of v1 and v2). Tried to define a global variable in the function to store instantaneous solution, but that variable was also cleared. In this regard, can you please support me by suggesting proper method to achieve this simulation setup?

Thank you in advance. BR, Shovon

Answer 1

[Alan Weiss]

I am not sure that I completely understand the problem, but perhaps you could try to do the following. If the two electrodes have identical geometry, then perhaps you could model the two solutions v1 and v2 as being two components of a solution u on a single unified geometry. Instead of having the two components separated by a dimension z, as in your diagram, maybe you could just write the components as being in exactly the same place.

The reason I think this might not work is that I do not understand the physics of the flow from one electrode to the other, and perhaps you have some equation that describes this flow (though I did not understand it), and you cannot regard the two electrodes as being right next to each other. But, if you can do this modeling, then perhaps you could get the nonlinear solver to obtain the two solutions u = [v1,v2] simultaneously.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Shovon Goutam]

Hello Alan, Thank you for your answer. The problem as I stated might be little confusing because the tranverse current on Z direction is a simplification of complex elctrochemical processes in a battery. I can not use one equation for both geometry as they are of different material and thus the conductivity is different. Although their geometry is the same is x,y dimension. So I had to use two different equation. But as it is a simplified model, the term J (transvers current in Z direction) do not depend on the distance between the plate. Therefore, I made the geometry on the x,y plane side by side.

Good news is after trying several thing, I finally figured it out. The problem was, I was using u(small letter) and sd(small letter) in the input argument of the coefficient function q. u and sd contains the corresponding value of 1 domain at a time. So as long as the equation is evaluated on domain 1 it is cleared to store the values of domain 2.

What I did next, I put a breakpoint in the function and went in the core functions of the pdetool step by step. There I found out, if I use U or SD (capital letter), those matrix contains the value of all the domain. Voila! So I wrote the following function to achive my goal.

function f = ffunc(sd,SD,U)

    % here, sd= domain matrix of individual domain
    %------ SD= domain matrix of whole geomtery
    %--------U= solution U over 2 domain
    %first separate the solution according to domian

    u1=U(SD==1);            %solution of domain 1
    u2=U(SD==2);            %solution of domain 2

    %calculate the f coffcient matrx according to domain

    f1=-(3.69-u1-u2)/(0.1354*2.1e-05);      %domain 1
    f2=-(3.69-u1-u2)/(0.1354*1.2e-5);       %domain 2

    %as the pde tool solve pdes 1 after other
    %f fucntion should be according to the working pde or domain

    if sd(1,1)==1
      f=f1;
    elseif sd(1,1)==2
      f=f2;
    end
end

end % code end

Unfortunately, this information regarding u and U was not mentioned explicitly in pde toolbox user guide, to my knowledge. Or may be I missed,may be you can confirm.

Although every thing is fine from Syntax point of view. The solver still can not solve it. Because of initial value. The initial value of u1 and u2 have a huge difference. (ca. 3 for ~u1 nad 2e-3 for u2). Thus I am getting error 'step size tool small'. I guess this is because I gave the initial value as 3, which create a matrix of u0 with each value as 3.

I am thinking of putting a similar function as the coefficient function q. I hope intial value can also take a function. I have not tried it yet. That is my next step. If i am successful I have to simulate the whole process for different value of OCV and R. I do not know yet how will I achieve that. I can not use workspace variable in the pdetool. Or can not argument on pdetool function (I guess). Can you please give me head start on that? Should I make the model by using step by step pde solving (recommended workflow by programming).

N.B.- In case, if you are interested, I can post update here, so others can get help as well.

Thank you again. Regards, Shovon.