System of PDEs which is tricky for PDEPE

8 次查看（过去 30 天）
Matthew Hunt
Matthew Hunt 2018-5-23
回答： Matthew Hunt 2018-6-15
I have a system of PDEs, mainly diffusion equations of the form:
T_t-(k(x)(T_x)_x=a*E^2
c_t-(D(x)*c_x)_x=d*(E_x+c_x-T_x)
(epsilon*E)_x=-b*c
Where a and d are constants and _t,_x represent partial differentiation w.r.t. t and x respectively. In terms of pdepe, I would have: f=0 for the third equation. Would this cause errors in the code?
I'm also struggling to see how I can input the boundary conditions.
Any suggestions?
  5 个评论
显示 3更早的评论
Matthew Hunt
Matthew Hunt 2018-5-24
That is the boundary condition, you can ignore one of them.
Bill Greene
Bill Greene 2018-5-24
Do you mean T=0 at the boundaries? What about the initial conditions on the three variables? At t=0, the boundary and initial conditions must be consistent and also satisfy your differential equations.

请先登录，再进行评论。

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

回答（5 个）

Torsten
Torsten 2018-5-23
编辑：Torsten 2018-5-23
"pdepe" is designed to solve systems of parabolic-elliptic pdes. Your third pde is hyperbolic in nature. Thus "pdepe" is not suited to solve your system.
You will have to discretize your equations in space and solve the resulting system of ordinary differential equations using ODE15S.
Look up "method-of-lines" for more details.
Best wishes
Torsten.
  2 个评论
Matthew Hunt
Matthew Hunt 2018-5-23
Cheers,
I had a feeling that this might be the case. You will note however that it is possible to solve the hyperbolic equation completely thus giving E as a function of c essentially. Can I drop that equation somehow and use the integral of c as part of the code for is it too hard to do that?
Torsten
Torsten 2018-5-23
Since "pdefun" is called only for a single value of x, you don't have the complete vector for c available. Thus no chance to calculate E, I guess.

请先登录，再进行评论。

Torsten
Torsten 2018-5-24
You could try to use "pdepe" with the third equation differentiated:
(epsilon*E)_xx=-b*c_x
with Dirichlet boundary condition
E = E0
at one end of the interval and
(epsilon*E)_x + b*c=0
at the opposite end.
Best wishes
Torsten.
  3 个评论
显示 1更早的评论
Torsten
Torsten 2018-5-24
But the E_x and T_x terms in the diffusion equation for c will be put in the source term s for this equation, not in the f-term.
Matthew Hunt
Matthew Hunt 2018-5-24
That still won't work as it won't fit the boundary equation for E.

请先登录，再进行评论。

Matthew Hunt
Matthew Hunt 2018-5-25
So I managed a work around to get pdepe to work with my system and the code I'm using is:
%This is the simplistic model of a lithium ion battery model which I
%cobbled together knowing very little about batteries. The constants and so
%forth come from experiment.
%Define some global constants and parameters:
global rho; rho=1;global c_th; c_th=1;
global a_1; a_1=1;global a_2; a_2=1;global b; b=1; global c_1; c_1=1;global c_2; c_2=1;
global I_app; I_app=1.5;
global k;global D; global epsilon;
k=1; %Thermal conductivity, most likely a function.
epsilon=1; %Electrical conductivity, most likely a function
D=1; %lithium ion diffusion, most likely a function
m=0;
x=linspace(0,1,100);
t=linspace(0,1,50);
%Initial conditions;
global T_0; global c_0; global phi_0;
T_0=x./x;
c_0=smooth_step(x,0,1/3);
E_0=I_app/(epsilon)-(b/epsilon)*cumtrapz(x,c_0);
phi_0=cumtrapz(x,E_0);
sol = pdepe(m,@battery_GE,@battery_ic,@battery_bc,x,t);
T = sol(:,:,1);
E = sol(:,:,2);
c = sol(:,:,3);
figure;
surf(x,t,T);
title('Temperature');
xlabel('Distance x');
ylabel('Time t');
figure;
surf(x,t,E);
title('Electric Field');
xlabel('Distance x');
ylabel('Time t');
figure;
surf(x,t,c);
title('Li ion concentration');
xlabel('Distance x');
ylabel('Time t');
function [CC, FF, SS]=battery_GE(x,t,u,DuDx)
CC=[rho*c_th; 0; 1];
FF=[a_1*k; epsilon; c_1*D].*DuDx;
SS=[a_2*sigma*DuDx(2); -b*u(3);c_2*(DuDx(3)-(b*u(3)/epsilon)-DuDx(1))];
function u0=battery_ic(x)
u0=[T_0;phi_0;c_0];
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
T_L=1;
T_R=2;
pl=[T_L; I_app; 0];
ql=[0; sigma/epsilon 1];
pr=[T_R; -I_app; 0];
qr=[0; sigma/epsilon 1];
I get the following error, which I don't understand because I made T_0 a global variable:
Undefined function or variable 'T_0'.
Error in battery_model>battery_ic (line 54)
u0=[T_0;phi_0;c_0];
Error in pdepe (line 229)
temp = feval(ic,xmesh(1),varargin{:});
Error in battery_model (line 24)
sol = pdepe(m,@battery_GE,@battery_ic,@battery_bc,x,t);
I'm not sure what is going on now.
  2 个评论
Matthew Hunt
Matthew Hunt 2018-5-25
function y=smooth_step(x,a,b)
%This is a smooth approximation of the Step function done via the erf
%function
z=zeros(1,length(x));
if (a==x(1))
m=find(x>b,1,'first');
z(1:m-1)=1;
y=min(z,0.5*(1+erf(-30*(x-b))));
elseif(b==x(end))
n=find(x<=a,1,'last');
z(n+1:end)=1;
y=min(z,0.5*(1+erf(30*(x-a))));
else
y=0.5*min(1+erf(30*(x-a)),1+erf(-30*(x-b)));
end

请先登录，再进行评论。

Precise Simulation
Systems of PDEs might also be easier to solve with the FEATool FEM Toolbox which features a GUI and easy syntax for defining custom PDEs and equations. This example of heat transport and diffusion might be a good start, otherwise you could perhaps use the built-in convection and diffusion physics mode and modify it according to your equations.
Matthew Hunt
Matthew Hunt 2018-6-15
I have some values which allow the code to run BUT I have oscillations in the diffusion which I should not get. I am assuming that this is a numerical instability, is it caused by too large a timestep?

类别

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

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by