-
Notifications
You must be signed in to change notification settings - Fork 145
reactiveWetting
Completed a simulation with the second set of my own parameters. The X2 formulation was used. The movie is attached. The initial fractions are 1/3 for solid, liquid and vapor, respectively. The initial densities are
The equilibrium densities are
and the final fractions are 1/6, 1/2 and 1/3 for solid, liquid and vapor, respectively. The final simulation equilibrium densities were,
In the [attachment:wiki:reactiveWetting:X2-equilibrium-displaced.mov] the solid-liquid interface moves right across the domain, leaving only a small amount of liquid. It then takes a large number of steps for the solid-liquid interface to move back to the correct position. In the movie the solid
liquid interface still seems to be not at equilibrium (the solid fraction is ~0.2), although the densities are pretty much correct. This is most probably due to a small energy barrier solid and liquid phases. The simulation took about 70000 sweeps to get to this stage. This is probably too much for a 2D run.
I have also tried running Bill's numbers using the updated X2 formulation (i.e with the bug fix). I still get the same issues. The velocity equation has the large residual, which decreases, but starts out so big that it takes a large number of sweeps to get to the required tolerance. The problem is that after a good number of steps, the density develops oscillations. Is this due to not having a high enough tolerance? I don't know. The equations do decrease in tolerance so it is rather confusing.
Note: the above simulation was run with the wrong relaxation factors. They were set to 0.2 and 0.5 rather than 0.5 and 0.2 for the velocity and density relaxation factors respectively. Not sure if this will make much difference, may slow convergence.
* The above material is one where the density difference is quite large but the fractions are similar for all phases. * Probably good idea to create a material were the fractions have large differences but the densities are fairly similar across all three phase. This may not be easy. * Next job is to evaluate the surface tensions for all three interfaces.
The solid liquid problem converges very quickly. 1866 sweeps and the velocities are less than 1e-7. The equilibrium densities should be
they turn out to be
That will do. (Data dir: dataX2-solid-liquidThu-Jul-26-17:32:09-2007)
Note: Jon just had a good idea. I should plot the rho1, rho2 plot through the interface to figure out the path. This may shed some light on why the solid-vapor interface always seems to form a liquid.
Numbers for the liquid-vapor interface should be:
The turn out to be,
The data directory is dataX2-liquid-vaporThu-Jul-26-17:58:47-2007 For solid-vapor they shoud be
They turn out to be
The data directory is dataX2-solid-vaporThu-Jul-26-17:35:45-2007 From the simulation it is clear that the solid-vapor interface has a liquid interface.
The surface tensions are as follows using the two methods
* liquid-vapor * method1: 0.0555091295871 * method2: 0.0549549862314 * solid-vapor * method1: 0.0860515746592 * method2: 0.0854607675124 * solid-liquid * method1: 0.0305386895437 * method2: 0.0305067688696
As anticipated, the solid-liquid and liquid-vapor sum is about equal to the solid-vapor surface tensions for both methods. The density space plot Image(densitySpace.png, 60px) demonstrates that the solid-vapor path goes close to the liquid equilibrium point. This gives an indication that if we can separate the solid and liquid equilibrium points there may be less chance of the path going via the liquid point.
Today I am going to try to create a material that has a greater distance between the solid and liquid equilibrium points. The plan is to vary the A's and leave the e's alone for now. I found a new material. It has the following equilibrium conditions,
I am running the equilibrium condition v-l-s just to check it, then I'll run the 3 surface tension calculations and plot the path for this material. Frome the simulation the equilibrium densities are
when velocity is below 1e-7. The solid-liquid equilibrium values are:
aand solid vapor are
and liquid vapor are
Again, even with this material, in which the solid and liquid are further apart on the phase diagram, the solid-vapor surface tension is the addition of solid-liquid and liquid-vapor surface tensions. The plot of the paths on the phase diagram [attachment:wiki:reactiveWetting:phaseSpace.png] shows that the vapor-liquid path goes close to the liquid.
Jim's idea is to plot the excess free energy (f - f_inf - mu * (c - c_inf)). That will show where the humps and saddles are in phase space. The plan is to parametrize along each of the three phase region lines with parameter t. The contour plot will then have t and phase as the axes. The plots along the [attachment:wiki:reactiveWetting:vapor-liquid.png] line, the [attachment:wiki:reactiveWetting:vapor-solid.png] line and the [attachment:wiki:reactiveWetting:liquid-solid.png] line show the excess energies for various values of phase and parameterized density.
The linked plots above are line plots, here are the corresponding contour plots for [attachment:wiki:reactiveWetting:vapor-liquid-contour.png], [attachment:wiki:reactiveWetting:vapor-solid-contour.png] and [attachment:wiki:reactiveWetting:liquid-solid-contour.png].
Jon pointed out that the [attachment:wiki:reactiveWetting:vapor-solid-contour.png] has two minimums. One of the minimums lies along phi=0. Thus, it it seems better for the path to stay at phi = 0 while the density varies from the vapor to the liquid density and then change phi rapidly across the small change in density between the liquid and solid equilibrium densities.
The basic idea is to create a scatter plot of minimum excess energy in the three phase region. Firstly, I create parameterized lines between solid and liquid, solid and vapor and solid and any point on the liquid-vapor line. For each line I find the poisition of the minimum for various values of the phase field. This will give (rho1, rho2, f_excess_min) for a whole bunch of scattered points.
Jim noticed that the [attachment:wiki:reactiveWetting:vapor-solid-contour.png] has a jump and thus the path cannot go along the vapor solid line. He also noticed that the free energy around the solid point is extremely steep. The solution to this may be to increase W while maintaining the interface thickness
}}} . This should make the jump between the liquid and the solid higher. From the [attachment:wiki:reactiveWetting:liquid-solid-contour.png], the barrier is fairly low.
When,
the contour plot change so that there is a greater barrier between ths solid and the liquid. Image(liquid-solid-contour-W-1000.0.png, 60px) Image(vapor-solid-contour-W-1000.0.png, 60px) Image(vapor-liquid-contour-W-1000.0.png, 60px) Remembering tha the equilibrium conditions for this system are
The three phase equilibrium is
The veloicty was less than 1e-7. The solid-vapor equilibrium conditions for nx=200 are
It is clear from the denisty plots that an intermediate phase has formed. (dir: dataX2-params3-vapor-solidWed-Aug--1-17:52:14-2007). If the resolution is increased to nx=400. The equilibrium values become,
The paths on the phase diagram for both nx=200 and nx=400 correspond. Image(phaseSpace-W10000.png, 60px)
New system with a gentle free energy has equilibirum values of
For the three phase equilibrium problem the simulation gives
The solid-vapor equilibrium values are
The solid-liquid equilibrium values are
The liquid-vapor equilibrium values are
The paths on the phase diagram don't look much better Image(phaseSpaceGentle.png, 60px) This material has not improved matters.
The following material has a very gentle solid free energy and the solid and liquid equilibrium points are well spaced. The parameters for this material are:
The equilibrium values for the three phase simulation are:
The equilibrium values for the solid-vapor simulation are:
The equilibrium values for the solid-liquid case are
The equilibrium values for the liquid-vapor case are
Awesome! From this, we are looking at a contact angle of ~73 degrees. It is clear from the phase diagram that the solid-vapor path goes no where near the liquid equilibrium point. Image(goodPhaseSpace.png, 60px)
I have tried to run a number of 2D jobs, but the density is increasing over 3.0 in the solid. This is resulting in the following error
Looking at the movie for data set "dataX2-params8-2D-largeWed-Aug--8-11:21:46-2007", the effect seems to be 1D as a wave of density is caused in the solid away from the liquid. Why does this not happen in the 1D calculations? The 1D calculations reach 2.88 for solid liquid case and 2.83 for the three phase case. The sampling for these two cases is only every 10 steps, so it may have got higher. Also, the 2D case could have had a wave that concentrated in a corner and thus increased above 3. For now I've corrected all cases of "molecularWeight - density * vs" to stay above 1e-20. I've launched a new 2D case with this condition. Still get the same problem even when I prevent the log argument becoming negative, which makes no sense. Running a job with a higher viscosity may help damp the density waves in the solid.
|| viscosity || max density || || 1 || 2.81 || || 10 || 2.81 || || 100 || 2.71 ||
Increasing the viscosity in the solid seems to keep the maximum density under control. Let's get it on! BTW Higher viscosity doesn't seem to stiffen things up very much. After about 4300 the max velocity is under 1e-5.
2D simulation is looking great. I need to get the following done:
* movie with velocity vectors * add a write up about the analytical solution and compare with simulation * measure the solid, liquid and vapor proportions at the end and the beginning of the simulations * How does this material behave with temperature? Bill's question.
Bill thinks that we need to figure out what actually separated the trajectories on the phase diagram. Bill also wants to know how the solid behaves with temperature. His material became soft as the temperature was reduced.
Create some tables for dimensionless parameters for the system I am using and for some of the previous systems, including Bill's old system and new system.
The 2D movie is finally ready [attachment:wiki:reactiveWetting:param8-2D.mov].
This model consists of two caps of constant curvature. The solid-vapor interface is assumed to be flat. The law of signs gives,
$ \cos{\theta_1} = \frac{ \left( \gamma_{sv}^2 + \gamma_{lv}^2 - \gamma_{sl}^2 \right) }{ 2 \gamma_{sv} \gamma_{lv} } $
$ \cos{\theta_2} = \frac{ \left( \gamma_{sv}^2 - \gamma_{lv}^2 + \gamma_{sl}^2 \right) }{ 2 \gamma_{sv} \gamma_{sl} } $
Define,
$ A_i = \frac{ \left( \theta_i - \cos{\theta_i} \sin{\theta_i} \right) }{ 2 \sin{\theta_i} } $
$ A = A_1 + A_2 $
and
$ a_i = A \rho_i^l - A_i \rho_i^s $
$ b_i = \rho_i^s W $
where W is the width of the box. The solid height is given by,
$ h = \frac{ Q_1 a_2- Q_2 a_1 } { b_1 a_2 - a_1 b_2 } $
and the radius of the drop is given by,
$ R = \sqrt{ \frac{ Q_1 b_2 - Q_2 b_1 } { a_1 b_2 - b_1 a_2 } } $
where Q_1 and Q_2 are the initial mass of each component.
For the numbers used in the simulation the analytical model is prediciting a height for the solid of 0.731 and a drop radius of 0.801. The simulation has not yet reached equilibrium and is giving a number of 0.755 using the 0.5 phase field position at the solid-vapor interface.
The data set is "dataX2-params8-2D-largeFri-Aug-10-16:02:11-2007".
From the simulation, the radius of the drop is about 0.7. This differs quite a lot from the expected analytical result of 0.8.
From the simulation the equilibrium values are
(0.0250175, 0.0665361, 2.0668057, 0.3346652, 0.3024389, 2.3836282)
and from the phase diagram.
(0.0237341, 0.0665500, 2.0616042, 0.3366000, 0.2879295, 2.3999297)
A choice for dimensionless numbers: \tilde{F} }}} Using this scale change it is clear that the equilibrium is determined from , \;\;\; \frac{B}{R T} }}} From the free energy there is one other parameter that is concerned with the shape of the interface, }
I added an updated [attachment:wiki:reactiveWetting:Untitled.mov] showing the fraction contour plot.
Anyway, back to parameter analysis. Dimensionless chemical potentials, , \;\;\; \mu_2^{\text{NC}} = \frac{R T}{m} \tilde{\mu}_2^{\text{NC}}, \;\;\; }}} The diffusion equation can then be written, {M_{\phi} \epsilon_{\phi} m^2} \tilde{\partial}_j \left( X_1 X_2 \tilde{\partial}_j \left[ \tilde{\mu}_2^{\text{NC}}] \right) }}} where the time scale is given by X / V and
This introduces two more dimensionless parameters,
}}} The phase field equation becomes,
+ \tilde{v}_j \tilde{\partial}_j \phi = \tilde{\partial}_j^2 \phi - \frac{\partial \tilde{F}}{\partial \phi} }}} The flow equation becomes, {M_{\phi} \epsilon_{\phi} m} \tilde{\partial}_j \left( \tilde{\partial}_j \tilde{v}_i + \tilde{\partial}_i \tilde{v}_j \right) - \frac{T \bar{v}}{m M_{\phi}^2 \epsilon_{\phi}} \left[\tilde{\rho}_1] }}}
From above there are 6 parameters concerned with the classical equilibrium and a further 4 parameters concerned with the non-classical equilibrium. There are also three parameters concerned with dynamics.
We have a new list of tasks in order of importance:
* Reduce Mbar by up to 5 orders of magnitude * Revisualize the rho2 / rho plot just showing the liquid composition differences. * Run the 2D calculation with displaced liquid composition away from the solid (more reactive) * lower epsionl_phi and W keeping the ratio constant (This method does not effect the vapor liquid surface energy so that will probably have to be lowered as well, plottingh the paths on the phase diagram will help figure this out) * metrics such as dynamic contact angle.
Here are the actual dimensionless parameters for the current set of numbers,
Bill has identified the important number to be reduced with this set of parameters, }
Bill just gave me some new numbers:
The problem with these numbers is that the free energy does not seem to have a turning point for densities in the range rho1 = (0,100), rho2=(0,100).
Bill came up with a new set of numbers on Friday. The python equilibrium calculations eventually agreed with his partial densities. It was mostly a tolerance issue. To get agreement, it required extracting the exact numbers from the matematica notebook, not the numbers that Bill wrote down. Anyway, here are the numbers
These numbers give the following dimensionless numbers,
These numbers give a length scale of about .3 mm and the magnitude of the terms are all near 1 excepting some of the free energy terms that are near 10.
I've been using sympy do the symbolic manipulation for the dimensionless parameters and then to output in tex format. It's very powerful. Here are the numbers that Bill suggested. I've chosen a length, time and velocity scale that bring terms to order 1 in the phase field equation. Unfortunately, this leaves the other equations with large disparities in parameter sizes.
The following shows the dimensionless parameter for each term after division by the the transient prefactor for each equation. The X dependence of each term gives a clue to the type of term, whether convection, diffusion, 4th order or source. In the exercise above we chose \mathrm{\tau} \mathrm{\epsilon_1} {\mathrm{X}}^{(-4)}= 1.021\times 10^{04} \newline \mathrm{\bar{M}} \mathrm{\tau} \mathrm{\epsilon_2} {\mathrm{X}}^{(-4)}= 1.021\times 10^{04} \newline \mathrm{A_1} \mathrm{\bar{M}} \mathrm{\tau} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 3.146\times 10^{01} \newline \mathrm{A_2} \mathrm{\bar{M}} \mathrm{\tau} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 8.482\times 10^{01} \newline \mathrm{\bar{M}} \mathrm{\tau} \mathrm{e_1} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -8.061\times 10^{01} \newline \mathrm{\bar{M}} \mathrm{\tau} \mathrm{e_2} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -4.030\times 10^{01} \newline \mathrm{\bar{M}} R \mathrm{\tau} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-1)}= 8.994\times 10^{00} \newline
\text{equation}: \phi \newline V \mathrm{\tau} {\mathrm{X}}^{(-1)}= 1.000\times 10^{00} \newline \mathrm{\tau} \mathrm{\epsilon_{\phi}} \mathrm{M_{\phi}} {\mathrm{X}}^{(-2)}= 1.000\times 10^{00} \newline \mathrm{A_1} \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 1.891\times 10^{-01} \newline \mathrm{A_2} \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 5.100\times 10^{-01} \newline \mathrm{\tau} \mathrm{e_1} \mathrm{M_{\phi}} {\mathrm{\bar{\rho}}}^{2} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -4.846\times 10^{-01} \newline \mathrm{\tau} \mathrm{e_2} \mathrm{M_{\phi}} {\mathrm{\bar{\rho}}}^{2} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -2.423\times 10^{-01} \newline R \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} \mathrm{log}\left({\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{v_a}}^{(-1)} (\mathrm{m}- \mathrm{\bar{\rho}} \mathrm{\bar{v}})\right)= -7.079\times 10^{-01} \newline B \mathrm{\tau} \mathrm{M_{\phi}} \mathrm{m} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{T}}^{(-1)} {\mathrm{v_s}}^{(-2)}= 9.995\times 10^{-01} \newline B \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 1.001\times 10^{00} \newline 2 B \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{T}}^{(-1)} {\mathrm{v_s}}^{(-1)}= 2.000\times 10^{00} \newline W \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{T}}^{(-1)}= 1.507\times 10^{-02} \newline \text{equation}: \overrightarrow{v} \newline V \mathrm{\tau} {\mathrm{X}}^{(-1)}= 1.000\times 10^{00} \newline \mathrm{\tau} \mathrm{\mu_l} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)}= 2.553\times 10^{03} \newline \mathrm{\tau} \mathrm{\mu_s} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)}= 1.276\times 10^{18} \newline 2 B \mathrm{\tau} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{v_s}}^{(-1)}= 1.659\times 10^{05} \newline 2 B \mathrm{\tau} \mathrm{m} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-2)} {\mathrm{X}}^{(-1)} {\mathrm{v_s}}^{(-2)}= 1.659\times 10^{05} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{e_1} {V}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{m}}^{(-2)}= -4.021\times 10^{04} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{e_2} {V}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{m}}^{(-2)}= -2.010\times 10^{04} \newline R \mathrm{\tau} \mathrm{T} {V}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{m}}^{(-1)}= 4.486\times 10^{03} \newline R \mathrm{\tau} \mathrm{T} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{\bar{v}}}^{(-1)}= 5.198\times 10^{03} \newline \mathrm{\tau} \mathrm{\epsilon_{\phi}} \mathrm{T} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-3)}= 8.297\times 10^{04} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{\epsilon_1} \mathrm{T} {V}^{(-1)} {\mathrm{X}}^{(-3)}= 5.092\times 10^{06} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{\epsilon_1} \mathrm{T} {V}^{(-1)} {\mathrm{X}}^{(-3)}= 5.092\times 10^{06} \newline }}}
New parameters are as follows:
The equilibrium densities are = 2945.61370003, \;\;\; \rho_1^{\text{vapor}} = 13.2482193942, \;\;\; \newline \rho_2^{\text{solid}} = 1809.0511119, \;\;\; \rho_2^{\text{liquid}} = 4323.60642661, \;\;\; \rho_2^{\text{vapor}} = 148.486645391, }}} The dimensionless numbers are as follows. \mathrm{\tau} \mathrm{\epsilon_1} {\mathrm{X}}^{(-4)}= 1.000\times 10^{-00} \newline \mathrm{\bar{M}} \mathrm{\tau} \mathrm{\epsilon_2} {\mathrm{X}}^{(-4)}= 1.000\times 10^{-00} \newline \mathrm{A_1} \mathrm{\bar{M}} \mathrm{\tau} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 5.185\times 10^{-02} \newline \mathrm{A_2} \mathrm{\bar{M}} \mathrm{\tau} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 8.482\times 10^{-02} \newline \mathrm{\bar{M}} \mathrm{\tau} \mathrm{e_1} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -6.046\times 10^{-02} \newline \mathrm{\bar{M}} \mathrm{\tau} \mathrm{e_2} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -4.030\times 10^{-02} \newline \mathrm{\bar{M}} R \mathrm{\tau} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)} {\mathrm{m}}^{(-1)}= 8.994\times 10^{-03} \newline }}} \mathrm{M_{\phi}} {\mathrm{X}}^{(-2)}= 1.000\times 10^{02} \newline \mathrm{A_1} \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 3.182\times 10^{01} \newline \mathrm{A_2} \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 5.206\times 10^{01} \newline \mathrm{\tau} \mathrm{e_1} \mathrm{M_{\phi}} {\mathrm{\bar{\rho}}}^{2} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -3.710\times 10^{01} \newline \mathrm{\tau} \mathrm{e_2} \mathrm{M_{\phi}} {\mathrm{\bar{\rho}}}^{2} {\mathrm{m}}^{(-2)} {\mathrm{T}}^{(-1)}= -2.474\times 10^{01} \newline R \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} \mathrm{log}\left({\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{v_a}}^{(-1)} (\mathrm{m}- \mathrm{\bar{\rho}} \mathrm{\bar{v}})\right)= -7.226\times 10^{01} \newline B \mathrm{\tau} \mathrm{M_{\phi}} \mathrm{m} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{T}}^{(-1)} {\mathrm{v_s}}^{(-2)}= 1.020\times 10^{02} \newline B \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{m}}^{(-1)} {\mathrm{T}}^{(-1)}= 1.021\times 10^{02} \newline 2 B \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{T}}^{(-1)} {\mathrm{v_s}}^{(-1)}= 2.042\times 10^{02} \newline W \mathrm{\tau} \mathrm{M_{\phi}} {\mathrm{T}}^{(-1)}= 1.538\times 10^{01} }}} }^{(-1)} {\mathrm{X}}^{(-2)}= 2.553\times 10^{00} \newline \mathrm{\tau} \mathrm{\mu_s} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-2)}= 1.276\times 10^{15} \newline 2 B \mathrm{\tau} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{v_s}}^{(-1)}= 1.694\times 10^{04} \newline 2 B \mathrm{\tau} \mathrm{m} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-2)} {\mathrm{X}}^{(-1)} {\mathrm{v_s}}^{(-2)}= 1.693\times 10^{04} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{e_1} {V}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{m}}^{(-2)}= -3.079\times 10^{03} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{e_2} {V}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{m}}^{(-2)}= -2.052\times 10^{03} \newline R \mathrm{\tau} \mathrm{T} {V}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{m}}^{(-1)}= 4.580\times 10^{02} \newline R \mathrm{\tau} \mathrm{T} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-1)} {\mathrm{\bar{v}}}^{(-1)}= 5.306\times 10^{02} \newline \mathrm{\tau} \mathrm{\epsilon_{\phi}} \mathrm{T} {V}^{(-1)} {\mathrm{\bar{\rho}}}^{(-1)} {\mathrm{X}}^{(-3)}= 8.297\times 10^{03} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{\epsilon_1} \mathrm{T} {V}^{(-1)} {\mathrm{X}}^{(-3)}= 5.092\times 10^{04} \newline \mathrm{\bar{\rho}} \mathrm{\tau} \mathrm{\epsilon_1} \mathrm{T} {V}^{(-1)} {\mathrm{X}}^{(-3)}= 5.092\times 10^{04} }}} As can be seen from above the time step size was chosen at 1e-9. We can go a lot higher then this as we are using an implicit scheme. That was simply chosen above so that the diffusion terms are of the same order of magnitude as the convection terms (the velocity scale is calculated from the time scale so automatically adjusts). When the viscosity in the solid is set to be equal to the liquid, after awhile, I get highly unstable behavior in the solid. The magnitude of the viscous term is way lower than the source terms, so this is hardly suprising. When the visocity is raised by five orders of magnitude, everything is cool. For the time being I am going to stick with a viscosity of 2e+2. I have a hunch that I'll be able to lower this in the liquid and vapor and still get converged solutions. It is generally the case in the SIMPLE algorithm that the relaxation factor has to be reduced as the Reynold's number is increased. For the moment, I don't want to bother with the viscosity issue as the surface tension issue is more critical.
With dx = 1e-8 and L = 2e-6 the solid-liquid system gave the following equilibrium results,
The liquid-vapor system gives
The maximum time step for the liquid-vapor system is much lower than the maximum time step for the solid-liquid system. The solid-vapor system gives the following results
The solid vapor calculation was fairly inaccurate. I don't think I let this run for long enough and maybe the resolution needs to be higher. However, the surface tensions will not change appreciably. The phase space plot for these three equilibrium calculations is below. Image(phaseSpace-billsNewParams2.png, 60px) The following three plots show the free energy excess along the three dotted black lines parameterized using the local denisty but also as a function of phase field, which creates a contour plot. Image(liquid-vapor-contour-billsNewParams2.png, 60px) Image(solid-liquid-contour-billsNewParams2.png, 60px) Image(solid-vapor-contour-billsNewParams2.png, 60px) The following plots are perhaps more useful than the last three. They show the free energy excess plotted for phase field values of 0.0, 0.5 and 1.0 on top of the equilibrium paths. Image(freeEnergyExcess0.0.png, 60px) Image(freeEnergyExcess0.5.png, 60px) Image(freeEnergyExcess1.0.png, 60px)
We lowered B a lot to get the following numbers
We got the following phase diagram. Image(freeEnergyExcess.png, 60px) The surface tensions are
Thus, we have a non-zero contact angle. The contact angle is ~34 degrees assuming a flat solid surface.
The previous numbers only give an acute non-zero contact angle; the contact angle needs to be bigger. Also, the B value is far to low. In the following set of numbers the B value has been reduced less, epsilon_1 has been reduced and the solid liquid equilibrium positions have been separated on the phase diagram. This increases the surface tension of solid-liquid boundary and reduces the surface tension of the solid-vapor boundary. The numbers are as follows:
The results for each of the three 1D simulations are
The following should be observed from the results above:
* The solid-vapor system has not yet converged. * The difference between the surface tensions is unlikely to be improved even when convergence is achieved in the solid-vapor case * Most of the variance in the solid-vapor system is in the rho_1 direction and epsilon_1 is smaller, the grid density is not sufficient to resolve this. * Still get a converged solution just lacks accuracy, this may be sufficient for now * Convergence is relatively slow, maybe slightly too slow for 2D calculations. Remains to be seen. Need to examine time step and number of sweeps at each step. The above calculation takes 10 sweeps per step. This could well be reduced. * The above calculations suggest a contact angle of 70 +/- 5 degrees. * One problem for the 2D simulations could be the size of the box. It could be too small for adequate resolution of the drop.
The phase space plot is as follows. It is clear why there is now a finite contact angle. Image(freeEnergyExcess2.png, 60px)
* Figure out why 2D simulation is not running anymore * launch new 2D job * get 1D simulation running faster * larger time steps * less sweeps * profile * 2D simulation * profile * implement and test trilinos to reproduce Max's speed ups.
We have changed the free energy from + W \phi^2 \left( 1 - \phi \right)^2 }}} to + \tilde{W} \rho \phi^2 \left( 1 - \phi \right)^2 }}} The phase field and momentum equations are changed slightly. The new numbers I am running with are as follows
The surface tensions for each of the three 1D simulations are
This new formualtion has decreased the surface tension on the solid vapor interface.
New free energy and some new numbers from Bill.
$ \text{parameters} \newline v_a= 1.000\times 10^{00} \newline \mu_s= 2.000\times 10^{04} \newline M_{\phi}= 1.000\times 10^{04} \newline m= 1.180\times 10^{-01} \newline v_s= 1.507\times 10^{-05} \newline \bar{v}= 1.300\times 10^{-05} \newline \epsilon_{\phi}= 1.000\times 10^{-09} \newline T= 6.500\times 10^{02} \newline A_1= 2.826\times 10^{04} \newline A_2= 5.644\times 10^{04} \newline \bar{M}= 1.000\times 10^{-07} \newline B= 2.542\times 10^{04} \newline R= 8.314\times 10^{00} \newline W= 1.274\times 10^{05} \newline e_1= -4.560\times 10^{-01} \newline e_2= -4.560\times 10^{-01} \newline \mu_l= 2.000\times 10^{04} \newline \epsilon_2= 1.000\times 10^{-16} \newline \epsilon_1= 1.000\times 10^{-17} \newline $
The equilirium results were for the solid-liquid (dataWed-Sep-19-11:43:08-2007, out.2145.0)
$ \rho_1^s = 7554.8 (7489.1), \;\;\; \rho_1^l = 776.41 (767.67), \;\;\; \rho_2^s = 360.74 (349.29), \;\;\; \rho_2^l = 6618.0 (6586.67), \;\;\; \gamma_{sl}^1 = 26.014, \;\;\; \gamma_{sl}^2 = 25.931, \;\;\; \max \left( |v| \right) = 1.5649 \times 10^8 $
for the solid-vapor (dataWed-Sep-19-13:37:27-2007, out.2147.0)
$ \rho_1^s = 7299.8 (7489.1), \;\;\; \rho_1^v = 9.0488 (8.6488 ), \;\;\; \rho_2^s = 540.53 (349.29), \;\;\; \rho_2^v = 76.528 (74.207), \;\;\; \newline \gamma_{sv}^1 = 24.572, \;\;\; \gamma_{sv}^2 = 20.135 $
for the liquid-vapor (dataWed-Sep-19-15:58:59-2007, out.2162.0)
$ \rho_1^l = 656.31 (767.67)), \;\;\; \rho_1^v = 7.4282 (8.6488 ), \;\;\; \rho_2^l = 6698.1 (6586.67), \;\;\; \rho_2^v = 76.626 (74.207), \;\;\; \newline \gamma_{lv}^1 = 18.148, \;\;\; \gamma_{lv}^2 = 18.582 $
The solid-vapor and liquid-vapor velocites are far too big for good convergence. The approximated contact angle is (94.6, 108.17).
This page is getting too long and will be continued [wiki:ReactiveWetting2].