# Why is my Cahn-Hilliard simulation separating out so finely?

I am trying to simulate the Cahn-Hilliard equation using Python, but the 2 fluids aren't separating into big blobs, as desired, under any conditions. I'm setting up (what I think is) an orthogonal mesh with uniform concentration values of 0.5 plus random noise values within ±0.05. $$\gamma = 0.5$$. $$D$$ = 50. Time steps of 1 'second'. Total time of 10 'seconds'.

The first image shows the state at $$t$$ = 1s. The second image shows it at $$t$$ = 2s, after which the system ceases to evolve. When zoomed in, Image 2 looks a bit like what I expected but its too fine; I'm hoping for bigger clumps as in Expecting:

After 1 s:

After 2 s:

Expecting:

To find the second derivative I am using (for $$\frac{d^2}{dx^2}$$):

mesh_out[y][x] = mesh_in[y][x+2]) - (2 * mesh_in[y][x]) + mesh_in[y][x-2]


After each time step if the concentration at a point is greater than 1 is set it back to 1. Likewise, if it's less than -1 I set it to -1.

I have tried messing about with values of the time step, $$D$$, $$\gamma$$, initial concentration but they seem to make no difference. I am not using any interpolation - should I be? (If so can you point me to any relevant information?) Otherwise, I'm not sure why this isn't producing the results I'm after. My suspicion is that I should be using interpolation but my not sure why this would make a difference. Do I need an energy minimising function of some sort? I've seen people mentioning the Gibbs free energy, for example. Any help/suggestions would be appreciated.

I haven't posted any code as I don't know which bits will be helpful.

• It strikes me as a little curious that you're using a finite difference formula that gets mesh points 2 cells away for the 2nd derivative rather than 1, but without seeing more of your code there's only so much we can help. The linear part of the Cahn-Hilliard equation is just the biharmonic heat equation -- I would try writing a solver for that first and verifying against an analytical solution obtained via Fourier analysis. – Daniel Shapero May 8 at 0:10
• How are you performing the time integration? – whpowell96 May 8 at 2:03
• That stencil looks suspicious to me. Have you double checked if the stencils you use for derivatives check out? You could fill your field with a sine(xy) and compare the results of your stencils against the (analytical) derivative. – MPIchael May 8 at 11:00
• Seems like you are missing a factor of $1/\Delta x^2$ in your difference formula. – David Ketcheson May 9 at 10:34
• What I am asking is how are you solving the IVP forward in time. This problem is stiff so some care may need to be taken in algorithm choice but I would suggest fixing the other issues people have mentioned as they seem more likely to be the culprit – whpowell96 May 12 at 20:39