Solving 2D Convection Diffusion Equation

174 views (last 30 days)
Aimi Oguri on 14 Nov 2019
Edited: Aimi Oguri on 5 Dec 2019
I want to solve the above convection diffusion equation.
First, I tried to program in 1D, but I can't rewrite in 2D. I refered to here.
Can anybody help me?
function ConvectionDiffusion
m = 0;
x = linspace(0,0.025,10); %0~0.025の間を10個に等間隔でとるベクトル
t = linspace(0,10);
DATP = 2.36*10e-10; %m^2/s
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);
h = 2*10e-4; %m
vvar = h*tau/6*mu;
y = 6*10e-6;
vy = 6*vvar*y/h*(1-y/h); %velocity in x direction
function [g,f,s] = pdefun(x,t,c,DcDx)
v = vy;
g = 1;
f = D*DcDx;
s = -v*DcDx;
end
function c0 = icfun(x)
c0 = 0;
end
function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
pl = -10*D;
ql = 1;
pr = 0;
qr = 1;
end
end

Ravi Kumar on 14 Nov 2019
You might be able to setup the 2-D version using PDE toolbox. Take look at the following code. Please verify I have accuratly captured the problem statement.
function convDiff
model = createpde;
L = 0.025; %L
h = 2E-4; %h
vb = 1; % Define the correct value.
gdm = [3;4;0;L;L;0;0;0;h;h];
g = decsg(gdm,'R1',('R1')');
geometryFromEdges(model,g);
generateMesh(model,'Hmax', h/4);
pdegplot(model,'EdgeLabels','on');
D = 2.36E-10;
specifyCoefficients(model,'m',0,'d',0,'c',D,'a',0,'f',@fcoef);
applyBoundaryCondition(model,'dirichlet','Edge',4,'u',0);
applyBoundaryCondition(model,'neumann','Edge',1,'g',@gcoef);
applyBoundaryCondition(model,'neumann','Edge',3,'g',0);
R = solvepde(model)
pdeplot(model,'XYData',R.NodalSolution)
function g = gcoef(location,state)
Vmax = 0.8E-6;
Km = 0.457;
k = Vmax/Km;
mu = 9.45E-4;
tauW = 6*mu*vb/h;
a1 = 9.281E-11;
a2 = 1.505E-9;
a3 = 5.624;
S = a1 + a2*tauW/(a3+tauW);
g = k.*state.u - S;
end
function f = fcoef(location,state)
yb = location.y/h;
v = 6*vb*yb.*(1-yb);
f = v.*state.ux;
end
end
Regards,
Ravi

Aimi Oguri on 15 Nov 2019
What should I see if I want to know C at y=0?
Ravi Kumar on 15 Nov 2019
You can find all nodes on the edge 1, which is at y = 0 and plot the value of C by indexing into R.NodalSolution:
y0Nodes = findNodes(model.Mesh,'region','Edge',1);
figure
plot(model.Mesh.Nodes(1,y0Nodes),R.NodalSolution(y0Nodes))
Aimi Oguri on 16 Nov 2019
I understand.
Thank you so much!