Solving 2D Convection Diffusion Equation

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

Answer 1

[Ravi Kumar]

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]

I appreciate your help.

I would like to ask a question additionally.

What should I see if I want to know C at y=0?

[Ravi Kumar]

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]

I understand.

Thank you so much!

[Emily Condon]

@Ravi Kumar Is there a way to write a coefficient function for the coefficient of 'c'? For example, if I wanted to have a coefficient Dx for the d2u/dx2 part of c, and Dy for the d2u/dy2, is there a way to write a sort of dot product function Dcoef? I know that the state fields only are available for first order dervatives, eg. state.ux

Thank you in advance!