On Feb 26, 2008, at 7:16 AM, Etienne Rivard wrote:
After looking at the classes of various terms in the examples, I
couldn't find a way to represent $\frac{\partial}{\partial \xi}$.
The module examples.phase.impingement.mesh40x1.input has a similar
term but it uses the absolute, or magnitude of theta and this is
not exactly what I have. I have tried removing the getMag() part so
that I define theta.getGrad(), but that doesn't work, it's not the
right class.
It's clearer (to me, at least) to represent the spatial derivatives
consistently, either as $\frac{\partial}{\partial \xi}$ or as $\nabla
$, but not mixed. My preference is for $\nabla$, because it
highlights the vector relationships that are tripping you up, I
think. E.g.,
\[ \xi^2 \frac{\partial}{\partial t} (cy_F) - \xi^3 \frac{\dot{R}}{R}
\cdot \nabla (cy_F) + \frac{1}{R} \nabla \cdot
(\xi^2 c \vec{v} y_F) = \frac{1}{R^2} \nabla \cdot (\xi^2 c \bar{D}
\nabla y_F) \]
\[ \xi^2 \frac{\partial}{\partial t} (cy_F \theta) - \xi^3
\frac{\dot{R}}{R} \cdot \nabla (cy_F \theta) +
\frac{1}{R} \nabla \cdot (\xi^2 c \vec{v} y_F \theta) = \frac{1}
{R^2}
\nabla \cdot (\xi^2 c \tilde{D} y_F \nabla \theta) + \frac{1}{R^2}
\nabla \cdot (\xi^2 c \tilde{D} \theta \nabla y_F) \]
\[ \xi^2 \frac{\partial}{\partial t} (cy_F \Psi) - \xi^3
\frac{\dot{R}}{R} \cdot \nabla (cy_F \Psi) +
\frac{1}{R} \nabla \cdot (\xi^2 c \vec{v} y_F \Psi) = \frac{1}
{R^2} \nabla
\cdot (\xi^2 c \hat{D} y_F \nabla \Psi) + \frac{1}{R^2} \nabla \cdot
(\xi^2 c \hat{D} \Psi \nabla y_F) \]
so that you can see that the second term, resulting from the material
derivative, is a velocity vector $\xi \frac{\dot{R}}{R}$ dotted with
a field gradient $\nabla (cy_F)$, etc. When you multiplied
theta.getGrad() times the scalar quantity xi**3 * Rdot / R, you got a
vector source, which is incompatible with your scalar governing
equation. As far as FiPy is concerned, even in 1D, a scalar field and
a vector field are different, even if they hold the same values. So,
the first step is to make a vector variable to hold your velocity,
and then dot that with the gradient. Assuming you're using FiPy 1.2,
the simplest (but not best) way is:
xi = mesh.getCellCenters()[...,0]
velocity = VectorCellVariable(mesh=mesh)
velocity[...,0] = xi * Rdot / R
and then represent the second term of the first equation as
xi**2 * velocity.dot(c * yF.getGrad())
I didn't set the "velocity" to xi**3 * Rdot / R because it isn't
dimensionally accurate, but you might find that it makes things a
little easier.
Although the dot form will work, we don't recommend it. It's a first-
order difference and, depending on magnitude, it can destabilize your
equations. In my own work, I've made the transformation $\vec{v}
\cdot \nabla \phi \rightarrow \nabla \cdot \vec{v} \phi + \phi \cdot
\nabla \cdot \vec{v}$, which allows you to use an implicit convection
term and an implicit source term, instead of an explicit source. To
do this, you'll need to evaluate "velocity" on the faces:
xiF = mesh.getFaceCenters()[...,0]
velocityF = VectorFaceVariable(mesh=mesh)
velocityF[...,0] = xiF * Rdot / R
and then represent the second term as
PowerLawConvectionTerm(coeff=c * xiF**2 * velocity) -
ImplicitSourceTerm(coeff=c * (xiF**2 * velocity).getDivergence())
Please let me know if I've misunderstood what you're trying to do, or
if this doesn't work for you.
If you want to transform your results in $\xi$ back into the $r$
frame, we did some similar transformations between laboratory and
reference frames in a recent paper
@Article{Boettinger:2007p33,
author = {W. J Boettinger and J. E Guyer and C. E Campbell and G. B
McFadden},
journal = {P R Soc A},
pages = {3347--3373},
volume = {463},
year = {2007},
month = {Jan},
language = {English},
keywords = {Similarity Solution, Diffusion, Plane, Deformation,
Interdiffusion, Migration, Boltzmann-Matano, Phase Field, Dictra},
doi = {10.1098/rspa.2007.1904},
URL = {http://journals.royalsociety.org/content/n37t3346g7g17787/},
}
We included the heavily-commented FiPy script that we used for the
phase field portion in the electronic supplementary material.
Let me know if you can't obtain it and I'll send you a copy.
