3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey Oct 8 '17 at 4:43
2
$\begingroup$

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 *)

enter image description here

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]]

enter image description here

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

enter image description here

(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.

$\endgroup$
  • 1
    $\begingroup$ Thank you for your help, @bbgodfrey! I will leave one additional piece of information here for anyone who might come upon this thread with a similar problem. The reason one can't specify Neumann boundary conditions as I did is here and I quote: "It is not possible to derive (unambiguously) from the Neumann value with which PDE equation the value should be associated. Making a NeumannValue part of a PDE equation solves this problem without ambiguity." $\endgroup$ – Leonardo Oct 9 '17 at 2:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.