Concentration dependent diffusivity in diffusion equation
- To: mathgroup at smc.vnet.net
- Subject: [mg103999] Concentration dependent diffusivity in diffusion equation
- From: Chris Scullard <scullard at uchicago.edu>
- Date: Thu, 15 Oct 2009 07:16:17 -0400 (EDT)
Hi all,
I am attempting to solve the diffusion equation for the concentration of
a compound (MgO), where the diffusivity depends on the concentration of
some other compound (SiO2). The concentration of SiO2 is itself
determined by a differential equation. I have "solved" the problem, but
the solution takes several hours when I think it ought to be much
faster. Being quite sure that my solution is just an ugly hack, I'm
hoping somebody has some advice on how to speed it up. The algorithm is
as follows:
Step 1:
Solve the 1-D diffusion equation for SiO2 concentration in a linear
domain with a step-function initial condition. About half of the domain
contains all the SiO2 and the other half has none initially. I use an
approximation to the step function because the solver does not seem to
like UnitStep[] in the initial conditions. Also, the diffusivity of SiO2
depends on its own concentration (given by the function DSiO2) so the
equation is not linear. I solve this part with the following:
D0 = 4*10^(-1);
D1 = 0;
alpha = 12;
DSiO2[XSiO2_] = D0*Exp[-alpha*(XSiO2/100 - 0.5)] + D1;
left = 0;
right = 6.05;
leftSiO2 = 71.68;
rightSiO2 = 55.67;
solSiO2 =
NDSolve[{D[u[x, t], t] ==
D[DSiO2[u[x, t]]*D[u[x, t], x],
x], (D[u[x, t], x] /. x -> left) ==
0, (D[u[x, t], x] /. x -> right) == 0,
u[x, 0] == (rightSiO2 - leftSiO2)/(1 + Exp[-100 (x - 3.025)]) +
leftSiO2}, u, {x, left, right}, {t, 0, 15},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
MaxPoints -> 500000}}]
This part is very fast, taking only a few seconds and the result is
reasonable.
Step 2:
Solve for MgO concentration, where the diffusivity depends on the
previously-found SiO2 concentration, given by the function DMgO:
SiO2model[x_, t_] = u[x, t] /. solSiO2
D0MgO = 4*10^(-1);
D1MgO = 0 ;
alphaMgO = 14.5;
DMgO[x_, t_] =
D0MgO*Exp[-alphaMgO*(First[SiO2model[x, t]]/100 - 0.5)] + D1MgO;
leftMgO = 0.85;
rightMgO = 5;
Timing[solMgO =
NDSolve[{D[u[x, t], t] ==
D[DMgO[x, t]*D[u[x, t], x], x], (D[u[x, t], x] /. x -> left) ==
0, (D[u[x, t], x] /. x -> right) == 0,
u[x, 0] == (rightMgO - leftMgO)/(1 + Exp[-100 *(x - 3.025)]) +
leftMgO}, u, {x, left, right}, {t, 0, 15},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
MaxPoints -> 50000}}]]
This step takes hours. I assume there is some problem with the way I am
incorporating the previous solution into the new equation, because I
don't think this should take so long. I would appreciate any input.
Thanks,
Chris