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.

Improve this question
$\endgroup$
8
  • 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, 2020 at 0:10
  • $\begingroup$ How are you performing the time integration? $\endgroup$
    – whpowell96
    May 8, 2020 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, 2020 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, 2020 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, 2020 at 20:39

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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