Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I was trying to solve an diffusion-reaction rate model for the oxidation of Silicon (based on the nice paper by Almeida et.al.) with the help of NDSolve. In this model the boundary conditions are a bit unusual compared to the way Mathematica expects them. The command I tried looks like the following:

res = NDSolve[{
 D[nO2[x, t], {t, 1}] == D0 Exp[-(ED/(kB T))] D[nO2[x, t], {x, 2}] - 
                         K0 Exp[-(EK/(kB T))] nO2[x, t] nSi[x, t],
 nO2[0.0, t] == nO20[1, T],
 D[nSi[x, t], {t, 1}] == -K0 E^(-(EK/(kB T))) nO2[x, t] nSi[x, t],
 nSi[x, 0.0] == 1.0
 } //. {
 D0 -> 1.291 10^11, ED -> 2.22,
 K0 -> 2.022 10^7, EK -> 1.42, 
 nSi0 -> 0.91 10^15, 
 kB -> QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"/"ElementaryCharge"]]],
  T -> 700 + 273.15
 }, {nO2, nSi}, {x, 0, 3000}, {t, 0, 60 60}][[1]]

where Mathematica complains about the way how the boundary condition for the function nO2[x,t] has been defined and drops an NDSolve::ivone error.

Playing around with this boundary condition I came up with a setup which gave solutions to the functions nO2[x,t] and nSi[x,t], however the choice of the x coordinate for the upper boundary where nO2 is zero is causing a strong influence on the solution as one can see when NDSolve is run for two different settings for the x coordinate:

res = NDSolve[{
 D[nO2[x, t], {t, 1}] == D0 Exp[-(ED/(kB T))] D[nO2[x, t], {x, 2}] - 
                         K0 Exp[-(EK/(kB T))] nO2[x, t] nSi[x, t],
 nO2[x, 0.0] == Piecewise[{{0.161/T, x == 0}}, 0],
 nO2[3000.0, t] == 0.0,
 nO2[0.0, t] == nO20[1, T],
 D[nSi[x, t], {t, 1}] == -K0 E^(-(EK/(kB T))) nO2[x, t] nSi[x, t],
 nSi[x, 0.0] == 1.0
 } //. {
 D0 -> 1.291 10^11, ED -> 2.22,
 K0 -> 2.022 10^7, EK -> 1.42, 
 nSi0 -> 0.91 10^15, 
 kB -> QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"/"ElementaryCharge"]]],
  T -> 700 + 273.15
 }, {nO2, nSi}, {x, 0, 3000}, {t, 0, 60 60}][[1]]
Plot[1 - (nSi /. res)[x, 60 60], {x, 0, 30}, PlotRange -> {0, 1.1}]

yields the following Graph:

nSiO2 with nO2===0 at x=3000

whereas the following code

res = NDSolve[{
 D[nO2[x, t], {t, 1}] == D0 Exp[-(ED/(kB T))] D[nO2[x, t], {x, 2}] - 
                         K0 Exp[-(EK/(kB T))] nO2[x, t] nSi[x, t],
 nO2[x, 0.0] == Piecewise[{{0.161/T, x == 0}}, 0],
 nO2[300.0, t] == 0.0,
 nO2[0.0, t] == nO20[1, T],
 D[nSi[x, t], {t, 1}] == -K0 E^(-(EK/(kB T))) nO2[x, t] nSi[x, t],
 nSi[x, 0.0] == 1.0
 } //. {
 D0 -> 1.291 10^11, ED -> 2.22,
 K0 -> 2.022 10^7, EK -> 1.42, 
 nSi0 -> 0.91 10^15, 
 kB -> QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"/"ElementaryCharge"]]],
  T -> 700 + 273.15
 }, {nO2, nSi}, {x, 0, 300}, {t, 0, 60 60}][[1]]
Plot[1 - (nSi /. res)[x, 60 60], {x, 0, 30}, PlotRange -> {0, 1.1}]

gives the following graph

nSiO2 with nO2===0 at x=300

It's obvious that one gets different results by moving the position of the boundary condition and thus using the boundary condition in this way is not an option. I would be interested in finding a possibility to specify the boundary condition for nO2 as shown in the first code snippet.

share|improve this question

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.