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 η:
dφ/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:
∂/∂η ∫_Ω [ .5 dt J(η)·∇η - η (φ - φ0) ] = 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 [0]
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 [2]
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 [0]
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.dφ = 'φ - φ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[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.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
256 z6U0ps8zM/LPjspVypDKUsuLzk3MgM3OzYnN7s/61KfP2zH4MADNhst3z7DMoBcvyQ==''')
257 with self.subTest('chemical-potential'):
258 self.assertAlmostEqual64(state['η'], '''
259 eNoBYgCd/+s1Bcp4ztI3gjYFyZk4YzVjyfA2AzdAMj032zfLNTE4fMm7yLnGisbqxZPJ2MsfyD81csiv
260 x+E5xDhjOJA3msZ1xZTFa8ddx/fG88eCx73H1MieM/c0WDihMUrLvMYZNpvIrWQ0sw==''')