How to solve dependent PDEs

I'm trying to solve the above set of PDEs for h and p which are both functions of r and t with the given initial and boundary conditions. I'm currently solving for F_f using gaussian quadrature for each t step. V is also a function of t that is solved using a simple first order ODE not shown here. After finding these two, I need to use them in the boundary conditions to solve the two dependent PDEs. I'm trying to use pdepe but am not sure if it's possible with this set of equations and conditions. Any advice would be very appreciated, thanks!

Answer 1

[Nicolas Schmit]

https://www.mathworks.com/help/matlab/math/partial-differential-equations.html

[Travis]

I've tried using pdepe but it tells me it is unable to meet integration tolerances. This is currently how I have the various pdepe functions set up. I'm not sure what exactly is causing the error.

Thanks

 function [c,f,s] = pde(r,t,u,DuDr)
    global mu1 sigmabar Rbar
    c = [3*mu1; 0]; 
    f = [u(1)^3*DuDr(2); DuDr(1)];
    s = [0; -4/Rbar+2*u(2)/sigmabar]; 
end
% --------------------------------------------------------------
function u0 = pdeic(r)
    global H00 Rbar
  u0=[H00+r^2/Rbar; 0];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
    global Ff alpha rm Rbar
    pl = [0; 0]; 
    ql = [ul(1)^-3; 1]; 
    pr = [ur(1)-alpha*Ff+ul(1)+rm^2/Rbar; ur(2)]; 
    qr = [0; 0]; 
end

[Torsten]

My advice is to discretize both PDEs in r and solve the resulting system of ordinary differential equations in t using ODE15S (Method-of-Lines - approach).

This way, you have maximum flexibility for the implementation of your boundary conditions (especially for h at r=r_m) and maximum debugging options.

As you write correctly, using a tool like pdepe, one is often usure about what might cause the problem.

Best wishes

Torsten.

[Travis]

Hi Torsten, thanks for the advice. This seems like a much more feasible method. To solve this I'm treating the above as a DAE since dp/dt is not explicitly given. This is then giving me a system of N DAEs where N is twice the length of my r domain. I'm attempting to solve that system of DAEs with ode15s but am getting the error that the DAE is of index greater than 1. Is there any more insight you might be able to provide?

Thanks,

Travis

Answer 2

[Torsten]

Interesting problem.

I'd try this way:

Multiplying the second equation (p = ...) by 2*pi*r and integrating both sides from r=0 to r=r_m gives

integral_{r=0}^{r=r_m} 2*pi*r*p dr = 2*sigma/R * pi*r_m^2 - sigma*pi*r_m*(dh/dr)|r=r_m.

Together with your boundary condition for h at r=r_m you get

(h(r_m)+h(0)+r_m^2/R)/alpha = 2*sigma/R * pi*r_m^2 - sigma*pi*r_m*(dh/dr)|r=r_m

thus

-r_m*dh/dr|r=r_m = ((h(r_m)+h(0)+r_m^2/R)/(alpha*sigma*pi) -2*r_m^2/R  (*)

Now given that h is known in

r1=0, r2=deltar, r3=2*deltar,...,rn=r_m,

you can calculate p for

r1'=deltar/2, r2'=3/2*deltar, ...,r_(n-1)'=r_m-deltar/2

via

p(ri') = 2*sigma/R - 1/2*sigma/ri' * (r_(i+1)*(dh/dr)|r=r_(i+1) - r_i*(dh/dr)|r=r_i)/deltar 

(dh/dr)|r=r_i = (h(r_(i+1))-h(r_(i-1)))/(2*deltar)
(dh/dr)|r=0 = 0
(dh/dr)|r=r_m = (from above expression (*))

Now for inner points

(dh/dt)|r=ri = 1/(3*mu*ri)*(r_(i+1)*h_(i+1)^3*(dp/dr)|r=r_(i+1)-r_(i-1)*h_(i-1)^3*(dp/dr)|r=r_(i-1))/(2*deltar)

with

(dp/dr)|r=ri = (p(ri')-p(r_(i-1)'))/deltar

For r=r_m:

(dh/dt)|r=r_m = 1/(3*mu*r_m)*(r_m*h_m^3*(dp/dr)|r=r_m - r_(m-1)*h_(m-1)^3*(dp/dr)|r=r_(m-1))/deltar

and

(dp/dr)|r=r_m = (p(r_m-deltar/2)-p(rm))/(deltar/2) = 2*p(r_m-deltar/2)/deltar

For r=0:

h(0) is explicitly given by

h(0)=-h(rm)+alpha*F_f - r_m^2/R.

So you only need to discretize the first PDE for h and solve it using ODE15S.

Given h on the grid r1=0, r2=deltar,...,rn=r_m, you can deduce p on the grid r1'=deltar/2, r2'=3*deltar/2,...,r(n-1)=r_m-deltar/2 and calculate dh/dt at r=ri. Integrating (dh/dt)|r=r_i using ODE15S gives h.

Best wishes

Torsten