NDSolveValue: The dependent variable in ... in the boundary condition DirichletCondition[...] needs to be linear

Question

I'm solving the following convection-diffusion steady-state problem with 1st order slip conditions at the walls:

velocity:

$\qquad\mu\frac{d^2u}{dy^2}=\frac{dp}{dx}$

$\qquad u(x,y=0)=u_w+\beta_u\lambda\left.\frac{du}{dy}\right|_{y=0}$

$\qquad\left.\frac{du}{dy}\right|_{y=h/2}=0$

temperature:

$\qquad u\frac{\partial T}{\partial x} = \alpha \left(\frac{\partial^2T}{\partial x^2}+\frac{\partial^2T}{\partial y^2}\right)$

$\qquad T(x,y=0)=T_w+\beta_T\lambda\left.\frac{\partial T}{\partial y}\right|_{y=0}$

$\qquad T(x,y=h)=T_w+\beta_T\lambda\left.\frac{\partial T}{\partial y}\right|_{y=h}$

$\qquad T(x=0,y)=T_0$

The movement equation has an analytical solution, whereas the temperature equation does not. Here's a minimum non-working example in Mathematica:

Clear["Global`*"]

(* numerical values *)
μ = 1.8*10^-5;
h = 10*10^-6;
L = 10*h;
λ = 3.39*10^-8;
G = 0.1/L;
α = 2.1*10^-5;
βu = 1.1739;
βT = 1.8922;
Tw = 60;
T0 = 20;
uw = 0;

(* 1st order solution *)
uSol = u == DSolve[{
        μ*D[u[y], y, y] + G == 0,
        u[0] == uw + βu*(λ*D[u[y], y] /. y -> 0),
        (D[u[y], y] /. y -> h/2) == 0
    }, u, y][[1]][[1]][[2]][[2]]
Plot[{uSol[[2]]} , {y, 0, h}, AxesLabel -> {"y", "u"}]

Tsol = NDSolveValue[{
        u*D[T[x, y], x] - α*(D[T[x, y], x, x] + D[T[x, y], y, y]) == 0 /. u -> uSol[[2]],
        T[x, 0] == Tw  + βT*(λ*D[T[x, y], y] /. y -> 0),
        T[x, h] == Tw  + βT*(λ*D[T[x, y], y] /. y -> h),
        T[0, y] == T0
    }, T, {x, 0, L}, {y, 0, h}]
DensityPlot[Tsol[x, y], {x, 0, L/10}, {y, 0, h}, ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, FrameLabel -> Automatic]

DSolve will handle the velocity equation just fine. However, NDSolveValue gives me the following error messages:

NDSolveValue::fembdnl: The dependent variable in T==60+6.41456*10^-8 (T^(0,1))[x,0] in the boundary condition DirichletCondition[T==60+6.41456*10^-8 (T^(0,1))[x,0],y==0.] needs to be linear.

NDSolveValue::dsvar: 7.142857142857143`*^-10 cannot be used as a variable.

NDSolveValue::dsvar: 7.15`*^-7 cannot be used as a variable.

General::stop: Further output of NDSolveValue::dsvar will be suppressed during this calculation.

I will also note that if I try to solve these equations with no-slip conditions (0th order) instead, by replacing these:
$\qquad u(x,y=0)=u_w+\beta_u\lambda\left.\frac{du}{dy}\right|_{y=0}$
$\qquad T(x,y=0)=T_w+\beta_T\lambda\left.\frac{\partial T}{\partial y}\right|_{y=0}$
with these:
$\qquad u(x,y=0)=u_w$
$\qquad T(x,y=0)=T_w$
which results in the following code:

(* 0th order solution *)
uSol = u == DSolve[{
        μ*D[u[y], y, y] + G == 0,
        u[0] == uw ,
        (D[u[y], y] /. y -> h/2) == 0
    }, u, y][[1]][[1]][[2]][[2]]
Plot[{uSol[[2]]} , {y, 0, h}, AxesLabel -> {"y", "u"}]

Tsol = NDSolveValue[{
        u*D[T[x, y], x] - α*(D[T[x, y], x, x] + D[T[x, y], y, y]) == 0 /. u -> uSol[[2]],
        T[x, 0] == Tw ,
        T[x, h] == Tw ,
        T[0, y] == T0
    }, T, {x, 0, L}, {y, 0, h}]
DensityPlot[Tsol[x, y], {x, 0, L/10}, {y, 0, h}, ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, FrameLabel -> Automatic]  

Mathematica will solve them just fine.

Any help would be appreciated.

Answer 1

With constants specified as in the question, uSol can be determined by

uSol = DSolveValue[{μ*D[u[y], y, y] + G == 0, u[0] == uw + βu*(λ*D[u[y], y] /. y -> 0), 
    (D[u[y], y] /. y -> h/2) == 0}, u[y], y]
Plot[uSol, {y, 0, h}, AxesLabel -> {"y", "u"},ImageSize -> Large, 
    LabelStyle -> Directive[Black, Bold, Medium]]

(* 0.0000110542 + 277.778 y - 2.77778*10^7 y^2 *)

This intermediate result is no different from that in the question, but the code is a bit simpler. The error messages quoted in the question arise from the first two boundary conditions for T. For instance,

T[x, 0] == Tw  + βT*(λ*D[T[x, y], y] /. y -> 0)

should be written as a generalized Neumann boundary condition. So recast it as

α D[T[x, y], y] - uSol T[x, y] == α (T[x, y] - Tw)/(βT*λ) - uSol T[x, y]

evaluated at y == 0. This is the standard form given in "Details" of the NeumannValue documentation. Hence, this boundary condition becomes

NeumannValue[-uSol T[x, y] + α (T[x, y] - Tw)/(βT*λ), y == 0]

and Tsol is determined by

Tsol = NDSolveValue[{uSol*D[T[x, y], x] - α*(D[T[x, y], x, x] + D[T[x, y], y, y]) == 
    NeumannValue[-uSol T[x, y] + α (T[x, y] - Tw)/(βT*λ), y == 0] - 
    NeumannValue[-uSol T[x, y] + α (T[x, y] - Tw)/(βT*λ), y == h],
    T[0, y] == T0}, T[x, y], {x, 0, L}, {y, 0, h}];
Plot3D[Tsol, {x, 0, L/10}, {y, 0, h}, PlotRange -> All, AxesLabel -> {x, y, T}, 
    ImageSize -> Large, LabelStyle -> Directive[Black, Bold, Medium]]

The value of the solution is determined primarily by T0 near x == 0 and Tw elsewhere, with the transient oscillations due to the large deference between T0 and Tw. The asymmetry between the solution near y == 0 and y == h occurs, because D[T[x, y], y], as specified in the question, points into the computational region at y == 0 and out at y == h. Reversing the sign of NeumannValue at y == h would symmetrize the solution.

Addendum: Resolution Issue

Increasing the resolution and reversing the sign of NeumannValue at y == 0 eliminates the oscillations seen above.

Tsol = NDSolveValue[{uSol*D[T[x, y], x] - α*(D[T[x, y], x, x] + D[T[x, y], y, y]) == 
    -NeumannValue[-uSol T[x, y] + α (T[x, y] - Tw)/(βT*λ), y == 0] 
    -NeumannValue[-uSol T[x, y] + α (T[x, y] - Tw)/(βT*λ), y == h],
    T[0, y] == T0}, T[x, y], {x, 0, L}, {y, 0, h}, Method -> 
    {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 10^-7}}}]

which still runs rapidly and yields

(The plotting option PlotRange -> {20, 76} was used for the last figure to ease comparison with the second plot above.) Note that not reversing the sign of NeumannValue at y == 0 yields a large but steady oscillation at that boundary.