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
