Posts in category reactiveWetting

Reactive Wetting Equations in Cylindrical Coordinates

$\partial_j = \left( \partial_r, \partial_z \right) $

Note

We need to remove the radial distance $r$ and embed it in the term's coefficient. For a general equation in cylindrical coordinates that are axisymmetric the general equation can be written,

$ \frac{\partial \phi}{\partial t} + \frac{1}{r} \partial_r \left( r u_r \phi \right) + \partial_z \left( u_z \phi \right) = \frac{1}{r} \partial_r \left( r \Gamma \partial_r \phi \right) + \partial_z \left( \Gamma \partial_z \phi \right)$

Multiply through by $r$ to get,

$ \frac{\partial \left( r \phi \right)}{\partial t} + \partial_r \left( r u_r \phi \right) + \partial_z \left( r u_z \phi \right) = \partial_r \left( r \Gamma \partial_r \phi \right) + \partial_z \left( r \Gamma \partial_z \phi \right)$

which is,

$ \frac{\partial \left( r \phi \right) }{\partial t} + \partial_j \left( r u_j \phi \right) = \partial_j \left( r \Gamma \partial_j \phi \right) $

Thus, in general, it is very easy to recast equations in cylindrical axi-symmetric coordinates with fipy. The fourth order term in the concentration equations may be an issue and the correction equation also has a fourth order term.

Fourth Order Term

$\frac{1}{r} \partial_r \left( r \Gamma_1 \partial_r \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \partial_z \left( \Gamma_1 \partial_z \phi \right) \right] \right) + \partial_z \left( \Gamma_1 \partial_z \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \partial_z \left( \Gamma_1 \partial_z \phi \right) \right] \right) $

This can be rewritten as

$\frac{1}{r} \partial_r \left( r \Gamma_1 \partial_r \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \frac{1}{r} \partial_z \left( r \Gamma_1 \partial_z \phi \right) \right] \right) + \frac{1}{r} \partial_z \left( r \Gamma_1 \partial_z \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \frac{1}{r} \partial_z \left( r \Gamma_1 \partial_z \phi \right) \right] \right) $

There is an awkward $\frac{1}{r}$ that is not next to a $\Gamma$ and thus can not just multiply a regular coefficient. The outer $\frac{1}{r}$ gets dealt with because we multiply all the terms in the equation by $r$. The above can be rewritten as,

$\frac{1}{r} \partial_j \left( \Gamma_1 \left[\partial_j - \frac{\hat{r}_j}{r} \right] \partial_k \left[ r \Gamma_2 \partial_k \phi \right] \right)$

Unfortunately, The last term needs to be added as a source term. If we had a third order convection term then it could be added implicitly.

Continuity

$ \frac{\partial \left(r \rho \right)}{\partial t} + \partial_j \left( r u_j \rho \right) = 0 $

Velocity

$ \frac{\partial \left(r \rho u_i \right)}{\partial t} + \partial_j \left( r u_j u_i \rho \right) = \partial_j \left( r \mu \left[ \partial_j u_i + \partial_i u_j \right] \right) - r \rho_k \partial_i \mu_k + \epsilon_k T r \rho_k \partial_i \left( \frac{1}{r} \partial_j \left(r \partial_j \rho_k \right) \right) $

Concentration Equation

$ \frac{\partial \left(r \rho X_2 \right)}{\partial t} + \partial_j \left( r u_j X_2 \rho \right) = \partial_j \left( \frac{ r \bar{M}}{T} X_1 X_2 \left[ \left(\frac{ \left( A_2 - A_1 \right) }{m} - \frac{ \left( e_2 - e_1 \right) }{m^2} \rho \right) \partial_j p + \left( 1 - p \right) \frac{ \left( e_2 - e_1 \right)}{m^2} \partial_j \rho \right] \right) $

$ + \partial_j \left( \frac{ r \bar{M} R }{m} \partial_j X_2 \right) - \partial_j \left( \bar{M} X_1 X_2 \left[ \epsilon_1 + \epsilon_2 \right] \left[\partial_j - \frac{\hat{r}_j}{r} \right] \partial_k \left[ r \rho \partial_k X_2 \right] \right) $

$+ \partial_j \left( r \bar{M} X_1 X_2 \partial_j \left[ \frac{1}{r} \partial_k \left( r \left[ \epsilon_1 X_1 - \epsilon_2 X_2 \right] \partial_k \rho \right) \right] \right)$

The last term can remain in its original form as it is a convection term.

Alternate Method

The formulation above could lead to some coding issues and maybe other difficulties concerning the fourth order source terms. A more expedient and relatively trivial method that is probably no less accurate is to define a CylindricalGrid2D object that has its mesh volumes and face areas modified accordingly. I just realized that the discretiztion is exactly equivalent in the finite volume method when using a cylindrically modified grid and rewriting the equations.

Reactive Wetting Paper

Outline

  1. Introduction
    1. Boiler plate - motivation
    2. Summary of physical parameters coinciding with typical applications (water-ice-salt, solder alloys)
    3. Thermodynamics
  2. Model
    1. Transport
    2. Thermodynamics
    3. Surface Tensions
  3. Parameters Choices - Numerics
    1. Governing equations summary.
    2. Numerical Approach
      1. Solution algorithm
      2. Parasitic Currents
    3. One dimensional liquid vapor analysis
      1. Dimensionless form of the pure-liquid vapor
      2. Convergence of a one-dimensional pure liquid-vapor system
        1. $\delta$,
        2. $Re$
        3. $M$,
        4. $Ca$
        5. $\kappa$ liquid
      3. Convergence of a one-dimensional binary liquid vapor system
        1. $Pe$, $Pe \rightarrow \infty$, What is the meaning?
        2. $Re > 1$ issues.
    4. Two dimensional Binary Solid-Liquid-Vapor
      1. Obtaining Plausible contact angles
      2. $\kappa$ solid
      3. Approximation of a solid with $\mu \rightarrow \infty$, density trapping.
      4. Comparison with analytical solution in equilibrium
  4. Results
    1. Extracting interface locations and junctions with phase field methods
      1. Ridge detection
      2. Wheeler method
      3. Villanueva method
    2. $\theta_{\text{contact}}$ versus $t$.
    3. Cox comparisons - look at later cox work with surfactants.
    4. $R\left(t\right)$ for different initial conditions.
    5. Pure Hydrodynamics with and without reactive wetting ($\theta_{\text{liquid}}$ vs. $\vec{v}$)
    6. Pretty pictures
    7. Compare with Gustav
      1. Steady state versus transient
      2. Compare streamlines
    8. $\theta$ versus diffusion
      1. $\bar{M}\ll 1$ versus $\bar{M}=0$
      2. $\bar{M}_{\phi}=0$
    9. Compare with S. Troian
    10. Spherical cap deviation
    11. Best measures of angle and shape
    12. Compare with Jim and Bill's right angle assumption

Copies of notes

Warren's notes

Wheeler's notes

openmp/weave timings.

A matrix multiplication in weave really scales well with openmp. The code is here. The observed speedup is almost perfect with two threads.

This code, involving large array multiplications of size N, has the following speedups with two threads.

N Speedup
1E7 1.51
1E6 1.37
1E5 1.39
1E4 1.0

It should be noted that the number of loops in python increased inversely with the size of the array.

The question remains whether we can get speed ups for smaller arrays typically used in FiPy.

Streamlines

How to make a contour plot of the streamlines for a flow field in 2D?

Given a velocity field $\vec{u}$ we need a field $\psi$ that is constant in directions that are tangent to the velocity field. This leads to the equation,

$$ \vec{u} \cdot \nabla \psi = 0$$
The relationships
$$ u = \frac{ \partial \psi }{ \partial y} $$
and
$$v = -\frac{ \partial \psi }{ \partial x} $$
satisfy the equation. From the above equations we can say,
$$ \nabla^2 \psi = \frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} $$
Solving this gives the stream-function. If we have no inflow conditions on regular domain walls, then we can choose $\psi=0$ on these walls. Since these walls have no circulation, we have made the choice that $\psi=0$ represents no circulation, which is sensible. The circulation orientation is given by the sign of $\psi$.

See an example of streamlines

FiPy code to compute streamlines

2D simulations integrated with subversion

ID status condor ID $\Delta t_{\text{max}}$ $R$ $\sigma_{\text{CFL}}$ $\bar{M}_S$ $\bar{M}_L$ $\mu_S$ $\mu_L$ Jim's movie Bill's movie
1 unstable $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
2 running 12095.0 $1\times 10^{-4}$ $5\times10^{-7}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
3 stopped $1\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
4 running 12096.0 $1\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
5 unstable $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
6 running 12097.0 $1\times10^{-6}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
7 stopped $1\times 10^{-4}$ $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ movie movie
8 running 12100.0 $1\times 10^{-4}$ $1\times10^{-6}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ movie movie
9 unstable $1\times 10^{-4}$ $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
10 running poole 5545 $1\times 10^{-4}$ $2\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
11 running poole 9411 $1\times 10^{-4}$ $4\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
12 stopped $1\times 10^{-4}$ $2\times10^{-6}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
13 stopped $1\times 10^{-4}$ $4\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie

2D simulations

Table of simulations

ID status condor ID Initial T $\Delta t_{\text{max}}$ $R$ $\sigma_{\text{CFL}}$ $\bar{M}_S$ $\bar{M}_L$ $\mu_S$ $\mu_L$ data directory movie link source files
H2D-restart H2D 550 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ movie
I2D-restart I2D 550 $1\times 10^{-8}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ movie
J2D-restart running 9099.0 J2D 600 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Feb-29-12:25:35-2008 movie files
K2D-restart running 9098.0 K2D 600 $1\times 10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Feb-29-12:25:35-2008 movie files
L2D unstable 9077.0 650 $5 \times 10^{-7}$ $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-27-16:59:24-2008 movie files
L2D1 unstable 9087.0 650 $5 \times 10^{-7}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-28-11:27:18-2008 movie files
M2D running 9088.0 650 $1 \times 10^{-6}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-28-11:28:02-2008 movie files
N2D swapping 9089.0 650 $2 \times 10^{-6}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-28-11:31:03-2008 movie files
O2D unstable 9090.0 650 $5 \times 10^{-7}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataFri-Feb-29-01:27:24-2008 movie files
O2D1 running 9108.0 650 $5 \times 10^{-7}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataFri-Feb-29-14:35:04-2008 movie files
P2D running 9111.0 650 $1 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataSat-Mar--1-07:01:06-2008 movie files
Q2D swapping 9112.0 650 $2 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataSun-Mar--2-22:30:18-2008 movie files
R2D unstable 9083.0 650 $5 \times 10^{-7}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-07:43:56-2008 movie files
R2D1 running 9093.0 650 $5 \times 10^{-7}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-11:37:49-2008 movie files
S2D running 9094.0 650 $1 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-11:38:37-2008 movie files
T2D swapping 9095.0 650 $2 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-11:41:06-2008 movie files

Upper bound on CFL condition

The simulations below do not seem to have a bounded $\Delta t$. They do, however, have an upper bound for $\sigma_{\text{CFL}}$ that decreases dramatically as $\mu_L$ is reduced. Of course, at lower $\mu_L$ the upper bound for $\Delta t$ is probably due to some other numerical factor unrelated to $\sigma_{\text{CFL}}$.

Matrix of simulations

ID status $\Delta t_{\text{max}}$ $L$ $f_V$ $f_L$ $f_S$ $\sigma_{\text{CFL}}$ $\bar{M}_S$ $\bar{M}_L$ $\mu_S$ $\mu_L$ data directory movie link
N finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-14:55:59-2008 movie
N3 unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $2\times 10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-27-12:27:04-2008 movie
O finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataMon-Feb-25-17:58:15-2008 movie
P unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:30:26-2008 movie
Q finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:57:53-2008 movie
Q1 finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataWed-Feb-27-14:51:24-2008 movie
R unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-10:13:19-2008 movie
S unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-3}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-11:45:36-2008 movie

2D simulations with reduced temperature

Matrix of 2D simulations with reduced temperature

ID status condor ID Initial T $\Delta t_{\text{max}}$ $R$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
A2D finished 6952.0 650 $1\times 10^{-7}$ $3.33 \times10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-16-00:48:07-2008 movie
B2D unstable 7672.0 A2D 550 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Jan-25-18:17:34-2008 movie
C2D finished 7828.0 650 $1\times 10^{-7}$ $6.66 \times10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-11:53:40-2008 movie
D2D finished 7836.0 650 $1\times 10^{-7}$ $1.33 \times10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-12:19:30-2008 movie
E2D finished 7848.0 650 $1\times 10^{-7}$ $1 \times10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-13:52:02-2008 movie
F2D finished 7859.0 650 $1\times 10^{-7}$ $5 \times10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-14:09:54-2008 movie
G2D finished 7869.0 650 $1\times 10^{-7}$ $2 \times10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-14:15:37-2008 movie
H2D finished 8541.0 F2D 550 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataThu-Feb-14-12:09:33-2008 movie
I2D finished 8582.0 F2D 550 $1\times 10^{-8}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Feb-15-14:55:08-2008 movie
J2D finished 8803.0 F2D 600 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataTue-Feb-19-10:20:01-2008 movie
K2D finished 8805.0 F2D 600 $1\times 10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataTue-Feb-19-10:20:56-2008 movie

Most of these jobs ran for a good period of time, but were brought down with server problems on 02/14/08 and 02/22/08. When they restarted the condor system began them from scratch, I killed them and will restart any useful ones from the latest point in the data directory above. However, it is probably best to keep processors for the new set of 2D jobs as these will probably have unbounded time steps.

Efficiency

It is interesting to compare time versus time steps for some of the 2D jobs and the new 1D jobs. $N$ represents number of time steps.

Simulations E2D, F2D and G2D all lie on the same curve as expected. Elapsed time in simulations N and O is orders of magnitude greater than in the 2D simulations for an equivalent number of time steps. This is not the case in simulation O as it has a more restrictive CFL condition (see Table).

1D simulation with vapor

Previously, we've managed to get the 1D liquid-solid system working with all the desired numbers (apart from the liquid viscosity) with only a CFL time step stability criterion. We now need to include the vapor and see how the time step is affected. I'm assuming from previous experience and intuition that adding the vapor will require a far stricter stability limit.

Simulation With Vapor

As per usual my intuition proves to be false. The time step size is completely unrestricted. I'm extremely excited. We may be able to do the 2D simulations in orders of magnitude less real time. There are also some extra simulations to look at initial liquid fraction and domain size. It is clear from simulation N that diffusion based interface motion is occurring at these parameter values.

Diffusion Driven Interface Motion

Examining simulations N, O, Q, it is evident that interface motion is driven by diffusion. Notice also that the pressure equilibrates sooner as the viscosity is reduced and it seems that the concentration gradients are more pronounced. In simulation Q, the pressure is equilibrated much earlier, before significant interface motion occurs and concentration gradients form.

Stability

To maintain stability for lower $\mu_l$, $\sigma_{\text{CFL}}$ is being reduced. Simulations P, R, S were all unstable. The lowest value of $\mu_l$ that was successfully computed is $10^2$, simulation Q. The following figure shows simulation time against number of required steps ($N$). It is clear that as $\sigma_{\text{CFL}}$ the run times are significantly increased.

Interface Motion

Interface motion plot for simulation N, O and Q. At least between $\mu_l=10^4$ and $\mu_l=10^2$ the solid liquid interface motion is mostly unaffected by $\mu_l$.

Matrix of simulations

ID status $\Delta t_{\text{max}}$ $L$ $f_V$ $f_L$ $f_S$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
H finished None $10^{-6}$ 0 1/2 1/2 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:46:23-2008 movie
L finished None $10^{-6}$ 1/3 1/3 1/3 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-11:44:28-2008 movie
M finished None $10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-12:08:13-2008 movie
N finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-14:55:59-2008 movie
O finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataMon-Feb-25-17:58:15-2008 movie
P unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:30:26-2008 movie
Q running None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:57:53-2008 movie
R unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-10:13:19-2008 movie
S unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-3}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-11:45:36-2008 movie

Realistic Solid Behavior

The following plot includes the curves for each of the 1D simulations (see Table).

It is clear that $\bar{M}_s$ has more influence on the interface motion than $\mu_s$ during melting.

The next step is to try simulation H with the vapor added in 1D to figure out a satisfactory time step.

Also, it would be interesting to know how $\mu_l$ influences the interface motion.

Finite Size Effects

Having examined simulations I and K and compared them with B and D (see Table), it seems that there may be finite size effects in these simulations. The interface motions are way off. For example simulation B and I use the same parameter set yet the interface motion is over two orders of magnitude slower for case I. Some simulations need to be performed to ascertain whether there is an eventual convergence of the interface motion as the domain is expanded.

Table of Simulations

ID status $\Delta t_{\text{max}}$ $L$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
B finished None $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^4$ $10^4$ dataWed-Feb-20-14:09:53-2008 movie
I finished None $10^{-5}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{4}$ $10^4$ dataThu-Feb-21-18:08:24-2008 movie
J finished None $2\times10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{4}$ $10^4$ dataFri-Feb-22-15:14:20-2008 movie

Chat with Jim

After a chat with Jim we both came to the realization that curves on a liquid fraction against time plot for varying domain sizes do not have to correspond at all. That raises the question about how to determine whether we have finite size effects. There does not seem to be an easy way to do this. Certainly plotting interface position as a function of time may be more enlightening. The velocity of the interface at early times should be independent, but at long times I don't think that is true.

Absolute Interface Position

It is clear from the figure that all our fears are unfounded and this whole episode was a waste of time. Anyhow, it is interesting to notice that there is more "ringing" in the larger domain.

High solid viscosity

Jim wants to understand why we cannot model a rigid solid interface with a high viscosity. He needs the story for his talk. Suggestion is to run a simulation with a very high viscosity in the solid and a liquid that is not at equilibrium to induce melting of the solid

Initial Condition

The initial condition is a phase fraction of 1/2 liquid and 1/2 solid. The solid is at the three phase equilibrium values, $\rho_1^s$ and $\rho_2^s$. The initial densities in the liquid region are given by,

$ \rho_1 = \rho_1^l + \lambda \left( \rho_1^s - \rho_1^l \right) $

$ \rho_2 = \rho_2^l + \lambda \left( \rho_2^s - \rho_2^l \right) $

In this case $\lambda=-0.1$. This induces melting of the solid. For the numbers that we use in the example the final liquid fraction is given by 0.55. Using this approach there is a lower bound for $\lambda$ and thus an upper bound on the final liquid fraction.

Regular solid viscosity

When $\mu_s=1000$ the solid behaves much like a liquid. The movie is given here for $\Deltat=1\times10^{-6}$ (dir: dataWed-Feb-20-10:37:23-2008). In this simulation there is no clear distinction between interface motion due to pressure equilibration and motion due to the concentration gradient. In simulations where the temperature is reduced there is a clear distinction between these regimes of interface motion. One issue is that the total interface motion is equivalent to the interface width. It is definitely a good idea to run simulations with a small interface width compared with the total interface motion.

Interestingly enough, we can run this simulation with no upper bound on time step other than CFL. The parasitic velocities eventually give an upper bound less than 1. Here is the movie (dir: dataWed-Feb-20-14:09:53-2008)

High solid viscosity

When $\mu_s=1\times10^{12}$ the simulation shows similar interface motion. Here is the movie (dir: dataWed-Feb-20-14:36:18-2008). A simulation has also been run with a bounded $\Delta t$ (dir: dataWed-Feb-20-10:47:36-2008)

Liquid Fraction plot

The following plot is the liquid fraction against time for the whole matrix of simulations. The simulation with $\Delta t < \Delta t_{\text{CFL}}$, $\bar{M}_s=10^{12}$ and $\Delta t<\Delta t_{\text{CFL}}$ shows rather strange long time behavior. $\Delta t\approx5$ at long times, which could be too large.

Low diffusion in the solid

To get closer to real solid behavior, the value of $\bar{M}_s$ can be reduced. In this movie, $\bar{M}_s=10^{-15}$ (dir: dataWed-Feb-20-18:22:20-2008).

Matrix of simulations

ID status $\Delta t_{\text{max}}$ $L$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
A finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^4$ $10^4$ dataWed-Feb-20-10:37:23-2008 movie
B finished None $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^4$ $10^4$ dataWed-Feb-20-14:09:53-2008 movie
C finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-20-10:47:36-2008 movie
D finished None $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-20-14:36:18-2008 movie
E finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-15}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-20-18:22:20-2008 movie
F finished None $10^{-6}$ 0.1 $10^{-15}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:16:51-2008 long short
G finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:40:48-2008 movie
H finished None $10^{-6}$ 0.1 $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:46:23-2008 movie
I finished None $10^{-5}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{4}$ $10^4$ dataThu-Feb-21-18:08:24-2008 movie
K finished None $10^{-5}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-18:11:13-2008 movie

2008/02/20/16.43

Test