Convergence of a Robin boundary approximation for a Cahn--Hilliard system with dynamic boundary conditions
Patrik Knopf, Kei Fong Lam
TL;DR
This paper proves the existence and uniqueness of solutions for a generalized Cahn--Hilliard system with dynamic boundary conditions, including a Robin boundary approximation, and analyzes convergence and error estimates with supporting numerical experiments.
Contribution
It introduces a Robin boundary approximation for a Cahn--Hilliard system with dynamic boundary conditions and establishes convergence and error estimates for the approximation.
Findings
Existence and uniqueness of weak solutions for the Robin boundary model.
Weak convergence of solutions as the regularization parameter tends to zero.
Numerical experiments reveal complex dynamics not seen in the original model.
Abstract
We prove the existence of unique weak solutions to an extension of a Cahn--Hilliard model proposed recently by C.~Liu and H.~Wu (2019), in which the new dynamic boundary condition is further generalised with an affine linear relation between the surface and bulk order parameters. As a first approach to tackle more general and nonlinear relations, we investigate the existence of unique weak solutions to a regularisation by a Robin boundary condition. Included in our analysis is the case where there is no diffusion for the surface order parameter, which causes new difficulties for the analysis of the Robin system. Furthermore, for the case of affine linear relations, we show the weak convergence of solutions as the regularisation parameter tends to zero, and derive an error estimate between the two models. This is supported by numerical experiments which also demonstrate some non-trivial…
| - | - | - | ||||
| - | - | - | ||||
| - | - | - | ||||
| - | - | - | ||||
| - | - | - | ||||
| - | - | - | ||||
| - | - | - |
| - | - | - | ||||
| - | - | - |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Convergence of a Robin boundary approximation for a Cahn–Hilliard system with dynamic boundary conditions
Patrik Knopf 111Department of Mathematics, University of Regensburg, 93053, Germany ([email protected])
Kei Fong Lam 222Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong ([email protected])
Abstract
We prove the existence of unique weak solutions to an extension of a Cahn–Hilliard model proposed recently by C. Liu and H. Wu (2019), in which the new dynamic boundary condition is further generalised with an affine linear relation between the surface and bulk phase field variables. As a first approach to tackle more general and nonlinear relations, we investigate the existence of unique weak solutions to a regularisation by a Robin boundary condition. Included in our analysis is the case where there is no diffusion for the surface phase field, which causes new difficulties for the analysis of the Robin system. Furthermore, for the case of affine linear relations, we show the weak convergence of solutions as the regularisation parameter tends to zero, and derive an error estimate between the two models. This is supported by numerical experiments which also demonstrate some non-trivial dynamics for the extended Liu–Wu model that is not present in the original model.
This is a preprint version of the paper. Please cite as:
P. Knopf and K.F. Lam, Nonlinearity 33(8): 4191-4235, 2020
https://doi.org/10.1088/1361-6544/ab8351
**Key words. ** Cahn–Hilliard equation; Dynamic boundary conditions; Penalisation via Robin boundary conditions; Gradient flow
AMS subject classification. 35A01, 35A02, 35A35, 35B40
1 Introduction
Since its introduction by physicists [13, 14, 23], the Cahn–Hilliard equation with dynamic boundary conditions has been the subject of many analytical and numerical studies, see, for instance, [25, 27] for a recent overview. The modification from the standard zero Neumann boundary conditions to dynamic boundary conditions allows for the inclusion of effective short-range interactions between the boundary of the domain and the mixture contained within. One application in which these interactions cannot be neglected lies in the modelling of moving contact lines in binary fluid mixtures [16, 17, 28], where instead of the fixed contact angle of imposed by the typical Neumann conditions for the Cahn–Hilliard component, it has been proposed to employ a dynamic boundary condition that allows some transport on the domain boundary to better emulate the dynamics of the moving contact line.
The term dynamic boundary conditions in fact encapsulates a class of boundary conditions, exhibiting a common structure involving a surface time-dependent partial differential equation (PDE) containing the normal derivative of the bulk quantity. Many such dynamic boundary conditions can be derived as the variation of suitable surface free energies. To fix ideas, let be an arbitrary positive real number and let be a bounded domain with sufficiently smooth boundary . We denote the normal derivative as with the unit outer normal vector on . For fixed , the bulk Ginzburg–Landau free energy functional is defined as
[TABLE]
where is a double well potential with two equal minima and denotes the phase field variable. The Cahn–Hilliard equation is the gradient flow of the above free energy:
[TABLE]
typically furnished with an initial condition in . The condition arises naturally from taking the variation of , while the condition is chosen to guarantee conservation of mass. To account for interactions with the domain boundary , we prescribe a surface free energy functional
[TABLE]
where is a fixed constant, is a surface potential function, and denotes the surface gradient on . By performing variation of with respect to suitable inner products, various dynamic boundary conditions replacing can be obtained. For instance, the dynamic boundary condition of Allen–Cahn type
[TABLE]
proposed in [13, 14, 23] can be regarded formally as the -gradient flow of , where is the Laplace–Beltrami operator and the bulk effects are captured by the term . Another dynamic boundary condition, that was proposed in [20] and studied in [4, 7], is of Cahn–Hilliard type,
[TABLE]
for nonnegative constants , and arises from taking the subdifferential of the total energy with respect to a suitable scalar product. In the recent work [25], Liu and Wu employed the energetic variational approach to derive a new type of dynamic boundary condition for the Cahn–Hilliard equation. This approach combines the least action principle and Onsager’s principle of maximum energy dissipation to ensure that the resulting model fulfils the physical constraints of mass conservation, dissipation of energy and balance of forces. Our present interest stems from the observation that we can extend the model and derivation of Liu and Wu to the following Cahn–Hilliard system with dynamic boundary conditions:
[TABLE]
Here, , , and are fixed constants, and are the bulk and surface phase field variables, and and are the corresponding bulk and surface chemical potentials. For applications, we may consider as the difference in volume fractions between two species of materials in the bulk domain, and as the difference in volume fractions between two species of materials restricted to the domain boundary. Introducing the variable allows us to express (1.1) equivalently as
[TABLE]
which coincides with the model of [9, 25] when and . Our consideration of the affine linear transmission condition as oppose to the classical relation is motivated in part by observations of some non-trivial dynamics that is not present in the original model with and , see Section 7 below. For example, if and , then the surface phase is the opposite of the bulk phase . If now and have equal minima at , we immediately see that is a solution to (1.2) which cannot be achieved for the usual dynamic boundary condition.
Associated to (1.2) is the energy functional
[TABLE]
In the appendix, we outline the derivation of (1.2) by the energetic variational approach with the help of . In Section 3, we show that (1.2) (and hence (1.1)) can be interpreted as the gradient flow of with respect to a suitable inner product. This allows us to employ a natural implicit time discretisation to prove, akin to the first author’s previous work [19], the existence of a unique global weak solution to (1.1).
A further extension would be to replace the classical relation with the more general situation for some continuous function , and we attempt to tackle this through the well known method of penalising this Dirichlet-like condition by a Robin boundary approximation. More precisely, for , we consider the system
[TABLE]
Formally, as , we recover the desired relation . In the case with and , this strategy has been successful for the Allen–Cahn system with rigorous convergence results and error estimates derived in [8]. In this work we aim to deduce analogous assertions for the more difficult Cahn–Hilliard model (1.3).
Associated to (1.3) is the energy functional
[TABLE]
for , and likewise, (1.3) can be interpreted as the gradient flow of with respect to a suitable inner product. Thus, a natural implicit time discretisation is used to show the existence of a unique global weak solution to (1.3) for a large class of nonlinear relation , see (). While one may expect that, as a regularisation, solutions to (1.3) exhibit better regularity properties than solutions to the limit system (1.1), in fact the opposite is true. This means we encounter more difficulties, especially in the case , than in the study of the limit model (1.1). The chief reason is that the surface variable is not just given by the trace of the bulk variable , but rather coupled through the Robin boundary condition . Therefore, the regularity of cannot be deduced by the trace theorem with the help of as it was done in [19, 25]. We have to apply further techniques and restrictions to overcome these difficulties.
Aside from existence and uniqueness of global weak solutions to (1.1) and (1.3), we investigate the rigorous convergence as of solutions for the case and we derive an error estimate of the form (under the same initial data for )
[TABLE]
which implies that the affine linear transmission condition is attained from the approximation solutions to (1.3) at a linear rate in . This is supported by numerical experiments performed for selected values of in the affine linear case. Moreover, we have also observed the same rate of convergence in numerical experiments for two examples of nonlinear relations .
The structure of this paper is as follows: in Section 2 we outline several preliminary results essential for our investigations and present the main results of this work. The gradient flow structure for both Cahn–Hilliard systems is demonstrated in Section 3, while the existence and uniqueness of weak solutions to (1.1) and (1.3) is established via an implicit time discretisation in Section 4 and 5, respectively. In Section 6, we discuss the convergence of (1.3) to (1.1) for the case and in Section 7 we present finite element discretisations of (1.2) and (1.3), numerical evidence in support of the error estimate (1.4), as well as a comparison between the dynamics of solutions to (1.2) for different values of and .
2 Preliminaries and Main results
Notation.
Throughout this paper we use the following notation: For any and , the standard Lebesgue and Sobolev spaces defined on are denoted as and , along with the norms and . For the case , these spaces become Hilbert spaces and we use the notation . A similar notation is used for Lebesgue and Sobolev spaces on . For any Banach space , we denote its dual space by and the associated duality pairing by . If is a Hilbert space, we denote its inner product by . We define
[TABLE]
as the spatial mean of , where denotes the -dimensional Lebesgue measure of . The spatial mean for and can be defined analogously. For brevity, we also use the notation
[TABLE]
for all and .
In the sections 2-6 we set the parameters in both models (1.1) and (1.3) as their values have no impact on the analysis we will carry out.
Assumptions.
We take with to be a bounded domain with Lipschitz boundary . 2.
We assume that the potentials and are bounded from below by constants such that
[TABLE]
and exhibit a decomposition in the form
[TABLE]
satisfying
. 2.
and are convex nonnegative functions. 3.
For any , there exist positive constants and such that
[TABLE] 4.
and are Lipschitz continuous. 3.
The function belongs to and there exist and a positive constant such that for all ,
[TABLE]
Remark 2.1**.**
The assumption (\mathrm{A2}$$\mathrm{A2iv}) implies there exist positive constants and such that for all ,
[TABLE]
Remark 2.2**.**
For our notion of weak solutions to (1.1) and (1.3), cf. [19, Defn. 4.1], it suffices to have a Lipschitz boundary , where the definition of tangential gradient on a Lipschitz surface can be found in [2, Defn. 3.1]. For the derivation of the error estimate (1.4), we will assume a smooth boundary in order to employ classical elliptic regularity results, which also is consistent with the set-up of [25].
Preliminaries.
For or , let and . Then, the unique solvability of the problems
[TABLE]
induces the solution operators , and , , satisfying
[TABLE] 2.
In this paper, the dual spaces and are endowed with the norms
[TABLE]
which are equivalent to the standard norms on these spaces. The spaces and are endowed with the same norms, i.e.,
[TABLE]
for and . 3.
Suppose that along with solutions and . Then, it holds that in the sense of distributions, which can be inferred by testing the Neumann problem for with a test function , where and are arbitrary, integrating by parts in time and employing the uniqueness of solutions. With this characterisation one also arrives at
[TABLE]
and analogously, for with , the solutions and satisfy in the sense of distributions and the equality
[TABLE] 4.
For any we define
[TABLE]
For , there exist and such that
[TABLE]
Hence, we can define the inner product
[TABLE]
for functions . The induced norm is given by . If actually lies in , the functions and can be expressed by
[TABLE]
We remark that is also a norm on . However, is not complete with respect to this norm. 5.
For fixed we define the Hilbert space
[TABLE]
endowed with the inner product and induced norm
[TABLE]
For any , we also define the subsets
[TABLE]
For , there exist and such that
[TABLE]
Thus, we can define an inner product on by
[TABLE]
Its induced norm is given by . If lies in , the functions and can be expressed by
[TABLE]
It holds that is also a norm on but is not complete with respect to this norm. 6.
A closer examination of the proof of the trace theorem (cf. [10]) shows that there exists a constant depending only on and such that for all ,
[TABLE]
The main results of this paper are listed as follows.
Theorem 2.1** (Global existence and uniqueness for (1.2)).**
Let , , , and be arbitrary and suppose that () and () hold. Then, for any initial condition satisfying and , there exists a unique weak solution to (1.1) in the following sense:
- (i)
The solution has the following regularity
[TABLE] 2. (ii)
For all , and with , and for a.e. , the following weak formulation
[TABLE]
and the initial condition are satisfied. 3. (iii)
For all , the energy inequality is satisfied
[TABLE]
In addition, if the boundary is smooth and there exist positive constants and such that for all ,
[TABLE]
then for any it holds that . Moreover, if , then .
Theorem 2.2** (Global existence and uniqueness for (1.3)).**
Let , , and be arbitrary. Suppose that ()-() hold, and in the case , the following assumptions hold additionally
* and there exist a real number , and constants and such that*
[TABLE] 2.
* for all where and are constants.*
[TABLE]
Then, for any initial condition with and satisfying and , there exists a unique weak solution to (1.3) in the following sense:
- (i)
The solution has the following regularity
[TABLE] 2. (ii)
For all , , , , and for a.e. , the following weak formulation
[TABLE]
and the initial conditions and are satisfied. 3. (iii)
For all , the energy inequality is satisfied
[TABLE]
Theorem 2.3** (Convergence for and error estimates).**
Suppose that (), () and () hold and let be as defined in (2.7). For arbitrary , , , , and , let , and denote the weak solution to (1.3) on with initial data satisfying
[TABLE]
Then, along a nonrelabelled subsequence,
[TABLE]
where is the weak solution to (1.2) on with initial data in the sense of items (i) and (ii) of Theorem 2.1.
Moreover, if (2.6) holds, then for any there exists a positive constant independent of and the solution variables such that
[TABLE]
For the case , the above estimate (2.14) holds with the modification where we remove the and estimates for .
3 Gradient flow structure
In this section we show that the Cahn–Hilliard systems (1.1) and (1.3) can be expressed as gradient flows of the energy functionals and with respect to suitable inner products.
3.1 Limit system (1.2)
The first variation of with respect to in the direction with , where we recall the subspace defined in () with dual space , reads as
[TABLE]
We introduce the solution operator which is defined by
[TABLE]
for . Then, for any , the functions , , and are well-defined, and so we can define the inner product
[TABLE]
with norm .
We seek to identify the gradient of at associated to this inner product, which satisfies for all with . In particular, this yields for the identity
[TABLE]
Following ideas in [19, Sec. 4.2], for arbitrary with and arbitrary nonnegative not identically zero, we define
[TABLE]
so that
[TABLE]
Defining the constants
[TABLE]
that are independent of and plugging the choice (3.3) into (3.1) yields the following identifications after a short calculation via the fundamental lemma of calculus of variations:
[TABLE]
Thus, the gradient flow
[TABLE]
can now be expressed as
[TABLE]
Choosing as in (3.3) we obtain the identification
[TABLE]
where the constants and independent of are defined as
[TABLE]
Using the definition of and , the system comprising of (3.6) and (3.7) is exactly the Cahn–Hilliard system (1.2).
3.2 Robin system (1.3)
The first variation of with respect to in direction reads as
[TABLE]
Recalling the inner product defined in (2.5), we seek the gradient of at associated to this inner product, which satisfies for all . In particular, this yields for the identity
[TABLE]
For arbitrary test functions , , we observe that . Defining constants
[TABLE]
independent of and plugging the choice of into (3.2) yields the following identifications after a short calculation:
[TABLE]
Thus, the gradient flow
[TABLE]
can now be expressed as
[TABLE]
Choosing and for arbitrary and , we obtain the identifications
[TABLE]
Using the definition of and , the system comprising of (3.11) and (3.13) is exactly the Cahn–Hilliard system (1.3).
4 Proof of Theorem 2.1
We define
[TABLE]
and it is easy to see that the assumption () implies that satisfies [19, (P2)], namely, is bounded from below and admits a decomposition where is convex, nonnegative and satisfies for any ,
[TABLE]
for some positive constant depending only on , while is Lipschitz continuous
[TABLE]
where is the Lipschitz constant of .
Since the main ideas have been outlined in [19], we simply sketch the crucial steps of the proof.
Step 1: Implicit time discretisation.
Let be arbitrary and denote a time step size. For , we define functions recursively by the following construction. The zeroth order iterate is taken to be the initial data, i.e., . If is already constructed, then we choose to be a minimiser of the functional
[TABLE]
on the set . We point out that may attain the value , since the -regularity for may not guarantee that is finite if the potentials and grow sufficiently fast. However, any minimiser to will have finite energy due to its construction. By following the proof of [19, Lem. 5] we find that is bounded from below for all and via the direct method of the calculus of variations, we infer the existence of a global minimiser to on the space . Hence, the existence of is guaranteed. Employing the convexity of and , we proceed as in [18] to obtain the Euler–Lagrange equation
[TABLE]
for any with . Setting
[TABLE]
so that by (), and are, respectively, the unique solutions to the Poisson problems
[TABLE]
it holds that
[TABLE]
Furthermore, for any arbitrary nonnegative not identically zero, and any arbitrary with , we define with where and are the constants defined in (3.2), as well as constants
[TABLE]
that are independent of . Then the functions
[TABLE]
satisfy
[TABLE]
for all , and \tilde{\eta}\in\mathcal{V}^{\kappa}\color[rgb]{0,0,1}\cap L^{\infty}(\Omega) with . This leads to an implicit time discretisation of (1.2).
Hence, the collection describes a time-discrete approximate solution. For and , we define the piecewise constant extension
[TABLE]
for , and the piecewise linear extension
[TABLE]
for any and .
Step 2: Uniform estimates.
Following [19, Sec. 4.4 and 4.5] we immediately infer the following uniform estimates: There exists a positive constant independent of , and , such that for any ,
[TABLE]
where is defined in (2.7). Then, from the uniform boundedness of and , we also deduce the boundedness of in and in .
Step 3: Passing to the limit.
Proceeding as in [19, Sec. 4.6], we obtain the existence of limit functions in the limit possessing the regularities stated in item (i) of Theorem 2.1. Items (ii) and (iii) of Theorem 2.1 can also be shown analogously to [19, Sec. 4.7], and the uniqueness of the weak solution follows as in [19, Sec. 5].
Step 4: Higher regularity.
We now establish the -regularity assertion by following the arguments outlined in [3, Sec. 4]. The new hypothesis (2.6) implies that , and so the third equality in item (ii) of Theorem 2.1 can be made to hold for all . Now, we choose and obtain the elliptic equation
[TABLE]
holding in the sense of distributions in , with a right-hand side belonging to . Elliptic regularity theory (see for instance [8, Thm. A.2] with , , , , , ) yields that for ,
[TABLE]
Since and thus , by a variant of the trace theorem involving elliptic operators (see for instance [8, Thm. A.1] with , , ), we have . Then, using (4.2) and the third equality in item (ii) of Theorem 2.1 we infer the surface elliptic equation
[TABLE]
holding in the sense of distributions on , with the right-hand side belonging to . Hence, by elliptic regularity theory, we have , and when revisiting the elliptic equation (4.2) now with a more regular trace , we have for ,
[TABLE]
For the case we take note of the Sobolev embedding , so that it follows from (4.3) that . Then, by elliptic regularity theory (see [8, Thm. A.2] with , , , , , ), we obtain . By the trace theorem, and hence, due to (2.6), it holds that . From (4.3) we now have and by [8, Thm. A.2] (with ) it holds that . This finishes the proof of Theorem 2.1.
5 Proof of Theorem 2.2
The proof is divided into seven steps. Throughout this section, the symbol will denote generic positive constants independent of , and that may change their value from line to line.
Step 1: Implicit time discretisation.
Let be arbitrary and denote a time step size. For , we define functions recursively by the following construction. The zeroth order iterate is taken to be the initial data, i.e., . If is already constructed, we choose to be a minimiser of the functional
[TABLE]
on the set
[TABLE]
As in section 4, if the potentials and grow sufficiently fast, the functional may attain the value , but minimisers, whose existence are addressed in Step 2, will have finite energy. At this moment we briefly describe the construction of time-discrete solutions. Convexity of and allow us to proceed as in [18] to infer that satisfies
[TABLE]
for all , where
[TABLE]
For arbitrary and , we define
[TABLE]
as well as functions
[TABLE]
where the constants are given by
[TABLE]
Then, following the previous section, we infer from (5.2) that the quadruplet , , , satisfies the equations
[TABLE]
for all , , , . This system is an implicit time discretisation of (1.3). Hence, the collection is a time-discrete approximate solution. For and , we define the piecewise constant extension
[TABLE]
for and the piecewise linear extension
[TABLE]
for any and .
Step 2: Existence of a minimiser.
We can use the direct method in the calculus of variations to show that the functional has at least one minimiser in the set . First, we observe that is bounded from below:
[TABLE]
for all . Hence, the infimum exists and we can find a minimising sequence such that
[TABLE]
which leads to the estimates
[TABLE]
for all . As , we can use (5.10) and Poincaré’s inequality to conclude that is a bounded sequence in . Thus, there exists a function such that in as along a nonrelabelled subsequence. This implies that and directly follows.
For the surface variable we have two cases: If we proceed similarly as above and infer the existence of a function with in and . This implies that the limit belongs to . In the case we recall that for all . By (5.10) and the continuous embedding we get
[TABLE]
Furthermore, from (), we can use (5.10) to obtain the estimate
[TABLE]
Hence, there exists a function such that in . The weak convergence implies that and directly follows. This means that .
It remains to prove that is a minimiser of the functional . The convergence in , continuity and nonnegativity of , and Fatou’s lemma yield
[TABLE]
If we can proceed similarly and show that
[TABLE]
In the case we recall that is continuous and convex (hence weakly lower semicontinuous) and that
[TABLE]
Hence, (5.11) holds true thanks to weak lower semicontinuity, and ensures that in all cases we have
[TABLE]
Step 3: Uniform estimates.
In the following, we exploit the gradient flow structure to establish uniform bounds on the piecewise constant extensions, i.e.,
[TABLE]
by following the approach in [19, Sec. 4.4]. Firstly, since is a minimiser of on the set , we infer the a priori estimate
[TABLE]
and subsequently, for all ,
[TABLE]
Since and , we can proceed similarly as in Step 2 and use the definition (5.7) to obtain the uniform bounds
[TABLE]
Next, testing (5.5) with arbitrary such that and using (5.15) gives
[TABLE]
From assumption (\mathrm{A2}$$\mathrm{A2iii}) with and Remark 2.1 we infer that
[TABLE]
Consequently, there exists some nonnegative constant depending on but not on , or , such that
[TABLE]
Proceeding as in [19, Sec. 4.3] we can use a generalised Poincaré inequality (see [1, p. 242]) to conclude that
[TABLE]
To establish a uniform bound on we recall from (5.7) that
[TABLE]
for all , where is fixed. From the definition of and , as well as the a priori estimate (5.13), we derive that
[TABLE]
It follows inductively that
[TABLE]
and in particular, with , we have
[TABLE]
Together with (5.18) we conclude that is uniformly bounded in .
To establish a uniform bound on we test (5.6) with . This gives
[TABLE]
upon using (5.14), an analogous estimate to (5.16) for , the fact that for and
[TABLE]
Thus, is uniformly bounded in , and by Poincaré’s inequality we conclude that is uniformly bounded in .
Step 4: Hölder-in-time estimates.
Via interpolation type arguments we can show that the piecewise linear extension is Hölder continuous in time, i.e.,
[TABLE]
To prove this assertion let be arbitrary where without loss of generality. Since and are continuous and piecewise linear in time, they are weakly differentiable and we infer from (5.3) and (5.4) that
[TABLE]
for all and all test functions and . Choosing , integrating with respect to from to and using the uniform bounds from Step 3, we obtain that
[TABLE]
In the case an analogous estimate can be established for . If , we recall that and choose \theta:=\mathcal{N}_{\Gamma}\big{(}\bar{v}(t_{1})-\bar{v}(t_{2})\big{)}. Proceeding similar as above and using integration by parts yields
[TABLE]
For any , we can find and such that . Then, an immediate consequence of (5.21) and the definitions (5.7) and (5.8) are the uniform estimates
[TABLE]
Moreover, by integrating (5.22) with respect to from [math] to for test functions and , it is easy to see that
[TABLE]
Step 5: Convergence assertions.
Thanks to the uniform estimates (5.12), (5.21), (5.23) and (5.24), we deduce the existence of functions
[TABLE]
such that for any , the following convergences hold:
[TABLE]
and in the case we additionally have and a.e. on . Since the proof follows almost analogously from the arguments in [19, Lem. 8], we omit the details.
Additionally, from the trace inequality in () and (5.12), we find that
[TABLE]
and so it holds that
[TABLE]
Step 6: Existence of weak solutions.
Due to Step 5, the limit possesses the regularity stated in item (i) of Theorem 2.2. To verify item (ii), we pass to the limit in the weak formulation. It is easy to see that (2.8) and (2.9) can be deduced from integrating (5.22) over for test functions and , and passing to the limit . The recovery of terms depending linearly on in (2.10) and (2.11) is straightforward, while the recovery of nonlinear terms involving and (if ) follows analogously as in [19, Sec. 4.7] via Vitali’s convergence theorem. In the following we sketch the details for the nonlinear terms involving and also for in the case .
For the case , for any , let denote a measurable subset of with . In virtue of (), (5.12), Hölder’s inequality and the Sobolev embedding for all , it is easy to see that
[TABLE]
Then Vitali’s convergence theorem ensures that and strongly in , which is sufficient to recover the right-hand sides of (2.10) and (2.11).
For the case , assumption () implies that is a constant, and hence the right-hand sides of (2.10) and (2.11) can be obtained easily. Meanwhile, to recover the term involving on the left-hand side of (2.11), we employ Minty’s trick. By assumption (), the surface potential is convex and nonnegative. We associate its derivative with the operator
[TABLE]
where and . This mapping is well-defined since,
[TABLE]
for all . In particular, since is bounded in , this shows that is uniformly bounded in , and by the Banach–Alaoglu theorem there exists such that in as along a nonrelabelled subsequence. Then, for any arbitrary test function , integrating (5.6) over and passing to the limit leads to
[TABLE]
The idea is to use Minty’s trick to show that . Convexity of implies that the operator is monotone. Furthermore, for arbitrary ,
[TABLE]
which implies that is continuous for arbitrary , i.e., is a hemicontinuous operator. Together with monotonicity, we infer that is maximal monotone. Hence, in order to prove the identification it suffices to show that
[TABLE]
Let be arbitrary. By the monotonicity of we see that
[TABLE]
In light of the convergence results stated in Step 5, as well as
[TABLE]
which follows from (5.23) and weak lower semicontinuity, respectively, we infer from taking the limit superior in (5.28) the inequality
[TABLE]
Using the identity (5.26) yields (5.27) and in particular, weakly in . Thus, we have derived item (ii) of Theorem 2.2.
For item (iii), proceeding similarly to the approach in Step 2 we find that
[TABLE]
for almost all . These results can be used to prove the energy inequality (2.12) by taking the limit inferior in (5.19).
Step 7: Uniqueness.
Let and denote weak two solutions to (1.3) corresponding to the same initial data. We define the difference of these solutions as
[TABLE]
For convenience, we use the notation
[TABLE]
For arbitrary , test functions and , we define , ,
[TABLE]
Then, and . From (2.8) with this particular , we have
[TABLE]
This implies that we can write
[TABLE]
for some constant . Choosing in turn gives
[TABLE]
where we used and . Proceeding similarly, we infer that
[TABLE]
and thus,
[TABLE]
For any , we define the projection as
[TABLE]
satisfying and , where denotes the characteristic function of the set . Then, in (2.10) we can consider the test function . Monotonicity of yields
[TABLE]
By the dominated convergence theorem, we can pass to the limit rigorously and obtain an analogous inequality where is replaced by . Lipschitz continuity of then yields
[TABLE]
For the case , proceeding similarly, in (2.11) we consider the test function and obtain
[TABLE]
Adding these two inequalities to (5.30) leads to
[TABLE]
where we have used the embeddings and for all and the regularities of the solutions. By the Gagliardo–Nirenberg inequality in two dimensions we find that
[TABLE]
Hence, applying Young’s inequality to absorb all terms involving by the left-hand side, and using (2.1), (2.2) and (2.5), we arrive at
[TABLE]
Uniqueness thus follows from Gronwall’s inequality.
For the case , recall that is affine linear and that is convex. Thus, the right-hand side of (5.32) arising from the Lipschitz continuity of is now zero. Adding (5.31) and (5.32) to (5.30) now gives
[TABLE]
owning to the fact and . By the trace inequality in (), Young’s inequality, and arguing as in the above case for terms involving the -norm of , we see that
[TABLE]
Hence, uniqueness follows from Gronwall’s inequality. This completes the proof of Theorem 2.2.
6 Proof of Theorem 2.3
In this section, we take for some and . The letter will denote generic positive constants independent of that may change their value from line to line. Moreover, the symbol will denote the space introduced in (2.7).
Weak convergence.
We first address the convergence of the initial data. From (2.13), () and the Poincaré inequality we see that
[TABLE]
Hence, there exists a nonrelabelled subsequence such that
[TABLE]
The last compactness statement and affine linearity of imply that strongly in and a.e. on , and so by uniqueness of strong and weak limits. In particular, we have that the limit of the initial data satisfies .
Similarly, from the energy inequality (2.12), the assumption (2.13), and Poincaré’s inequality, we immediate infer the following uniform estimates:
[TABLE]
Arguing as in Step 3 of the proof of Theorem 2.2, we use a test function with in the weak formulation (2.10) to obtain, in conjunction with (\mathrm{A2}$$\mathrm{A2iii}) and the previous uniform estimates, that
[TABLE]
Similar to the derivation of (5.18), we use the generalised Poincaré inequality (cf. [1, p. 242]) to infer that
[TABLE]
Consequently, choosing in (2.10) we have
[TABLE]
and together with (2.11) this yields
[TABLE]
As the right-hand side is bounded in we deduce by the Poincaré inequality that
[TABLE]
Lastly, the uniform boundedness of and implies that
[TABLE]
These uniform estimates allows us to deduce that, along a nonrelabelled subsequence, there exists a limit quadruplet such that
[TABLE]
Here we point out that the a.e. convergence of to on in the case is due to the fact that is affine linear and a.e. on . In particular, we have a.e. on and hence a.e. on . Taking arbitrary test functions in (2.8) and in (2.9), integrating over and passing to the limit yields
[TABLE]
We now take an arbitrary test function with in (2.10) and choose in (2.11). Adding these two equalities and integrating over leads to
[TABLE]
For the case , we can follow the arguments in Step 6 of the proof of Theorem 2.2 to pass to the limit in order to obtain
[TABLE]
Hence, substituting and in the above equalities shows that the triplet is a weak solution to (1.2) in the sense of item (ii) of Theorem 2.1.
For the case , we can employ exactly the same argument. Indeed, we have a.e. on and by Vitali’s convergence theorem strongly in . Let us point out that in Step 6 of the proof of Theorem 2.2, we do not have a pointwise a.e convergence of for the case , and hence we have to use Minty’s method to identify the limit of .
Lastly, to recovery the Hölder-in-time regularity of the limit functions, we fix arbitrary with , and consider the (time-independent) test function in (2.8) for . Then integrating over and passing to the limit gives
[TABLE]
which implies that . For , we use the same argument to obtain . Furthermore, for , we take in (2.9) for . Then integrating over and passing to the limit gives
[TABLE]
Hence, fulfills also item (i) of Theorem 2.1.
Error estimates.
We reuse the notation to denote the differences between and , where and . With the regularity assumed on , which is guaranteed by Theorem 2.1 under the condition (2.6), we see that the weak formulation of (1.1) in item (ii) of Theorem 2.1 can be expressed as
[TABLE]
holding for all , , , , and for a.e. . Then, it holds that
[TABLE]
where we have used the relation to obtain the differences in the right-hand sides of (6.3) and (6.4). We note that the hypothesis (2.13) implies and for all . Then, choosing and yields by application of (),
[TABLE]
As in the uniqueness proof of Theorem 2.2, we consider and and employ the monotonicity of and . Then passing to the limit in the resulting inequalities and use the Lipschitz continuity of and leads to
[TABLE]
where for the last inequality we employed the interpolation estimates in (), as well as the identification of the -norm in (). Adding (6.6) to (6.5) leads to the differential inequality for the case
[TABLE]
with a positive constant independent of . Then, the desired estimate (2.14) follows from first applying Gronwall’s inequality to (6.7) and then applying the interpolation inequalities from ().
For the case , the convexity of from () implies that the term involving on the right-hand side of the first inequality in (6.6) is not present. Although the surface gradient term on the left-hand side of (6.6) is lost, nevertheless we can proceed as before to infer the estimate (2.14) without the terms for .
7 Numerical approximation
In this section we present fully discrete finite element approximations of the systems (1.3) and (1.2) based on the implicit time discretisations in Section 5 and Section 4, respectively. Moreover, we display several two-dimensional numerical simulations. Let us mention that there are numerous contributions on the numerical approximation of the Cahn–Hilliard equation with Neumann/periodic boundary conditions (of which we only cite [11, 12, 21, 29, 30] and refer the reader to the references cited therein), and with dynamic boundary conditions [5, 6, 15, 22, 31]. As the focus of this paper is on the analytical aspects, we leave the investigation of efficient numerical algorithms for (1.3) to future research, and a recent contribution on the numerical analysis for the original model of [25] can be found in [26].
7.1 Finite element scheme
Let be a polygonal or a polyhedral approximation of the domain with . The faces of are -simplices with vertices lying on in the sense that is an interpolation of . We consider a quasi-uniform triangulation of consisting of closed simplices with . Associated to is the set of nodal points for and the set of continuous, piecewise linear nodal basis functions satisfying for . Let denote the space of affine linear polynomials on a closed, connected set . Then, we define the finite element space
[TABLE]
As observed in [5, Sec. 2] and [24, Sec. 3.2.1], the set of restrictions of the basis functions forms a basis for the space of piecewise linear polynomials on , which allows us to define
[TABLE]
We further introduce the characteristic function for the boundary nodal points as
[TABLE]
and the Lagrange interpolation operators and as
[TABLE]
For an arbitrary but fixed , let denote the number of time steps and the step size. Based on the implicit time discretisation in Section 5, we propose the following fully discrete finite element scheme: Given and , find and satisfying
[TABLE]
for all and . In the above, and denote the -inner product on and on , respectively. Meanwhile, and are the lumped integrations defined as
[TABLE]
Furthermore, for we set where is the -projection onto and analogously, where is the -projection onto . The use of lumped integration in (7.1) leads to diagonal mass matrices:
[TABLE]
Together with the symmetric stiffness matrices
[TABLE]
and the vectors of nonlinearities (using the Einstein summation convention in the arguments)
[TABLE]
for and we find that (7.1) with and leads to the following system of nonlinear equations to solve
[TABLE]
for our finite element solutions
[TABLE]
For the case , we have
[TABLE]
where is the vector with all entries equal to one, and the discrete nonlinear system (7.2) simplifies considerably.
In view of the choice of our finite element basis functions, we can derive a finite element scheme for the limit system (1.2) based on the implicit time discretisation (4.1): Given , find and satisfying
[TABLE]
for all and . Introducing the vector
[TABLE]
and recalling the matrices , , and , as well as the vector , we find that (7.3) with and leads to the following system of nonlinear equations to solve
[TABLE]
for our finite element solutions
[TABLE]
The subsequent nonlinear schemes (7.2) and (7.4) are solved using a Newton method with MATLAB built on previous work by D. Trautwein [31]. In the following, we always take with a standard Friedrichs–Keller triangulation and spatial step size . This yields 101 degrees of freedom along each axis.
7.2 Qualitative experiments
7.2.1 Dynamics of the limit system (1.1)
In the first set of experiment, we demonstrate the difference in dynamics for the limit system (1.1) for the case , i.e., the original model proposed in [25], and for nontrivial values of . We begin with the initial condition
[TABLE]
for . It is clear that and .
Effects of varying with .
We then solve (7.4) with the fixed parameters and potentials
[TABLE]
In Figure 1 we plot the projection of the discrete solution on the line at iteration for , as the solution seems to be constant in the orthogonal direction, as well as the differences between the profiles for and , and for and . Identical profiles are observed for the case of negative , and so we do not report them here.
It appears that the overall solution profile is rather robust with respect to the size of . We also report that the magnitude of the difference between the profiles for and is of order , while the magnitude of the difference between the profiles for and is of order . Further comparison between and yields a difference of order , which seems to suggest some form of convergence as .
Effects of varying with .
Next, we take the same parameter values for , and the same potentials as in (7.5), but now with and . In Figure 2 we plot the discrete solution at iteration for and . Similar profiles are observed for the other cases where and are positive/negative, and so we do not report them here. We mention that in preliminary testings the profiles for are similar to their positive counterparts but reflected along the line .
Compared to Figure 1, it appears that the parameter has a stronger influence on the solution profile than . Most interestingly, as increases in magnitude, the trace of the discrete solution on the domain boundary approaches zero, while the range in the bulk domain decreases, as seen in the bottom row of Figure 2. Whereas for , the solution profile exhibits some monotonicity along the -direction. The transition between the two different profiles occur approximately around .
Let us attempt an explanation of the numerical observations. The similarity of solution profiles for positive/negative values of and can be attributed to the symmetry of the potential and the initial condition . In choosing we had hoped that the boundary dynamics might dominate the bulk dynamics, and prioritise driving the trace of to the stable minima of , which in this case are . For the case and this is inconsistent with the conservation of boundary mass, and so we believe that the system instead favours to have the trace of to be close to zero. In turn the solution has to accommodate by moving away from the monotonic profile inherited from the initial condition to that observed in Figure 2. In particular, for our choice of parameter values, potentials and initial condition, it appears that the mass conservation mechanism imposes a stronger influence than the Cahn–Hilliard phase separation mechanism.
This competition between the two mechanism is more pronounced by considering an asymmetric initial condition , which we take with values uniformly distributed in such that and to three decimal points. We now consider the parameters and potentials to be
[TABLE]
The values of and are chosen to allow for the typical Cahn–Hilliard phase separation dynamics, which is observed in Figure 3, where we plot the discrete solution at iteration for and selected values of .
We point out two observations from Figure 3. The first observation is, unlike Figure 2, there is no symmetry corresponding to the case of positive and negative values of , as evident in the second and third columns. The second observation is that for , the trace of the discrete solution is nearly constant on the domain boundary, with an average value of 0.401 to three decimal points. In particular, this is in good agreement with the conservation of boundary mass ( to three decimal points), and is similar to the behaviour in the previous set of experiments. On the other hand, for small there is a remarkable difference in solution behaviour. In the second column of Figure 3 we have the plots for and , in which the Cahn–Hilliard phase separation dynamics can be clearly seen on the domain boundary. The black level lines in each of these two plots correspond to the stable minima of the surface potential , which are for and for . However, for and , we see that the trace of the discrete solution takes nearly constant values, similarly as in the case for large values. We further comment that the solution behaviour is insensitive to the sign of , since the stable minima are regardless of the sign of .
Based on Figure 3, if and denote the stable minima of with , we could postulate that Cahn–Hilliard phase separation is more likely to occur on the domain boundary if . Otherwise, it is likely that the trace of the solution is spatially constant with values close to .
However, in some simulations not reported here we have observed counterexamples to this conjecture. In particular, for , the trace of the discrete solution is nearly constant with values close to , even though in each of these cases. Therefore, at present we do not have a robust method to predict whether Cahn–Hilliard phase separation or near constant trace values will occur based on knowledge of the parameters and the initial data.
7.2.2 Dynamics of the Robin system (1.3)
Convergence as .
We now demonstrate numerical evidence in support of Theorem 2.3. The potentials are as described in (7.6), and for the following choice of parameters:
[TABLE]
In this section we report on the errors between the discrete solutions to the limit model (1.1) and the Robin model (1.3) as decreases. For fixed , let denote the error under consideration at , . Then, the experimental order of convergence is defined as
[TABLE]
For and , the -norm of a function is approximated as
[TABLE]
where is the discrete approximation of at iteration .
To demonstrate the convergence behaviour for , we consider two different discrete initial data as illustrated in Figure 4:
- (a)
For any node , we set . 2. (b)
For any node , attains a random value in . These random values are uniformly distributed.
In view of Theorem 2.3, we set . In particular, the discrete initial data do not depend on . In Figures 5 and 6, we display the discrete solution to initial data (a) and (b), respectively, at iteration for and three choices of along with the corresponding discrete solution to the limit system (1.2). Visually, we observe good agreement of the discrete solution at and the discrete solution to the limit system.
In Tables 1 to 5, we collect the error between the discrete solution to (1.3) and the discrete solution to (1.2) (where stands for ) measured in various norms for several choices of and for the initial data (a) and (b). In the tables, the labelling of the columns is to be understood as follows:
- •
The columns and display the error ,
- •
the column displays the difference ,
- •
the columns and display the error
measured in the respective norms on these spaces.
In all instances, a linear rate of convergence is observed. In view of Theorem 2.3, this seems to suggest that the error estimate for in is sharp, while the rate of convergence for the other terms in (1.4) can be improved further.
Simulations for nonlinear transmission relations.
In the last set of experiments we consider nonlinear relations for the system (1.3). We additionally assume that there exists an interval such that is a bijection. For the simulations (see Figure 7) we use once more the initial data (a) and (b) (as defined above) where we set . As we did not derive a limit system in the general nonlinear case, we use the solution to system (1.3) with as the reference solution instead. The discrete error measured in various norms (as described above) between and are displayed in Tables 6 and 7, where we again observe a linear rate of convergence.
8 Conclusion
In this work, we derived and analysed an extension of a recent Cahn–Hilliard model proposed by C. Liu and H. Wu with a new type of dynamic boundary conditions. The extension replaces the classical transmission condition with an affine linear relation between the bulk and surface order parameters. A natural time discretisation inspired by the associated gradient flow structure allow us to establish the existence of unique weak solutions, and numerical simulations demonstrate that under certain choices of potentials and initial conditions, there seems to be a competition between mass conservation and Cahn–Hilliard phase separation dynamics occurring on the domain boundary. This behaviour is not observed with the original model of Liu and Wu, and at present we have yet to find a criterion based on the initial data and parameter values that would allow us to predict which of the competing effects is dominant.
As a first step towards nonlinear relations in the dynamic boundary condition, we study a Robin boundary type regularisation. Similarly, the regularised system exhibits a gradient flow structure allows us to exploit a natural time discretisation to establish the existence of unique weak solutions. In comparison, for the case where there is no diffusion for the surface phase field variable, we have to place more restrictions on the surface potential to obtain analogous results. Meanwhile, for the case of affine linear relations, we can show the weak convergence of solutions to the Robin problem as the regularisation parameter tends to zero, and recover the unique solution to the extended Liu–Wu model. An error estimate is also derived and is supported by our numerical findings. Furthermore, our simulations also suggest that similar rates of convergence hold for nonlinear relations in the dynamic boundary condition, and we leave the analytical investigation to future research.
Acknowledgements.
The work of the second author was supported by a Direct Grant of CUHK, Hong Kong (project 4053288). The kind support of the Fudan Scholar program from Fudan University is also acknowledged. The authors would like to thank the anonymous referees for their careful reading and useful suggestions which have improved the quality of the manuscript.
Appendix A Energetical variational approach for derivation of (1.2)
The derivation of (1.2) using the energetically variational approach proceeds in three steps.
Mass conservation.
In the bulk domain , we assume there is a microscopic effective velocity satisfying the no-flux boundary condition on , and is a locally conserved (sufficiently regular) quantity satisfying the continuity equation
[TABLE]
Furthermore, we assume that is a microscopic effective tangential velocity field on and the boundary dynamics are characterised by a local mass conservation law for the trace of , which we denote again as :
[TABLE]
Then it is clear that and .
Energy dissipation.
The basic energy law considered is
[TABLE]
where and is the rate of energy dissipation defined as
[TABLE]
Force balance.
It remains to determine the microscopic velocities and in the continuity equations (A.1) and (A.2) via the energetic variational approach. We denote by the reference configuration with Lagrangian coordinate , and the deformed configuration with Eulerian coordinate , which are related via the bulk flow map defined as a solution to the ordinary differential system
[TABLE]
induced by the bulk microscopic velocity . Similarly, we introduce the reference boundary , the deformed boundary and the surface flow map induced by with Lagrangian and Eulerian surface coordinates and , respectively. We mention that the surface flow map is the solution to an analogous ordinary differential system with the microscopic velocity , and is not related to the bulk flow map. Now, let the total action functional be defined as the negative integral of over . By computations similar to [25, Sec. 7.1], the variation of with respect to the spatial variables in and on yields
[TABLE]
for arbitrary smooth vector fields and satisfying on . Here, and are defined as
[TABLE]
Invoking the principle of least action, the conservative force balance equations are identified by solving , which yields the bulk and surface conservative forces:
[TABLE]
Furthermore, via Onsager’s maximum dissipation principle for irreversible dissipative processes in the regime of linear response, taking the variation of the Rayleigh dissipation function with respect to the rate functions (i.e., the microscopic velocities and ) gives the dissipative force balance equations. In our case, they read as
[TABLE]
The classical Newton’s force balance law leads to the force balance relations in the bulk and on the boundary , which additionally allows us to infer an expression for the microscopic velocities:
[TABLE]
Eventually, in combination with (A.1) and (A.2) we conclude that
[TABLE]
which coincides with the system (1.2).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] H.W. Alt. Linear Functional Analysis - An Application-Oriented Introduction . Springer-Verlag London, 2016.
- 2[2] A. Buffa, M. Costabel and D. Sheen. On traces for H(curl, Ω Ω \Omega ) in Lipschitz domains. J. Math. Anal. Appl. 276 (2002) 845–867
- 3[3] L. Calatroni and P. Colli. Global solution to the Allen-Cahn equation with singular potentials and dynamic boundary conditions. Nonlinear Anal. 79 (2013) 12–27
- 4[4] L. Cherfils, S. Gatti and A. Miranville. A variational approach to a Cahn–Hilliard model in a domain with nonpermeable walls. J. Math. Sci. 189 (2013) 604–636
- 5[5] L. Cherfils and M. Petcu. A numerical analysis of the Cahn–Hilliard equation with non-permeable walls. Numer. Math. 128 (2014) 517–549
- 6[6] L. Cherfils, M. Petcu and M. Pierre. A numerical analysis of the Cahn–Hilliard equation with dynamic boundary conditions. Discrete Contin. Dyn. Syst. 27 (2010) 1511–1533
- 7[7] P. Colli and T. Fukao. Equation and dynamic boundary condition of Cahn–Hilliard type with singular potentials. Nonlinear Anal. 127 (2015) 413–433
- 8[8] P. Colli, T. Fukao and K.F. Lam. On a coupled bulk-surface Allen–Cahn system with an affine linear transmission condition and its approximation by a Robin boundary condition. Nonlinear Anal. 184 (2019) 116–147