I am trying to solve the system of equations (momentum conservation, heat, and polymerization equations) in the Couette problem. There are two planes, the bottom is fixed, the top is moving with the speed 1. The bottom plane is heated with the temperature $T_0$, the top plane with $T_s$. At the initial time, I set a linear function for the velocity and the temperature, $0$ in the whole domain for the polymerization degree. Shortly, polymerization degree means the quality of curing of some polymer (for example, glue): $0$ means the "glue is fresh", $1$ means "glue is totally cured". Thus, $y$ always lays between $0$ and $1$.
Above process can be formulated by these system of equation: $$\partial_t u=\partial_x(\mu(T,y)\partial_x u),$$ $$\qquad\partial_t T=\kappa\Delta T + H_{tr}\partial_t y,$$ $$\qquad\partial_t y= A(1-y)^n\exp(-E_a/T),$$ $$\mu(T,y)=\mu_0 \exp(E_\eta/T + \xi y)$$
The Wolfram Mathematica code is given as follows:
tmax = 10;
kappar = 0.8853475;
n = 1.2;
A = 120.9;
Ea = 8;
Eeta = 15.075;
Htr = 0.942895;
mu0 = 0.00213189;
ksi = 20;
T0 = 1; Ts = 400./300.;
mu[T_, y_] := mu0*Exp[Eeta/(T) + ksi*y];
mol[n_Integer, o_Integer] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> o, "Coordinates" -> N[Range[0, n]/n]}}
TinitLinear[x_] := T0 + (Ts - T0)*x;
sol = NDSolve[{D[u[t, x], t] ==
D[mu[T[t, x], y[t, x]]*D[u[t, x], x], x],
D[T[t, x], t] == kappar*D[T[t, x], x, x] + Htr*D[y[t, x], t],
D[y[t, x], t] == ((1 - y[t, x])^n)*A*Exp[-Ea/(T[t, x])],
y[0, x] == 0., T[0, x] == TinitLinear[x], u[0, x] == x,
T[t, 0] == T0, T[t, 1] == Ts, u[t, 0] == 0., u[t, 1] == 1.}, {u,
T, y}, {t, 0, tmax}, {x, 0, 1},
Method -> {"TimeIntegration" -> {"BDF", "MaxDifferenceOrder" -> 5},
"PDEDiscretization" -> mol[101, 4]}];
p1 = Plot3D[Evaluate[u[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, PlotRange -> All, AxesLabel -> {"t", "x", "u"}];
p2 = Plot3D[Evaluate[T[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, PlotRange -> All, AxesLabel -> {"t", "x", "T"}];
p3 = Plot3D[Evaluate[y[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, PlotRange -> All, AxesLabel -> {"t", "x", "\[Alpha]"}];
p4 = Plot3D[
Evaluate[mu[T[t, x] /. sol, y[t, x] /. sol]], {t, 0, tmax}, {x, 0,
1}, PlotRange -> All, AxesLabel -> {"t", "x", "\[Mu]"}];
GraphicsGrid[{{p1, p2}, {p3, p4}}, ImageSize -> 700]
Why the wolfram gives me such a warning? I set up a bunch of boundary conditions:
- Since the system of equations is the first order in time, then I set up one initial condition per $u, T, y$,
- Since there is a second-order derivative in space ($u, T$), then I set up 2 Dirichlet boundary conditions
- Since the polymerization equation is ODE, then no need to set up boundary conditions in space
Could you give me a hint, where is the wrong logic? Or is it a wolfram NDSolve
bug?
I searched similar problems in the StackExchange, but I found the cases where really insufficient boundary conditions: 1,2, 3. In my case: I set up all boundary and initial conditions, and I still get a suspicious error.
Update
I tried @xzczd's trick and set up the Dirichlet boundary condition:
y[t, 0] == 0
Such a trick get rid of the warning and mainly it gives a similar solution:
However, at t=10
and near x=0
, the place where I enforced the boundary condition, the solutions diverges:
Also, I tried y[t, 1] == 1
and it gives me similar poor results.
Note
At big time $t\approx 20$ the polymerization degree $y(t,x)$ should approach 1. The approximate time can be estimated if you consider the isothermal case near corresponding walls with $T=T_s$ or $T=T_0$.
$$ \int_0^{y_*} \frac{dy}{(1-y)^n}=A\tau\exp(-E_a/T) = -\frac{(1-y_*)^{1-n} - 1}{1-n}$$
- If you take $y_*=0.5$, then you will get $\tau=2.5, 18.3$ for hot and cold walls,
- If you take $y_*=0.9$, then you will get $\tau=9.7, 72$ for hot and cold walls.
y
:Derivative[0, 1][y][t, 1] == 0
. I get an error:NDSolve::bdord: Boundary condition (y^(0,1))[t,1] should have derivatives of order lower than the differential order of the partial differential equation.
$\endgroup$y[t,x]
for both planes (x=0
andx=1
). Do you mean to define ∂y/∂x=0 atx=0
andx=1
? I have tried it atx=1
: usingDerivative[0, 1][y][t, 1] == 0
$\endgroup$