I am trying to solve a coupled system of PDEs for 2 functions h[t,r]
and c[t,r]
, with initial conditions of h[0,r]
as a normal distribution over r
and c[0,r]
set to a fixed value of 0.42
(c[r,t]
is supposed to always be bounded between 0
and 1
), and fixed boundary values for h[t,r]
at two points r = 0.1
and r = -0.1
(solve domain boundaries). However, when I put the equations and boundary conditions into NDSolve
, I get the warning
Dif = 10^(-5);
cinitial = 0.42;
μ = 0.10;
H = 0.003;
k = 0.001;
boundary = 0.1;
u = H*0.025*D[c[t, r], r]/(μ*c[t, r]);
NDSolve[{k r (1 - c[t, r]) c[t, r] + h[t, r] Dif D[c[t, r], r] +
r D[h[t, r], r] Dif D[c[t, r], r] +
r h[t, r] Dif D[c[t, r], {r, 2}] ==
r h[t, r] (D[c[t, r], t] + u D[c[t, r], r]),
r D[h[t, r], t] + (r u D[h[t, r], r] + D[u, r]*h[t, r]) == k c[t, r] r,
h[0, r] == E^(-(r^2/(2 0.02^2)))/(0.02 10 Sqrt[2 π]),
c[0, r] == cinitial,
h[t, boundary] == E^(-(boundary^2/(2 0.02^2)))/(
0.02 10 Sqrt[2 π]),
h[t, -boundary] == E^(-(boundary^2/(2 0.02^2)))/(
0.02 10 Sqrt[2 π])}, {h[t, r], c[t, r]}, {t, 0,
10}, {r} ∈ ImplicitRegion[r^2 <= boundary^2, {r}]]
NDSolve::femcnsd: The PDE coefficient
0. -0.001 r c[r]+0.001 r c[r]^2-0.00001 h[r] c'[r]+(0.00075 r h[r] c'[r]^2)/c[r]-0.00001 r c'[r] h'[r]
does not evaluate to a numeric scalar at the coordinate{-0.1}
; it evaluated toIndeterminate
instead.
I am really not sure if the cause is my equations or NDSolve
. I have checked the equations and boundary conditions and was unable to find any errors, and I do not have enough understanding of NDSolve
and Mathematica. Is anyone able to identify the problem? Thanks!
I am using Mathematica 12.
Update: background info
The equations describe a droplet of water-alcohol mixture spreading on a surface due to a surface tension gradient which is the result of evaporation (concentration at the edge of droplet decreases faster, surface tension ). The first equation is the diffusion-convection equation for alcohol transport in the droplet, while the second equation is the continuity equation accounting for evaporative loss. u
is the radial velocity of each fluid element in the drop.
The model is 1D axisymmetric. The initial condition for h
is defined to be non-zero for all values of r
to make the solving easier, or an additional moving boundary will have to be added into the equation for the radius of spreading of the droplet. The boundary conditions for h
are equivalent to setting the bc that the height of the droplet at a radial coordinate approaching infinity to be 0
. I think it is also possible to include the condition D[c[r,t],r]
at r = 0.1
and -0.1
(again at the radial boundaries) to be 0
, if required.
With
to rewrite the code (to something likeWith[{u=u[t,x]},D[u,t]+D[u,x]==0]
) to make it easier to read. 3. The highest differential order ofc
inr
direction is2
, and that ofh
is1
, usually this suggests you need2
b.c. forc
and1
b.c. forh
, are you sure your b.c.s are correct? $\endgroup$