How do I solve this 1D transient convection-diffusion equation with the convection term coupled with transient boundary values?

The following equation is a non-dimensionlized 1D transient convection diffusion equation, where tau and eta are dimensionless time and y-axis from 0 to infinity for both. Because the boundary value phi_m is coupled in the convection term, I have some difficulty to use MATLAB PDE solver. Can anybody help me to figure out?

Answer 1

[Bill Greene]

I have created a PDE solver (pde1dm) that is similar to pdepe but includes some enhancements such as the capability to add some scalar equations to the set of PDE. In your example, one of these scalar equations can track the value of ϕ at the left end and this can be used in defining the PDE.

If you want to try pde1dm, it can be downloaded here. My MATLAB code for solving your example with pde1dm is shown below.

function matlabAnswers_2_21_2023
phi_p=1e-3;
phi_star=2.22;
Nx=20;
xInf=4;
x = linspace(0,xInf,Nx);
tf=1;
t = linspace(0,tf,200);
xOde = 0.0; % ODE at left end
%% solve pde
m = 0;
odeicf = @() ode_IC;
pdef=@(x,t,u,dudx,v) pde_F(x,t,u,dudx,v,phi_star);
pdebc=@(xl,ul,xr,ur,t) pde_BC(xl,ul,xr,ur,t, phi_p, phi_star);
[sol,odeSol] = pde1dm(m,pdef,@pde_IC,pdebc,x,t,@ode_F, odeicf,xOde);
figure; plot(x, sol(end,:)); grid; xlabel("eta");
title("phi at final time");
figure; plot(t, sol(:,1)); grid; xlabel("time");
title("phi at left end (phi_m)");
figure; plot(t, odeSol); grid; xlabel("time");
end

%% functions
function [c,f,s] = pde_F(x,t,u,dudx,v,phi_star) % PDE to be solved
phi_m=v;
c = 1;
f = dudx;
s = (phi_star-phi_m)*dudx;
end
% ---------------------------------------------
function u0 = pde_IC(x) % Initial conditions of PDE
u0=1;
end
% ---------------------------------------------
function [pl,ql,pr,qr] = pde_BC(xl,ul,xr,ur,t, phi_p, phi_star) % BCs
phi_m=ul;
pl=(phi_m-phi_star)*phi_p + (phi_star-phi_m)*phi_m;
ql=1;
qr=0;
pr=ur-1;
end

function R=ode_F(t,v,vdot,x,u,DuDx,f, dudt, du2dxdt) % ODE to be solved 
R= u-v;
end

function value = ode_IC() % initial condition of ODE
value = pde_IC(0);
end

[Albert Kim]

Dear Bill

I don't have any language to thank you enough.

Let me try to understand your codes and contact you.

This is my first time asking help from MATLAB user forum and got wonderful experience.

Albert

Answer 2

[Torsten]

MATLAB's "pdepe" doesn't allow you to access phi_m in all grid points.

So you will have to follow one of the ways described here:

https://de.mathworks.com/matlabcentral/answers/1912595-how-to-integrate-over-the-entire-bounds-of-a-solution-to-a-pde-at-a-given-time-step-using-pdepe?s_tid=srchtitle

And the gradient d(phi)/d(eta) in your PDE is also meant to be evaluated at 0, I guess, not at eta ?

[Albert Kim]

Dear Torsten

Thank you so much. I wanted to know, as you said, whether it is possible to solve it using MATLAB's default PDE solver, which seems to be for constant coefficients.

If you question is about the third equation from the bottom, my answer is yes. The gradient or slope of phy with eta should be evaluated at eta = 0. which is a Robin boundary condition.

I had a similar thought that I have to solve three ODEs simultaneously, instead of solving 1D transient PDE.

I truly appreciate your comment. Is there anyway that I can reach you through email?