cahnhilliard.py¶

In this script we solve the Cahn-Hilliard equation, which models the unmixing of two phases (φ=+1 and φ=-1) under the influence of surface tension. It is a mixed equation of two scalar equations for phase φ and chemical potential η:

/dt = -div(J(η))

η = ψ'(φ) σ / ε - σ ε Δ(φ)

along with constitutive relations for the flux vector J = -M ∇η and the double well potential ψ = .25 (φ² - 1)², and subject to boundary conditions ∂ₙφ = -σd / σ ε and ∂ₙη = 0. Parameters are the interface thickness ε, fluid surface tension σ, differential wall surface tension σd, and mobility M.

Cahn-Hilliard is a diffuse interface model, which means that phases do not separate sharply, but instead form a transition zone between the phases. The transition zone has a thickness proportional to ε, as is readily confirmed in one dimension, where a steady state solution on an infinite domain is formed by η(x) = 0, φ(x) = tanh(x / √2 ε).

The main driver of the unmixing process is the double well potential ψ that is proportional to the mixing energy, taking its minima at the pure phases φ=+1 and φ=-1. The interface between the phases adds an energy contribution proportional to its length. At the wall we have a phase-dependent fluid-solid energy. Over time, the system minimizes the total energy:

E(φ) = ∫_Ω ψ(φ) σ / ε + ∫_Ω .5 σ ε ‖∇φ‖² + ∫_Γ (σm + φ σd)

\                \                  \

mixing energy    interface energy   wall energy

Proof: the time derivative of E followed by substitution of the strong form and boundary conditions yields dE/dt = ∫_Ω η dφ/dt = -∫_Ω M ‖∇η‖² ≤ 0. □

Switching to discrete time we set dφ/dt = (φ - φ0) / dt and add a stabilizing perturbation term δψ(φ, φ0) to the doube well potential for reasons outlined below. This yields the following time discrete system:

φ = φ0 - dt div(J(η))

η = (ψ'(φ) + δψ'(φ, ψ0)) σ / ε - σ ε Δ(φ)

with the equivalent weak formulation:

∂/∂η ∫_Ω [ η (φ - φ0) + .5 dt J(η)·∇η ] = 0

∂/∂φ ∫_Ω [ E(φ) + δψ(φ, φ0) σ / ε - η φ ] = 0

For stability we wish for the perturbation δψ to be such that the time discrete system preserves the energy dissipation property E(φ) ≤ E(φ0) for any timestep dt. To derive suitable perturbation terms to this effect, we define without loss of generality δψ’(φ, φ0) = .5 (φ - φ0) f(φ, φ0) and derive the following condition for unconditional stability:

E(φ) - E(φ0) = ∫_Ω .5 (1 - φ² - .5 (φ + φ0)² - f(φ, φ0)) (φ - φ0)² σ / ε

- ∫_Ω (.5 σ ε ‖∇φ - ∇φ0‖² + dt M ‖∇η‖²) ≤ 0

The inequality holds true if the perturbation f is bounded from below such that f(φ, φ0) ≥ 1 - φ² - .5 (φ + φ0)². To keep the energy minima at the pure phases we additionally impose that f(±1, φ0) = 0, and select 1 - φ² as a suitable upper bound which we will call “nonlinear”.

We next observe that η is linear in φ if f(φ, φ0) = g(φ0) - φ² - (φ + φ0)² for any function g, which dominates if g(φ0) ≥ 1 + .5 (φ + φ0)². While this cannot be made to hold true for all φ, we make it hold for -√2 ≤ φ, φ0 ≤ √2 by defining g(φ0) = 2 + 2 ‖φ0‖ + φ0², which we will call “linear”. This scheme further satisfies a weak minima preservation f(±1, ±‖φ0‖) = 0.

We have thus arrived at the three stabilization schemes implemented here:

• nonlinear: f(φ, φ0) = 1 - φ²

• linear: f(φ, φ0) = 2 + 2 ‖φ0‖ - 2 φ (φ + φ0)

• none: f(φ, φ0) = 0 (not unconditionally stable)

The stab enum in this script defines the schemes in terms of δψ to allow Nutils to construct the residuals through automatic differentiation.

NOTE: This script uses dimensional quantities and requires the nutils.SI module, which is installed separate from the the core nutils.

78from nutils import mesh, function, solver, numeric, export, cli, testing
79from nutils.expression_v2 import Namespace
80from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse
81import numpy
82import itertools
83import enum
84import treelog as log
85
86
87class stab(enum.Enum):
88    nonlinear = '.25 dφ^2 (1 - φ^2 + φ dφ (2 / 3) - dφ^2 / 6)'
89    linear = '.25 dφ^2 (2 + 2 abs(φ0) - (φ + φ0)^2)'
90    none = '0'  # not unconditionally stable
91
92
93def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
94         wtensn: Tension, wtensp: Tension, nelems: int, etype: str, degree: int,
95         timestep: Time, tol: Energy/Length, endtime: Time, seed: int, circle: bool,
96         stab: stab, showflux: bool):
97    '''
98    Cahn-Hilliard equation on a unit square/circle.
99
100    .. arguments::
101
102       size [10cm]
103         Domain size.
104       epsilon [2mm]
105         Interface thickness; defaults to an automatic value based on the
106         configured mesh density if left unspecified.
107       mobility [1mL*s/kg]
108         Mobility.
109       stens [50mN/m]
110         Surface tension.
111       wtensn [30mN/m]
112         Wall surface tension for phase -1.
113       wtensp [20mN/m]
114         Wall surface tension for phase +1.
115       nelems 
116         Number of elements along domain edge. When set to zero a value is set
117         automatically based on the configured domain size and epsilon.
118       etype [square]
119         Type of elements (square/triangle/mixed).
120       degree 
121         Polynomial degree.
122       timestep [.5s]
123         Time step.
124       tol [1nJ/m]
125         Newton tolerance.
126       endtime [1min]
127         End of the simulation.
128       seed 
129         Random seed for the initial condition.
130       circle [no]
131         Select circular domain as opposed to a unit square.
132       stab [linear]
133         Stabilization method (linear/nonlinear/none).
134       showflux [no]
135         Overlay flux vectors on phase plot
136    '''
137
138    nmin = round(size / epsilon)
139    if nelems <= 0:
140        nelems = nmin
141        log.info('setting nelems to {}'.format(nelems))
142    elif nelems < nmin:
143        log.warning('mesh is too coarse, consider increasing nelems to {:.0f}'.format(nmin))
144
145    log.info('contact angle: {:.0f}°'.format(numpy.arccos((wtensn - wtensp) / stens) * 180 / numpy.pi))
146
147    domain, geom = mesh.unitsquare(nelems, etype)
148    if circle:
149        angle = (geom-.5) * (numpy.pi/2)
150        geom = .5 + function.sin(angle) * function.cos(angle)[[1, 0]] / numpy.sqrt(2)
151
152    bezier = domain.sample('bezier', 5)  # sample for surface plots
153    grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5)  # sample for quivers
154
155    φbasis = ηbasis = domain.basis('std', degree=degree)
156    ηbasis *= stens / epsilon  # basis scaling to give η the required unit
157
158    ns = Namespace()
159    ns.x = size * geom
160    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
161    ns.ε = epsilon
162    ns.σ = stens
163    ns.φ = function.dotarg('φ', φbasis)
164    ns.σmean = (wtensp + wtensn) / 2
165    ns.σdiff = (wtensp - wtensn) / 2
166    ns.σwall = 'σmean + φ σdiff'
167    ns.φ0 = function.dotarg('φ0', φbasis)
168    ns. = 'φ - φ0'
169    ns.η = function.dotarg('η', ηbasis)
170    ns.ψ = '.25 (φ^2 - 1)^2'
171    ns.δψ = stab.value
172    ns.M = mobility
173    ns.J_i = '-M ∇_i(η)'
174    ns.dt = timestep
175
176    nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=7)
177    nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=7)
178    nrg_wall = domain.boundary.integral('σwall dS' @ ns, degree=7)
179    nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=7)
180
181    numpy.random.seed(seed)
182    state = dict(φ=numpy.random.normal(0, .5, φbasis.shape))  # initial condition
183
184    with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps:
185        for istep in steps:
186
187            E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **state))
188            log.user('energy: {0:,.0μJ/m} ({1:.0f}% mixture, {1:.0f}% interface, {1:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E)))
189
190            state['φ0'] = state['φ']
191            state = solver.optimize(['φ', 'η'], nrg / tol, arguments=state, tol=1)
192
193            with export.mplfigure('phase.png') as fig:
194                ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]')
195                x, φ = bezier.eval(['x_i', 'φ'] @ ns, **state)
196                im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='bwr')
197                im.set_clim(-1, 1)
198                fig.colorbar(im)
199                if showflux:
200                    x, J = grid.eval(['x_i', 'J_i'] @ ns, **state)
201                    log.info('largest flux: {:.1mm/h}'.format(numpy.max(numpy.hypot(J[:, 0], J[:, 1]))))
202                    ax.quiver(*(x/'mm').T, *(J/'m/s').T, color='r')
203                    ax.quiver(*(x/'mm').T, *-(J/'m/s').T, color='b')
204                ax.autoscale(enable=True, axis='both', tight=True)
205
206    return state

If the script is executed (as opposed to imported), nutils.cli.run() calls the main function with arguments provided from the command line.

212if __name__ == '__main__':
213    cli.run(main)

Once a simulation is developed and tested, it is good practice to save a few strategic return values for regression testing. The nutils.testing module, which builds on the standard unittest framework, facilitates this by providing nutils.testing.TestCase.assertAlmostEqual64() for the embedding of desired results as compressed base64 data.

222class test(testing.TestCase):
223
224    def test_initial(self):
225        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
226                     stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
227                     etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
228                     endtime=parse('1h'), seed=0, circle=False, stab=stab.linear, showflux=True)
229        with self.subTest('concentration'):
230            self.assertAlmostEqual64(state['φ0'], '''
231                eNoBYgCd/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx
232                L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJMF4tqQ==''')
233
234    def test_square(self):
235        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
236                     stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
237                     etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
238                     endtime=parse('2h'), seed=0, circle=False, stab=stab.linear, showflux=True)
239        with self.subTest('concentration'):
240            self.assertAlmostEqual64(state['φ'], '''
241                eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk
242                ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''')
243        with self.subTest('chemical-potential'):
244            self.assertAlmostEqual64(state['η'], '''
245                eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0
246                xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''')
247
248    def test_mixedcircle(self):
249        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
250                     stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
251                     etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
252                     endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
253        with self.subTest('concentration'):
254            self.assertAlmostEqual64(state['φ'], '''
255                eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK