Discontinuous Galerkin approximation of linear parabolic problems with dynamic boundary conditions

0.42 MB·English

by  Paola F. Antonietti

Discontinuous Galerkin approximation of linear
parabolic problems with dynamic boundary
5
1
conditions
0
2

n P.F. Antonietti ∗ M. Grasselli † S. Stangalino‡ M. Verani§
a
J
January 21, 2015
0
2

]
A
N Abstract
.
h In this paper we propose and analyze a Discontinuous Galerkin
t method for a linear parabolic problem with dynamic boundary con-
a
m ditions. We present the formulation and prove stability and optimal
a priori error estimates for the fully discrete scheme. More precisely,
[
using polynomials of degree p≥1 on meshes with granularity h along

1 with a backward Euler time-stepping scheme with time-step ∆t, we
v
prove that the fully-discrete solution is bounded by the data and it
5
converges, in a suitable (mesh-dependent) energy norm, to the exact
6
7 solution with optimal order hp+∆t. The sharpness of the theoretical
4 estimates are verified through several numerical experiments.
0
.
1
1 Introduction
0
5
1 InthispaperwepresentandanalyzeaDiscontinuousGalerkin(DG)method
:
v forthefollowinglinearparabolicproblemsupplementedwithdynamicbound-
i
X
∗MOX-Dipartimento di Matematica, Politecnico di Milano, P.zza Leonardo Da Vinci
r 32, I-20133 Milano, Italy ([email protected]).
a
†Dipartimento di Matematica, Politecnico di Milano, P.zza Leonardo Da Vinci 32,
I-20133 Milano, Italy ([email protected]).
‡MOX-Dipartimento di Matematica, Politecnico di Milano, P.zza Leonardo Da Vinci
32, I-20133 Milano, Italy ([email protected]).
§MOX-Dipartimento di Matematica, Politecnico di Milano, P.zza Leonardo Da Vinci
32, I-20133 Milano, Italy ([email protected]).
1
ary conditions on Γ :
1

∂ u = ∆u+f, in Ω, 0 < t ≤ T,
 t

 ∂ u = −αu+β∆ u−λ∂ u+g, on Γ , 0 < t ≤ T,
n Γ t 1 (1)
periodic boundary conditions, on Γ , 0 < t ≤ T,
 2

 u = u , in Ω.
|t=0 0
Here the domain Ω and the subsets Γ ⊂ ∂Ω, i = 1,2, are depicted in
i
Figure1,∆ istheLaplace-Beltramioperator,∂ udenotestheouternormal
Γ n
derivative of u on Γ , g is a given function and α,β,λ are suitable non-
1
negative constants.
Dynamic boundary conditions have been recently considered by physi-
cists to model the fluid interactions with the domain’s walls (see, e.g.,
[11, 12, 19]). Despite the practical relevance of this kind of boundary con-
ditions from a modeling point of view and the intense research activity to
understand their analytical properties, see, e.g., [15, 29, 30], the study of
suitable numerical methods for their discretization is still in its infancy. To
the best of our knowledge, the only work along this direction is [5], where
the authors analyze a conforming finite element method for the approxima-
tion of the Cahn-Hilliard equation supplemented with dynamic boundary
conditions. Motivated by the flexibility and versatility of DG methods, here
weproposeandanalyzeaDGmethodcombinedwithabackwardEulertime
advancing scheme for the discretization of a linear parabolic problem with
dynamic boundary conditions. The main goal of the present work is the
numerical treatment of dynamic boundary conditions within the DG frame-
work. Here we consider just a linear equation. However, our results aim to
be a key step towards the extension to (non-linear) partial differential equa-
tions with dynamic boundary conditions, as, for example, the Cahn-Hilliard
equation. Inthiscontext,wementionDGmethodshavebeenalreadyproved
to be an effective discretization strategy for the Cahn-Hilliard equation as
shown in [18] where the authors constructed and analyzed a DG method
coupled with a backward Euler time-stepping scheme for a Cahn-Hilliard
equation in two-dimensions, cf. also [31].
The origins of DG methods can be backtracked to [24, 21] where they
have introduced for the discretization of the neutron transport equation.
Since that time, DG methods for the numerical solution of partial differen-
tialequationshaveenjoyedagreatdevelopment,seethemonographs[25,14]
for an overview, and [3] for a unified analysis of DG methods for elliptic
problems. In the context of parabolic equations, DG methods in primal
form combined with backward Euler and Crank-Nicholson time advancing
2
techniques have been firstly analyzed in [2, 26], respectively. DG in time
methods have also been studied for parabolic partial differential equations,
see, for example, [16, 8, 9, 20] and the reference therein; cf. also [28, 27] for
the hp-version of the DG time-stepping method.
The paper is organized as follows. In Section 2 we introduce some useful
notation and the functional setting. Section 3 is devoted to the introduction
and analysis of a DG method for a suitable auxiliary (stationary) problem.
These results will be then employed in Section 4 to design a DG scheme to
approximate the linear parabolic problem with boundary conditions and to
obtain optimal a priori error estimates for the fully discrete scheme. Finally,
in Section 6 we numerically assess the validity of our theoretical analysis.
2 Notation and functional setting
In this section we introduce some notation and the functional setting.
Let D ⊂ R2 be an open, bounded, polygonal domain with boundary
Γ = ∂D. On D we define the standard Sobolev space Hs(D), s = 0,1,2,...
(for s = 0 we write L2(D) instead of H0(D)) and endow it with the usual
inner scalar product (·,·) , and its induced norm (cid:107)·(cid:107) , cf. [1]. We
Hs(D) Hs(D)
also need the seminorm defined by |·| = ((cid:80) (cid:107)∂α(·)(cid:107) )1/2.
Hs(D) |α|=s L2(Ω)
We next introduce, on Γ, the Laplace-Beltrami operator. We first define
the projection matrix P = I − n ⊗ n = (δ − n n )2 , where n is the
ij i j i,j=1
outward unit normal to D, a⊗b = (a b ) is the dyadic product, and δ is
i j ij ij
the Kroneker delta. We define the tangential gradient of a (regular enough)
scalar function u : Γ → R as ∇ u = P∇u. The tangential divergence of
Γ
a vector-valued function A : Γ → R2 is defined as div (A) = Tr(cid:0)(∇A)P(cid:1),
Γ
being Tr(·) the trace operator. With the above notation, we define the
Laplace-Beltrami operator as ∆ u = div (∇ u).
Γ Γ Γ
We next introduce the following Sobolev surface space
Hs(Γ) = {v ∈ Hs−1(Γ) | ∇ v ∈ [Hs−1(Γ)]2}, s ≥ 1,
Γ
cf. [7], with the convention that H0(Γ) ≡ L2(Γ), L2(Γ) being the standard
Sobolev space of square integrable functions (equipped with the usual inner
scalar product (·,·) and the usual induced norm (cid:107)·(cid:107) ). We equipped
Γ L2(Γ)
3
the space Hs(Γ) with the following surface seminorm and norm
|v| = (cid:107)∇ v(cid:107) ∀v ∈ Hs(Γ), s ≥ 1,
Hs(Γ) Γ Hs−1(Γ)
(cid:113)
(cid:107)v(cid:107) = (cid:107)v(cid:107)2 +|v|2 ∀v ∈ Hs(Γ), s ≥ 1,
Hs(Γ) Hs−1(Γ) Hs(Γ)
respectively. In [17, Lemma 2.4] is proved that the above norm is equivalent
to the usual surface norm present in literature [22], which is defined in local
coordinates after a truncation by a partition of unity.
Next, for a positive constant λ, we introduce the space
Hs(D,Γ) = {v ∈ Hs(D) : λv ∈ Hs(Γ)}, s ≥ 0,
λ |Γ
and endow it with the norm
(cid:114)
(cid:16) (cid:17)
(cid:107)u(cid:107) = (cid:107)u(cid:107)2 +λ(cid:107)u (cid:107)2 .
Hλs(D,Γ) Hs(D) |Γ Hs(Γ)
As before, for s = 0 we will write Hs(D,Γ) instead of H0(D,Γ). Moreover,
λ λ
to ease the notation, when λ = 1, we will omit the subscript.
Finally, throughout the paper, we will write x (cid:46) y to signify x ≤ Cy,
where C is a generic positive constant whose value, possibly different at any
occurrence, does not depend on the discretization parameters.
3 The stationary problem and its DG discretiza-
tion
Let Ω = (a,b)×(c,d) ⊂ R2 be a rectangular domain and let Γ ,Γ be the
1 2
union of the top and bottom/left and right edges, respectively, cf. Figure 1.
WeconsiderthefollowingLaplaceproblemwithgeneralizedRobinboundary
conditions:

−∆u = f, in Ω,

∂ u = −αu+β∆ u+g, on Γ , (2)
n Γ 1
 periodic boundary conditions, on Γ ,
2
where α,β are positive constants, and f ∈ L2(Ω), g ∈ L2(Γ ) are given
1
functions.
Defining the bilinear form a(u,v) : H1(Ω,Γ )×H1(Ω,Γ ) → R as
1 1
a(u,v) = (∇u,∇v) +β(∇ u,∇ v) +α(u,v) ,
L2(Ω) Γ Γ L2(Γ1) L2(Γ1)
4
the weak formulation of (2) reads: find u ∈ H1(Ω,Γ ) such that
1
a(u,v) = (f,v) +(g,v) ∀v ∈ H1(Ω,Γ ). (3)
L2(Ω) L2(Γ1) 1
The following result shows that formulation (3) is well posed.
Theorem 3.1. Problem (3) admits a unique solution u ∈ H2(Ω,Γ ) satis-
1
fying the following stability bound:
(cid:107)u(cid:107) (cid:46) (cid:107)f(cid:107) +(cid:107)g(cid:107) . (4)
H2(Ω,Γ1) L2(Ω) L2(Γ1)
Moreover, if f ∈ Hs−2(Ω) and g ∈ Hs−2(Γ ), s ≥ 2, then u ∈ Hs(Ω,Γ )
1 1
and
(cid:107)u(cid:107) (cid:46) (cid:107)f(cid:107) +(cid:107)g(cid:107) . (5)
Hs(Ω,Γ1) Hs−2(Ω) Hs−2(Γ1)
Proof. The existence and uniqueness of the solution are proved in [17, The-
orem 3.2]. The proof of the regularity results is shown in [17, Theorem
3.3-3.4]. The same arguments used in [17, Theorem 3.3-3.4] apply also in
our case thanks to periodic conditions.
Remark3.2. Weobservethattheforthcominganalysisholdsinmoregeneral-
shaped domains and/or more general type of boundary conditions provided
that the exact solution of the differential problem analogous to (3) satisfies
a stability bound of the form of (4).
3.1 Discontinuous Galerkin space discretization
In this Section we present a discontinuous Galerkin (DG) approximation of
problem (3).
Let T be a quasi-uniform partition of Ω into disjoint open triangles T such
h
that Ω = ∪ T. We set h = max{diam(T), T ∈ T }. For s ≥ 0, we define
T∈Th h
the following broken space
Hs(T ) = {v ∈ L2(Ω) : v ∈ Hs(T,∂T), T ∈ T },
h |T h
where, as before, H0(T ) = L2(T ). For an integer p ≥ 1, we also define the
h h
finite dimensional space
Vp(T ) = {v ∈ L2(Ω) : v ∈ Pp(T), T ∈ T } ⊂ Hs(T ),
h |T h h
for any s ≥ 0. An interior edge e is defined as the non-empty intersection
of the closure of two neighboring elements, i.e., e = T ∩T , for T ,T ∈ T .
1 2 1 2 h
We collect all the interior edges in the set E0. Recalling that on Γ ⊂ ∂Ω we
h 2
5
impose periodic boundary conditions, we decompose Γ as Γ = Γ+ ∪Γ−,
2 2 2 2
cf. Figure 1 (left), and identify Γ+ with Γ−, cf. Figure 1 (right). Then
2 2
we define the set EΓ2 of the periodic boundary edges as follows. An edge
h
e ∈ EΓ2 if e = ∂T−∩∂T+, where T± ∈ T such that ∂T± ⊆ Γ±, cf. Figure 1
h h 2
(right). We also define a boundary edge e as the non-empty intersection
Γ1
between the closure of an element in T and Γ and the set of those edges
h 1
by EΓ1. Finally, we define a boundary ridge r as the subset of the mesh
h
vertexes that lie on Γ , and collect all the ridges r in the set RΓ1. Clearly,
1 h
the corner ridges have to be identified according to the periodic boundary
conditions (cf. Figure 1, right). The set of all edge will be denoted by E ,
h
i.e., E = E0∪EΓ1 ∪EΓ2.
h h h h
Figure 1: Example of a domain Ω and an admissible triangulation T (left).
h
On the right, we highlight the edges e ∈ EΓ2 with red lines.
h
For v ∈ Hs(T ), s ≥ 1, we define
h
(cid:88) (cid:88)
|v|2 = |v|2 , |v|2 = |v|2 .
Hs(Th) T∈Th Hs(T) Hs(EhΓ1) eΓ1∈EhΓ1 Hs(eΓ1)
Next, for each e ∈ E0 ∪EΓ2 we define the jumps and the averages of v ∈
h h
H1(T ) as
h
1
[v] = (v+)n++(v−)n− and {v} = (v++v−),
e e e e 2
where v± = v and n± is the unit normal vector to e pointing outward of
|T± e
T±. For each e ∈ EΓ1 we define
h
[v] = v n , {v} = v , v ∈ H1(T ).
e |e e e |e h
Analogously, for each r ∈ RΓ1, we set
h
1
[v] = (v+(r))n++(v−(r))n− and {v} = (v+(r)+v−(r)),
r r r r 2
6
where, denoting by e± the two edges sharing the ridge r, v±(r) = v (r)
|e±
and n± is the unit tangent vector to Γ on r pointing outward of e±. The
r 1
abovedefinitionscanbeimmediatelyextendedtoa(regularenough)vector-
valued function, cf. [3]. To simplify the notation, when the meaning will be
clear from the context, we remove the subscripts from the jump and average
operators. Adopting the convention that
(cid:88) (cid:88)
(v,w) = (v,w) , (ξ,η) = ξ(r)η(r)
Eh L2(e) RΓ1
h
e∈Eh r∈RΓ1
h
for regular enough functions v,w,ξ,η, we introduce the following bilinear
forms
(cid:88)
B (v,w) = (∇v,∇w) −([v],{∇w}) −([w],{∇v}) +σ([v],[w])
h T E0 E0 E0
h h h
T∈T
h
−([v],{∇w}) −([w],{∇v}) +σ([v],[w])
EΓ2 EΓ2 EΓ2
h h h
and
b (v,w) = (∇ v,∇ w) −([v],{∇ w}) −([w],{∇ v}) +σ([v],[w]) ,
h Γ Γ EΓ1 Γ RΓ1 Γ RΓ1 RΓ1
h h h h
for all v,w ∈ H2(T ). Here σ = γ, being γ a positive constant at our
h h
disposal. We then set
A (u,v) = B (u,v)+α (u,v) +β b (u,v). (6)
h h L2(Γ1) h
The discontinuous Galerkin approximation of problem (2) reads: find
u ∈ Vp(T ) such that
h h
A (u ,v ) = (f,v ) +(g,v ) ∀v ∈ Vp(T ). (7)
h h h h L2(Ω) h L2(Γ1) h h
In the following we show that the bilinear form A (·,·) is continuous and
h
coercive in a suitable (mesh-dependent) energy norm. To this aim, for w ∈
Hs(T ), we define the seminorm
h
1
|||w|||2 = |w|2 +σ(cid:107)[w](cid:107)2 + (cid:107){∇w}(cid:107)2
Bh H1(Th) L2(E0∪EΓ2) σ L2(E0∪EΓ2)
h h h h
and the norm
|||w|||2 = |||w|||2 +α(cid:107)w(cid:107)2
∗ Bh L2(Γ1)
β
+β|w|2 +βσ(cid:107)[w](cid:107)2 + (cid:107){∇ w}(cid:107)2 , (8)
H1(EΓ1) L2(RΓ1) σ Γ L2(RΓ1)
h h h
7
where we adopted the notation
(cid:88) (cid:88)
(cid:107)w(cid:107)2 = (cid:107)w(cid:107)2 , (cid:107)w(cid:107)2 = (cid:107)w(cid:107)2 .
L2(Eh) L2(e) L2(RΓ1) L2(r)
e∈Eh h e∈RΓ1
h
Reasoning as in [2], it is easy to prove the following result.
Lemma 3.3. It holds
A (v,w) (cid:46) |||v||| |||w||| ∀v,w ∈ H2(T ). (9)
h ∗ ∗ h
Moreover, for γ large enough, it holds
|||v|||2 (cid:46) A (v,v) ∀v ∈ Vp(T ). (10)
∗ h h
Proof. Let us first prove (9). The term B (·,·) can be bounded by Cauchy-
h
Schwarz inequality as in [2]. Also the term b (·,·) can be handled using the
h
Cauchy-Schwarz inequality:
(cid:12)
|b (v,w)| = (cid:12)(∇ v,∇ w) −([v],{∇ w})
h (cid:12) Γ Γ EΓ1 Γ RΓ1
h h
(cid:12)
−([w],{∇ v}) +σ([v],[w]) (cid:12)
Γ RΓ1 RΓ1(cid:12)
h h
(cid:16) 1 (cid:17)1/2
(cid:46) |v|2 +σ(cid:107)[v](cid:107)2 + (cid:107){∇ v}(cid:107)2 ×
H1(EΓ1) L2(RΓ1) σ Γ L2(RΓ1)
h h h
(cid:16) 1 (cid:17)1/2
|w|2 +σ(cid:107)[w](cid:107)2 + (cid:107){∇ w}(cid:107)2 ,
H1(EΓ1) L2(RΓ1) σ Γ L2(RΓ1)
h h h
and (9) follows employing the definition (8) of the norm |||·||| .
∗
We now prove (10). As before the term B (·,·) can be bounded as in [2]:
h
using the classical polynomial inverse inequality [6] we obtain
|||v|||2 (cid:46) |v|2 +σ(cid:107)[v](cid:107)2 (cid:46) B (v,v)
Bh H1(Th) L2(E0∪EΓ2) h
h h
for all v ∈ Vp(T ). The term b (·,·) can be estimated as follows:
h h
(cid:12) (cid:12)
b (v,v) ≥ |v|2 −2(cid:12)([v],{∇ v}) (cid:12)+σ(cid:107)[v](cid:107)2 .
h H1(EΓ1) (cid:12) Γ RΓ1(cid:12) L2(RΓ1)
h h h
Employing the arithmetic-geometric inequality we get:
(cid:12) (cid:12)
(cid:12)([v],{∇ v}) (cid:12) ≤ (cid:107)σ1/2[v](cid:107) (cid:107){σ−1/2∇ v}(cid:107)
(cid:12) Γ RΓ1(cid:12) L2(RΓ1) Γ L2(RΓ1)
h h h
1
≤ σ(cid:107)[v](cid:107)2 +4(cid:15)σ−1(cid:107){∇ v}(cid:107)2 ,
(cid:15) L2(RΓ1) Γ L2(RΓ1)
h h
8
for a positive (cid:15) > 0. Finally, estimate (10) follows using the polynomial
inverse inequality
h(cid:107){∇ v}(cid:107)2 (cid:46) |v|2 ∀v ∈ Vp(T )
Γ L2(RΓ1) H1(EΓ1) h
h h
and choosing γ sufficiently large.
Thefollowingresultshowsthatproblem(7)admitsauniquesolutionand
that the Galerkin orthogonality property is satisfied. The proof is straight-
forward and we omit it for sake of brevity.
Lemma 3.4. Assume that γ is sufficiently large. Then, the discrete solution
u of problem (7) exists and is unique. Moreover, formulation (7) is strongly
h
consistent, i.e.,
A (u−u ,v) = 0 ∀v ∈ Vp(T ). (11)
h h h
Forv ∈ Hs(Ω,Γ ),s ≥ 2,letIhvbethepiecewiseLagrangianinterpolant
1 p
of order p of u on T . Note that (Ihu) interpolates u on the set of degrees
h p |Γ1
of freedom that lie on EΓ1. By standard approximation results we get the
h
following interpolation estimate.
Lemma 3.5. For all v ∈ Hs(Ω,Γ ), s ≥ 2, it holds
1
|||v−Ihv||| (cid:46) hmin(s−1,p)(cid:107)v(cid:107) .
p ∗ Hs(Ω,Γ1)
Proof. Using the definition (8) of |||·||| norm and that Ihv(r) = v(r) for all
∗ p
r ∈ RΓ1, we get
h
|||v−Ihv|||2 = |||v−Ihv|||2 +α(cid:107)v−Ihv(cid:107)2 +β|v−Ihv|2 . (12)
p ∗ p Bh p L2(Γ1) p H1(EhΓ1)
Expanding the first term at right-hand side and using the multiplicative
trace inequalities
(cid:107)v(cid:107)2 (cid:46) h−1(cid:107)v(cid:107)2 +h|v|2 ,
L2(E ) L2(Ω) H1(Ω)
h
(cid:107)∇v(cid:107)2 (cid:46) h−1|v|2 +h|v|2 ,
L2(E ) H1(Ω) H2(Ω)
h
cf. [25], we get
|||v−Ihv|||2 = |v−Ihv|2 +σ(cid:107)[v−Ihv](cid:107)2
p Bh p H1(Ω) p L2(Eh0∪EhΓ2)
1
+ (cid:107){∇(v−Ihv)}(cid:107)2
σ p L2(E0∪EΓ2)
h h
(cid:46) h−2(cid:107)v−Ihv(cid:107)2 +|v−Ihv|2 +h2|v−Ihv|2 .
p L2(Ω) p H1(Ω) p H2(Ω)
Using standard interpolation estimates [23] we get the thesis.
9
Now we show that the discrete solution u of (7) converges to the weak
h
solution of (3).
Theorem 3.6. Let u ∈ Hs(Ω,Γ ), s ≥ 2, be the solution of the problem (3)
1
and let u be the solution of the problem (7). Then,
h
(cid:107)u−u (cid:107) +h|||u−u ||| (cid:46) hmin(s,p+1)(cid:107)u(cid:107) ,
h L2(Ω,Γ1) h ∗ Hs(Ω,Γ1)
provided γ is chosen sufficiently large.
Proof. By the triangular inequality we have
|||u−u ||| ≤ |||u−Ihu||| +|||Ihu−u ||| .
h ∗ p ∗ p h ∗
We first bound the second term on the right-hand side. Combining the
Galerkin orthogonality (11) with the continuity and the coervicity estimates
(9)-(10), we obtain:
|||Ihu−u |||2 (cid:46) A (Ihu−u ,Ihu−u )
p h ∗ h p h p h
= A (Ihu−u,Ihu−u )+A (u−u ,Ihu−u )
h p p h h h p h
(cid:46) |||Ihu−u ||| |||Ihu−u||| .
p h ∗ p ∗
Therefore,
|||Ihu−u ||| (cid:46) |||Ihu−u||| ,
p h ∗ p ∗
and
|||u−u ||| (cid:46) |||u−Ihu||| .
h ∗ p ∗
Then, using Lemma 3.5, we get
|||u−u ||| (cid:46) hmin(s−1,p)(cid:107)u(cid:107) . (13)
h ∗ Hs(Ω,Γ1)
For the L2 error estimate, we consider the following adjoint problem: find ζ
such that
(cid:26)
−∆ζ = u−u , in Ω,
h
∂ ζ = −αζ +β∆ ζ +(u−u ), on Γ ,
n Γ h 1
Asu−u ∈ L2(Ω,Γ ),usingTheorem3.1yieldsanuniqueζ ∈ H2(Ω,Γ )
h 1 1
satisfiying the following stability estimate
(cid:107)ζ(cid:107) (cid:46) (cid:107)u−u (cid:107) .
H2(Ω,Γ1) h L2(Ω,Γ1)
10