Posts for the month of December 2009

Simulations with van Leer convection term

The upwind convection term appeared to help the low concentration simulation work. However, the upwind convection term is inaccurate though. Let's use the van Leer convection term instead.

ID Notes Issues Status $N$ movie
267 base none benson, hobson, alice (stopped) 12
268 $X_1 = 0.1 X_1^{equ}$ none renfield, sancho, kato (stopped) 12
269 $\nu_f = 2 \times 10^{-2}$ crashes luggage (stopped) 12
270 $X_1 = 0.5 X_1^{equ}$ none luggage (stopped) 12
  • Posted: 2009-12-22 13:22 (Updated: 2010-01-14 10:47)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations with a negligible radial displacement

The radial displacement issue has been troubling me. Most of the simulations have had a radial displacement of dr / 10 and some others have dr. Some simulations don't run at all and an instability can clearly be seen to emanate from the r=0 location. I tested the base simulation with dr * 1e-10 just to figure out why I had trouble in the past. Needless to say the simulation seem to be running fine. I have just launched some more simulations with dr * 1e-10 to try and diagnose.

ID Notes Issues Status $N$ movie
264 base None benson, hobson, alice (stopped) 12
265 $X_1 = 0.1 X_1^{equ}$ worked with upwind convection term renfield, sancho, kato (stopped) 12
266 $\nu_f = 2 \times 10^{-2}$ crashed luggage (stopped) 12
  • Posted: 2009-12-17 12:32 (Updated: 2009-12-22 11:31)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Useful Simulations

ID Notes diff Issues Status $N$ movie
222 base diff benson, hobson, alice (stopped) 12
225 $\nu_f = 2 \times 10^{-2}$ diff luggage (stopped) 12
226 $\nu_f = 2 \times 10^{-1}$ diff crashed renfield, sancho, kato (stopped) 12
233 222: $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff cluster 7.0 (stopped) 4
240 222: $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed luggage (stopped) 12
252 243: $\nu_f = 2 \times 10^{-1}$, shift=dx, $\sigma_{\text{CFL}} = 0.01$ diff cluster 128 (stopped) 4
257 251: $X_1 = 1.0 X_1^{equ}$, $\bar{M}_f = 1 \times 10^{-6}$, shift=2*dx, $\sigma_{\text{CFL}}=0.05$ diff luggage (stopped) 12
  • Posted: 2009-12-15 15:32 (Updated: 2009-12-17 11:59)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Edwin Schwalbach Model

For the last three days Eddie and I have been working on his VLS growth model. We've been trying to implement it in FiPy?. We have successfully implemented a two phase, pure, 1D system. His notes are here

source:trunk/vls/Single_component_perturbation.pdf@1031

The code is here

source:trunk/vls@1031

Essentially, this is simply the van der Waals model with a free energy of the form

$ F_m = w \left( \rho - \rho_l \right)^2 \left( \rho - \rho_v \right)^2 $

We implement the model using the block matrix approach with $u$, $\rho$ and $\mu^{NC}$ as the base variables to "solve for". The $\mu^{NC}$ is linearized with,

$ \mu^{NC} = \mu_0^* + \rho \frac{\partial \mu}{\partial \rho}^* - \kappa \nabla^2 \rho $

where,

$\mu_0^* = \mu^* - \rho^* \frac{\partial \mu}{\partial \rho}^* $

Remember,

$F_m = \frac{m F_v}{\rho}$

and

$P = -F_v + \rho \mu$

and

$\frac{ \partial \mu}{\partial \rho} = \frac{1}{\rho} \frac{ \partial P}{\partial \rho}$

Explicitly,

$P = \frac{2 W \rho^2}{m} \left( \rho - \rho_l \right) \left( \rho - \rho_v \right) \left( 2 \rho - \rho_l - \rho_v \right) $

and

$\frac{\partial P}{\partial \rho} = \frac{2 W \rho}{m} \left( \rho \left( 2 \rho - \rho_l - \rho_v \right)^2 + 2 \left(3 \rho - \rho_l - \rho_v \right) \left( \rho - \rho_l \right) \left( \rho - \rho_v \right) \right) $

Source Terms

Take an equation say,

$ \frac{\partial \phi}{\partial t} = S \left( \phi \right) $

In general $ S \left( \phi \right) $ can be implemented explicitly (i.e. $S \left( \phi^* \right) $ where $*$ is the previous sweep value. Now, to implement this implicitly, we need to linearize,

$ S \left( \phi \right) = S \left( \phi^* \right) - \phi^* \frac{\partial S}{\partial \phi}^* + \phi \frac{\partial S}{\partial \phi}^* $

or

$S = S_c + \phi S_p$

where

$S_p = \frac{\partial S}{\partial \phi}^*$

In general $S_p \le 0$ for stability reasons, though one can often relax this somewhat if the source is not dominant. If the source is not negative, then $S$ should be evaluated entirely explicitly. The reason $S_p$ must be negative is the following. Take,

$ \frac{\partial \phi}{\partial t} = \lambda \phi $

then the discretized implicit form is

$ \phi = \frac{\phi^O}{1 - \lambda \Delta t} $

Say $\phi^O$ is positive, then this equation should increase the value of $\phi$, but if $ \lambda \Delta t > 1$ then we get a negative value of $\phi$. If $\lambda$ is negative, there are no stability criterion. Conversely, if $\lambda<0$ then the explicit scheme,

$ \phi = \left( 1 + \lambda \Delta t \right) \phi^O $

will drive the solution below $0$ for sufficiently negative $\lambda$. For very gentle intro,

clcik here

Pressure velocity coupling

We use the Rhie-Chow interpolation method to deal with this issue since $P$ and $u$ are not on staggered grids but collocated. Essentially, the Rhie-Chow interpolation scheme corrects the face velocities in the continuity equation by the use of a corrected velocity that approximates solving the momentum equation on the faces. Let's take the discretized form of the x-momentum equation:

$ a_P u_P = \sum_{nb} a_{nb} u_{nb} - V_P \rho \nabla_x \mu^{NC} $

where $P$ is a cell value. For small enough grid, we would like

$ a_f u_f - \sum_{nb} a_{nbf} u_{nbf} + V_f \rho_f \nabla_x \mu^{NC} = \overline{a_P u_P} - \overline{\sum_{nb} a_{nb} u_{nb}} + \overline{V_P \rho \nabla_x \mu^{NC}} $

that is to say that we would like the residual of the momentum equation evaluated at the face to be equal to the residual of the averaged residual over the neighboring cells. We assume that the $\sum_{nb}$ terms are unimportant and say that,

$ u_f = \frac{\overline{a_P u_P}}{a_f} - \frac{V_f}{a_f} \rho_f \nabla_x \mu^{NC} + \frac{\overline{V_P}}{a_f} \overline{ \rho \nabla_x \mu^{NC}} $

We approximate $a_f = \overline{a_P} $ and $V_f = \overline{V_P}$ the we get,

$ u_f = \overline{u_P} - \frac{\overline{V_P}} {\overline{a_P}} \left( \rho_f \nabla_x \mu^{NC} - \overline{ \rho \nabla_x \mu^{NC}} \right) $

The continuity equation can be written,

$ \frac{\partial \rho}{\partial t} + \nabla \cdot \left( u_f \rho \right) = 0 $

Substituting in gives a diffusion term for $\mu^{NC}$, which acts to stabilize the solution. See source:trunk/vls/input.py@1031#L78.

Parsitic Currents

Discretize using $\rho \nabla \mu$ rather than $\nabla P$. This induces dissipation assoicated with momentum. See source:trunk/reactiveWetting/paper/jacqmin.pdf@1031.

Movie of solution to equations

  • Posted: 2009-12-11 13:40 (Updated: 2009-12-11 16:04)
  • Author: wd15
  • Categories: (none)
  • Comments (0)