MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Mathematica 7.01 and Mac OS 10.6
  • Next by Date: About binary relations
  • Previous by thread: Re: Mathematica 7.01 and Mac OS 10.6
  • Next by thread: About binary relations