Concentration distribution in a phase separated mixture. Can't get the correct ODEs and boundary conditions

I wish to compute the equilibrium concentration distribution of a binary mixture that has phase separated. I start with writing the free energy as a functional depending of the concentration. I use the Beltrami identity to minimize the functional with the constraint that the total concentration is conserved. I end up with a differential equation that must be solved numerically.

Assuming I did the first part correctly, the problem is that I don't manage to find a set of differential equations with appropriate boundary conditions.

I use the following free energy density $f$:

$f(x)=\phi(x) \log(\phi(x)) + (1-\phi(x)) \log(1-\phi(x)) + \chi \phi(x)(1-\phi(x)) + \gamma (\phi'(x))^2$

with $x$ the space coordinate and $x \in [0,1]$, $\phi$ the concentration, $\chi$ the interaction parameter, and $\gamma$ a surface tension parameter. The first two terms account for the entropy of mixing, the third term the energy of interaction, and the last term accounts for the surface tension : it costs energy to have a variation in the concentration. This term has to be squared to conserve the symmetry left-right.

The total free energy is the functional $I$:

$I[\phi] = \int_{0}^{1} dx \space \space f(\phi(x),\phi'(x) )$

This functional has to be minimized with the constraint that the total concentration $\Phi_0$ is conserved:

$\int_0^1 \phi(x) dx = \Phi_0$

Therefore, using the Lagrange multiplier method, I have to minimize the functional $J$:

$J[\phi]= \int_{0}^{1} dx \space \space L(\phi(x),\phi'(x)) \quad$ with $\quad L= f - \lambda (\phi-\Phi_0)\quad$ and $\lambda$ a real constant.

$L$ has no explicit dependence in $x$ so I can use the Beltrami identity to find the extremum of $J$:

$L-\phi' \frac{\partial L}{\partial \phi'}=c \quad$ with $c$ a real constant.

This leads to the following differential equation:

$\phi'=\left( \frac{1}{\gamma} \left[ \phi(x) \log(\phi(x)) + (1-\phi(x)) \log(1-\phi(x)) + \chi \phi(x)(1-\phi(x)) + \lambda(\phi-\Phi_0)-c \right] \right)^{\frac{1}{2}}$

If I have done correctly so far, it's a matter of solving numerically a differential equation with the appropriate boundary conditions. But here comes my troubles.

So as I understand it, $\lambda, c$ appear here as parameters. They should be set in order to satisfy the boundary conditions and the mass conservation constraint. To enforce mass conservation, I can define a new function $h$:

$h(x)=\int_0^x \phi(s)ds \quad, \quad h'=\phi$

and the boundary condition $h(1)=\Phi_0$ enforce the mass conservation.

I intentionally omit to require $h(0)=0$ because it is already contained in the definition of $h$. So if I were to add this boundary condition, it would be redundant. Am I right?

To wrap up I now have a system of two first order differential equations and two parameters ($\lambda,c$)

$\phi'=...\phi,\lambda,c,... \\ h'=\phi$

To be consistent, I should find 4 boundary conditions. I only have $h(1)=\Phi_0$. I think I should impose no-flux boundary conditions : $\phi'(0)=\phi'(1)=0$, but I don't see how to impose it. Should I create another function $z$ like this?

$z=\phi' \\ z'=\phi''=......$

So I could impose $z(0)=z(1)=0$ ? Even-though the number of boundary conditions won't be enough.

I also see that I am missing one thing: the phase separated mixture can be concentrated in the left side and dilute in the right side. Or the opposite. This is one degree of freedom that I don't see how to handle.

• What, exactly, do you want to do? You write about a system undergoing phase separation, yet nowhere in your equations do you have time or dynamics. Why is it important that the global order parameter have a certain value? If you wish to study an interface (and say, compute the surface tension), you can just set up appropriate boundary conditions: At the far left you have phase 1, at the far right, phase 2. See e.g. borisv.lk.net/matsc597c-1997/phases/Lecture5/node5.html (if you don't expand the log terms, you'll have to do this numerically). Jan 11 '15 at 20:50
• Hi alarge thanks for helping. I want to find the equilibrium concentration distribution. Kinetic is not considered. I edited to clarify this point. "What, exactly, do you want to do?" I want to get the density profile. I am not studying something in particular like the interface. I just made this as an exercise and I try to solve it. You ask : "Why is it important that the global order parameter have a certain value?". I don't get it, where do say that? Jan 12 '15 at 11:15
• Well surely you have a profile in mind: For example, the free energy minimum would be everywhere a constant phase 0 or 1. By global order parameter I meant the total concentration. You fix this with a Lagrange multiplier (chemical potential of sorts). I don't exactly see a need for this in this particular problem, for if you want to see phase separation, you might as well fix the boundary conditions so that this is what you get (you don't have to, of course, but it's easier and most probably gets you what you really wanted, anyway, thus my previous question). Jan 12 '15 at 11:24
• I tried first not to fix the global concentration. The results I got where nonsense, like diverging concentrations... I concluded it was because I didn't fix the global concentration. Then I am not sure I understand your comment. For low interaction parameter, entropy dominates and the minimum free energy is attained with flat density profile. For high interaction parameter, energy dominates and the density profile that minimizes the free energy is not flat any more. I don't see why I should fix the boundaries, in my view they should rather be determined when solving the equations Jan 12 '15 at 12:16

This is not really an answer to your question of what is wrong with your equation, but rather probably the simplest way to solve the equation numerically.

Time dependent Ginzburg-Landau model B says $$\frac{\partial \varphi}{\partial t} = \nabla^2 \frac{\delta \mathcal{F}}{\delta \varphi}$$ Why use this equation? The stable state is $\frac{\delta \mathcal{F}}{\delta \varphi} = 0$, and the flow towards it, being diffusive, conserves the total order parameter.

This is trivial to implement numerically. I've attached code that does the simplest Euler integration using periodic boundary conditions (in Python). Note that some constants are slightly different from yours, so you ought to do the differentiation yourself just to make sure.

from scipy.ndimage import * # laplace operator

z = -5.
c = .01

TT  = 1000000
dt  = 1e-8
N   = 128

dx  = 1./(N-1)
dx2 = dx**2

x = ones(N) - 1e-1
x[:N/2] = 1e-1

def zfunx(x): # Functional derivative of F;
dFdm = 2*z*x + log(x/(1-x)) - c*laplace(x, mode='wrap')/dx2
return dFdm

for i in range(TT): # Euler integration
x = x + laplace(zfunx(x), mode='wrap')/dx2*dt


I do not recommend this method as it is slow to converge, but I hope this will get you started (look up semi-implicit Fourier space method for a more efficient algorithm). Here's what the solution looks like for this particular set of constants:

• Sorry I was working on something else for a while. Thank you for you answer. I try to understand it. I chose to minimize the free energy with the mass conservation constraint. Instead, you write $\nabla^2 \frac{\delta F}{\delta \phi}=0$ (I think you made a typo error) which I also understand. But the two things don't seem to be equivalent... or are they? Beside, in your code for zfunx(x) I don't quite understand : 2.z.x is the energy term, right? And c*laplace(x, mode='wrap')/dx2 for the surface tension? Why laplace? Jan 22 '15 at 20:05
• @David Yes, 2*z*x is the energy term. laplace is the second finite difference, with the (1 -2 1) kernel. Your concerns about uniqueness are valid (and I did have a typo), but I don't think they play a role: Suppose that $\nabla^2\mu=\nabla^2\frac{\delta\mathcal{F}}{\delta\varphi}=0$. This means that (in 1D) $\mu = ax+b$. I think we get from the boundary conditions that $a=0$, eliminating the problem. Jan 22 '15 at 20:49
• How do you impose the periodic boundary conditions? When I solve $\frac{d^2}{dx^2}\frac{\delta F}{\delta \phi}=0$, I get a 4th order differential equation. Then I try the periodic boundary conditions by requiring in Matlab (bvp4c) : $\phi(0)=\phi(1)$, $\phi'(0)=\phi'(1)$, $\phi''(0)=\phi''(1)$, $\phi'''(0)=\phi'''(1)$ Jan 23 '15 at 15:43
• But I get an error "Unable to solve the collocation equations -- a singular Jacobian encountered." Jan 23 '15 at 15:50
• @David It's all in the code: laplace wraps around (with mode = 'wrap'), so laplace(x)[0] = x[-1] - 2*x[0] + x[1] and laplace(x)[-1] = x[0] - 2*x[-1] + x[-2] (if you're not familiar with Python, -1 refers to the last element). As for your Matlab code, it's been years since I've last used Matlab (or bvp4c), so I can't really give you a quick answer there. Note that with the PBC, the problem does not have a unique solution: you can always shift the distribution left or right. This might be problematic for solvers, so you should specify one end as bulk phase 1 and the other as phase 2. Jan 23 '15 at 15:57