NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg59629] NDSolve
- From: Daniel Grady <d.c.grady at gmail.com>
- Date: Sun, 14 Aug 2005 04:38:13 -0400 (EDT)
- Reply-to: dcgrad at wm.edu
- Sender: owner-wri-mathgroup at wolfram.com
I've been having some trouble using NDSolve with a problem I've been
working on. I'm trying to simulate the temperature change in a solid
when it's hit by a laser pulse. The theory is that you keep track of
two different temperatures; one for the electron 'soup' in the
material, and another for the lattice in which the electrons are
embedded. When a laser pulse hits the material, all the energy goes
into the electrons, whose temperature rises very quickly, and then is
transferred to the lattice, and the heat then diffuses throughout the
lattice. The problems I've been running in to appear to be related to
the fact that the timescales that the laser pulse and diffusion
operate on are quite different; the duration of one laser pulse is
about 150 femtoseconds, and I'm interested in looking at the effects
of diffusion over the course of a millisecond or so.
Below I've pasted the model I've been using. There are two partial
differential equations; the first one is for the electrons and says
that electron specific heat times derivative of electron temperature
with respect to time equals the electron/lattice coupling term plus
the laser source term; the second one is for the lattice, and says
lattice specific heat time the derivative of lattice temperature with
respect to time equals the lattice diffusion term plus the
electron/lattice coupling constant. I'm only considering diffusion in
one dimension for the time being, as it's turned out to be plenty
complicated.
Because of the huge difference in timescales, the way I've been trying
to set this up is to first get a solution over a duration that's about
three times as long as the pulse itself; I then pull out the final
temperature profiles for the electrons and lattice, and use those as
initial conditions for another instance of NDSolve that gets a
solution over a millisecond duration.
Using the numbers included below, Mathematica will give me a precision
warning when it executes the first NDSolve command, although the
graphs it returns look very reasonable. However, when trying to
execute the second NDSolve command, it gives me precision warnings,
says boundary/initial conditions are inconsistent, and returns some
graphs that are clearly not accurate. I'm puzzled by this, since once
the pulse has ended, the temperature functions for both the electrons
and lattice shouldn't be behaving strangely at all; it should be a
smooth decrease. In addition, I'm pretty confident that the model
itself is good, as I have tried using smaller numbers for all of the
constants and everything works out beautifully.
I apologize for the length of this question, but this problem has been
driving me crazy for the last week or so, and no one on my campus uses
Mathematica. If anyone has any suggestions at all, I would greatly
appreciate them.
-D. Grady
The code follows. In this model,
id is the depth that the laser pulse penetrates the material
tau is the pulse duration
gamma, g, rho, cl, and lambda are all physical constants (for nickel,
in this case)
t0 is room temperature
dn is just an integer multiplier for the impact depth (id)
tn is the same as dn, but for pulse duration
\!\(\(\[Gamma] = 6*^3;\)\n
g = 8*^17; \n
\[Rho] = 89.08; \n
cl = 2.5*^4; \n
\[Lambda] = 90.9; \n
t0 = 300; \n\[IndentingNewLine]
amp = 1*^23; \n
\[Tau] = 150*^-15; \n\[IndentingNewLine]
id = 100*^-9; \n\[IndentingNewLine]
dn = 10; \n
tn = 3; \n\[IndentingNewLine]
sol = Flatten[\[IndentingNewLine]NDSolve[\[IndentingNewLine]{\[Gamma]*
e[z, t]*\[PartialD]\_t e[z,
t] \[Equal] \(-g\)*\((e[z, t] - l[z, t])\) +
If[\((0 \[LessEqual]
t \[LessEqual] \[Tau])\) && \((0 \[LessEqual]
z \[LessEqual] id)\), amp,
0], \[IndentingNewLine]\[Rho]*
cl*\[PartialD]\_t\
l[z, t] \[Equal] \[Lambda]*\[PartialD]\_\(z, z\)l[z, t] -
g*\((l[z, t] - e[z, t])\), \[IndentingNewLine]e[z, 0] \[Equal]
t0, l[z, 0] \[Equal] t0, \[IndentingNewLine]e[dn*id, t] \[Equal]
t0, l[dn*id, t] \[Equal]
t0, \[IndentingNewLine]\(\(Derivative[1, 0]\)[l]\)[0, t] \[Equal]
0}, \[IndentingNewLine]{e, l}, {z, 0, dn*id}, {t, 0,
tn*\[Tau]}\[IndentingNewLine]]\[IndentingNewLine]]; \n\
\[IndentingNewLine]
Plot3D[Evaluate[e[z, t] /. sol], {t, 0, tn*\[Tau]}, {z, 0, dn*id},
PlotRange \[Rule] All]; \n
\(Plot3D[Evaluate[l[z, t] /. sol], {t, 0, tn*\[Tau]}, {z, 0, dn*id},
PlotRange \[Rule] All];\)\n
\(e1[z_] := Evaluate[e[z, tn*\[Tau]] /. sol];\)\n
l1[z_] := Evaluate[l[z, tn*\[Tau]] /. sol]; \n
\(Plot[{e1[z], l1[z]}, {z, 0, dn*id}];\)\[IndentingNewLine] (*\
No\ Pulse\ *) \n
\(maxTime = .001;\)\[IndentingNewLine]
sol2 = Flatten[\[IndentingNewLine]NDSolve[\[IndentingNewLine]{\[Gamma]*
e[z, t]*\[PartialD]\_t e[z, t] \[Equal] \(-g\)*\((e[z, t] -
l[z, t])\), \[IndentingNewLine]\[Rho]*
cl*\[PartialD]\_t\
l[z, t] \[Equal] \[Lambda]*\[PartialD]\_\(z, z\)l[z, t] -
g*\((l[z, t] - e[z, t])\), \[IndentingNewLine]e[z, 0] \[Equal]
e1[z], l[z, 0] \[Equal]
l1[z], \[IndentingNewLine]e[dn*id, t] \[Equal] t0,
l[dn*id, t] \[Equal]
t0, \[IndentingNewLine]\(\(Derivative[1, 0]\)[l]\)[0, t] \[Equal]
0}, \[IndentingNewLine]{e, l}, {z, 0, dn*id}, {t, 0,
maxTime}\[IndentingNewLine]]\[IndentingNewLine]]; \
\[IndentingNewLine]\[IndentingNewLine]
Plot3D[Evaluate[e[z, t] /. sol2], {t, 0, maxTime}, {z, 0, dn*id},
PlotRange \[Rule] All]; \n
\(Plot3D[Evaluate[l[z, t] /. sol2], {t, 0, maxTime}, {z, 0, dn*id},
PlotRange \[Rule] All];\)\)