1
$\begingroup$

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 1 s

After 2 s:

After 2 s

Expecting:

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Daniel Shapero May 8 at 0:10
  • $\begingroup$ How are you performing the time integration? $\endgroup$ – whpowell96 May 8 at 2:03
  • 3
    $\begingroup$ 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. $\endgroup$ – MPIchael May 8 at 11:00
  • 1
    $\begingroup$ Seems like you are missing a factor of $1/\Delta x^2$ in your difference formula. $\endgroup$ – David Ketcheson May 9 at 10:34
  • 1
    $\begingroup$ 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 $\endgroup$ – whpowell96 May 12 at 20:39
1
$\begingroup$

I can recommend standard test for Cahn-Hilliard code for Python on

https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/cahn-hilliard/python/documentation.html#

and for Mathematica on

https://mathematica.stackexchange.com/questions/202446/solving-cahn-hilliard-equation-linearsolve-linear-equation-encountered-that-ha/202503#202503

Picture you are showing could be generated in this test. Note, in this test we solve pair of second order equation and not one of fourth order. Check this animation I made for the test

Figure 1

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.