How to solve dependent PDEs

5 次查看（过去 30 天）
Travis
Travis 2017-11-9
回答： Torsten 2017-11-14
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!

请先登录，再回答此问题。

回答（2 个）

Nicolas Schmit
Nicolas Schmit 2017-11-10
  3 个评论
显示 1更早的评论
Torsten
Torsten 2017-11-10
编辑：Torsten 2017-11-10
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
Travis 2017-11-11
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

请先登录，再进行评论。

Torsten
Torsten 2017-11-14
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

类别

在 Help Center 和 File Exchange 中查找有关 Electromagnetics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by