4
$\begingroup$

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]

It gives the following error: enter image description here

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: enter image description here

However, at t=10 and near x=0, the place where I enforced the boundary condition, the solutions diverges: enter image description here

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.
$\endgroup$
10
  • $\begingroup$ $\partial_x(\mu(T,y)\partial_x u)$ actually involves $\partial y/\partial x$, so you still need one b.c. for $y$. $\endgroup$
    – xzczd
    Jan 3 at 12:31
  • $\begingroup$ Dear @xzczd, when I try to set up homogeneous neumann boundary conditions for 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$
    – Eugene W.
    Jan 3 at 13:41
  • $\begingroup$ That's expected. For $\partial y/\partial x$ you can only set Dirichlet b.c. $\endgroup$
    – xzczd
    Jan 3 at 14:49
  • $\begingroup$ @xzczd What does "For ∂y/∂x you can only set Dirichlet b.c. "? It is not physical to set Dirichlet b.c. for y[t,x] for both planes (x=0 and x=1). Do you mean to define ∂y/∂x=0 at x=0 and x=1? I have tried it at x=1: using Derivative[0, 1][y][t, 1] == 0 $\endgroup$
    – Eugene W.
    Jan 3 at 14:58
  • $\begingroup$ As the warning said, "Boundary condition should have derivatives of order lower than the differential order of the partial differential equation." Just think about the simple $\frac{\partial u}{\partial t} - c\frac{\partial u}{\partial x}=0$, if you set a Neumann b.c. for it, it'll generally lead to contradiction between b.c. and PDE. (Of course, there exist exceptions e.g. this, but that's another story. ) So, you need one Dirichlet b.c. for $y$ at either $x=0$ or $x=1$. $\endgroup$
    – xzczd
    Jan 3 at 15:12
2
$\begingroup$

This problem can be solved with using additional boundary condition for y and by increasing order of equation for y by differentiation $y_t$ on x 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; 

Boundary condition for y[t,x] at x=0

sol0 = 
 DSolve[{y0'[t] == ((1 - y0[t])^n)*A*Exp[-Ea/(T0)], y0[0] == 0}, 
  y0[t], t];

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, x] == D[((1 - y[t, x])^n)*A*Exp[-Ea/(T[t, x])], x], 
    y[0, x] == 0., y[t, 0] == y0[t] /. sol0[[1]], 
    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" -> Automatic, 
     "PDEDiscretization" -> mol[137, 8]}];
    

Visualization (note, we use $\log_{10}\mu$)

p1 = Plot3D[Evaluate[u[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, 
   PlotRange -> All, AxesLabel -> {"t", "x", "u"}, Mesh -> None, 
   ColorFunction -> "Rainbow"];
p2 = Plot3D[Evaluate[T[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, 
   PlotRange -> All, AxesLabel -> {"t", "x", "T"}, Mesh -> None, 
   ColorFunction -> "Rainbow"];
p3 = Plot3D[Evaluate[y[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, 
   PlotRange -> All, AxesLabel -> {"t", "x", "\[Alpha]"}, 
   Mesh -> None, ColorFunction -> "Rainbow"];
p4 = Plot3D[
   Evaluate[Log[10, mu[T[t, x] /. sol, y[t, x] /. sol]]], {t, 0, 
    tmax}, {x, 0, 1}, PlotRange -> All, 
   AxesLabel -> {"t", "x", 
     "\!\(\*SubscriptBox[\(Log\), \(10\)]\)\[Mu]"}, Mesh -> None, 
   ColorFunction -> "Rainbow"];
GraphicsGrid[{{p1, p2}, {p3, p4}}, ImageSize -> 700]

Figure 1

We can compare sol with numerical solution computed with code @Eugene W.:

sol1 = 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" -> Automatic, 
     "PDEDiscretization" -> mol[137, 8]}];

In this case we have a message

NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution.

Nevertheless both solutions are consider for every t and x

Table[{Plot[{y[t, x] /. sol, y[t, x] /. sol1}, {x, 0, 1}, 
   AxesLabel -> {"x", "y"}, PlotLegends -> Automatic, PlotLabel -> t],
   Plot[{u[t, x] /. sol, u[t, x] /. sol1}, {x, 0, 1}, 
   AxesLabel -> {"x", "u"}, PlotLegends -> Automatic, PlotLabel -> t, 
   PlotRange -> All], 
  Plot[{T[t, x] /. sol, T[t, x] /. sol1}, {x, 0, 1}, 
   AxesLabel -> {"x", "T"}, PlotLegends -> Automatic, PlotLabel -> t, 
   PlotRange -> All]}, {t, 1, 10, 3}]

Figure 2

$\endgroup$
5
  • $\begingroup$ It is an interesting trick, but such a boundary condition gives a totally different solution. If you see the equation for $y$, then you will see that at big time the solution $y$ approaches 1 in the whole domain. $\endgroup$
    – Eugene W.
    Jan 3 at 17:39
  • $\begingroup$ @EugeneW. Do you mean that we need some additional equation for y[t,0]? $\endgroup$ Jan 5 at 7:15
  • $\begingroup$ @EugeneW. Please, see update to my answer. $\endgroup$ Jan 5 at 16:24
  • $\begingroup$ Thank you! it is really a very smart trick! Thank you! $\endgroup$
    – Eugene W.
    Jan 6 at 14:19
  • $\begingroup$ @EugeneW. You are welcome! $\endgroup$ Jan 6 at 14:36

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.