• No results found

Coupling requirements for multi-physics problems

N/A
N/A
Protected

Academic year: 2021

Share "Coupling requirements for multi-physics problems"

Copied!
21
0
0

Loading.... (view fulltext now)

Full text

(1)

Coupling Requirements for

Multi-physics Problems

Fatemeh Ghasemi and Jan Nordstr¨om

LiTH-MAT-R--2016/08--SE

(2)

Link¨

oping University

S-581 83 Link¨

oping

(3)

FATEMEH GHASEMI∗ AND JAN NORDSTR ¨OM

Abstract. We consider two hyperbolic systems on first order form of different size posed on two domains. Our ambition is to derive general conditions for when the two systems can and cannot be coupled.

The adjoint equations are derived and well-posedness of the primal and dual problems are dis-cussed. By applying the energy method, interface conditions for the primal and dual problems are derived such that the continuous problems are well posed.

The equations are discretized using a high order finite difference method on summation-by-parts form and the interface conditions are imposed weakly in a stable way, using penalty formulations. It is shown that one specific choice of penalty matrices leads to a dual consistent scheme and super-converging functionals.

By considering an example, it is shown that the correct physical coupling conditions are contained in the set of well posed coupling conditions. It is also shown that the correct convergence rates are obtained for both the solutions and functionals.

Key words. well posed problems, high order finite differences, stability, summation-by-parts, weak interface conditions, dual consistency

1. Introduction. Roughly speaking, a well posed initial boundary value prob-lem requires that a unique solution estimated in terms of data, exists. The most common procedure for showing well-posedness is the so called energy method. In this method, one multiplies the governing partial differential equations (PDEs) with the solution, integrate by parts and imposes a minimal number of boundary conditions, such that an energy estimate is obtained [21]. This procedure leads to a well posed problem. The same general knowledge is not wide-spread when it comes to the cou-pling of different multi-physics problems [15,19,18,27,8] at an interface. The reason for that is the more complex and to some extent more unclear nature of coupling conditions compared to boundary conditions [15,19,18,22,11].

Firstly, accuracy relations must exist such that a combination of variables for one set of PDEs at the interface is equal to another combination of variables for the other set. Secondly, the number of accuracy relations must fit both problems. Too many conditions ruin existence and too few ruin uniqueness. If the number of accuracy relations are too few, additional conditions requiring external data must be added. If the number of accuracy conditions are too many, only a subset can be used. Thirdly, the accuracy relations must be such that no artificial growth or decay is generated.

Sensitivity analysis for systems of PDEs are often used in optimization, param-eter estimation, model simplification, data assimilation, optimal control, uncertainty analysis and experimental design [13,12]. For problems involving a large number of sensitivity parameters, the adjoint method is the most efficient [14]. In this paper, we derive the adjoint equations and interface conditions for these types of problems, such that also the adjoint (or dual) problem is well posed and stable.

It has been shown that dual consistent summation-by-parts (SBP) discretizations lead to superconvergent functionals [9]. Additionally, in [2,3,4], it was shown that the dual boundary conditions depend on the primal boundary conditions and vice versa. The fact that the dual and primal boundary conditions are mutually dependent lead ∗Department of Mathematics, Computational Mathematics, Link¨oping University, SE-581 83

Link¨oping (fatemeh.ghasemi@liu.se).

Department of Mathematics, Computational Mathematics, Link¨oping University, SE-581 83

Link¨oping (jan.nordstrom@liu.se).

(4)

to new types of dual consistent boundary conditions. We will use the same technique and find dual consistent interface conditions.

The rest of the paper proceeds as follows. In Section2, the interface conditions, well-posedness and stability of the primal problem are derived. The dual problem, its interface conditions, well-posedness and stability are presented in Section3. We analyze dual consistency for the discrete problems in Section4. Section5contains an example of the developed theory and the numerical results for both the primal and the dual problems. Finally, in Section6 we summarize and draw conclusions.

2. The primal problem. We will consider interface conditions for the hyper-bolic systems (1) ut+ Aux= 0, −1 < x < 0, t > 0, u(x, 0) = f (x), −1 < x < 0, t = 0, and (2) vt+ Bvx= 0, 0 < x < 1, t > 0, v(x, 0) = g(x), 0 < x < 1, t = 0.

In (1)-(2), A, B are m × m and n × n symmetric constant matrices, respectively, and u, v are unknown vectors of sizes m and n. In general we have m 6= n. For clarity we often ignore the boundary conditions at x = ±1, and focus on the interface conditions at x = 0.

2.1. The interface conditions and well-posedness. The following definition of well-posedness is used.

Definition 1. Consider the coupled problem (1)-(2) with homogeneous boundary conditions and a minimal number of boundary and coupling conditions. The coupled problem is well posed if

d dt(kuk 2 2+ α p ckvk 2 2) ≤ 0,

where αpc is a positive free weight, kuk22=R0

−1u Tu dx and kvk2 2= R1 0 v Tv dx.

We apply the energy method by multiplying equations (1) and (2) with uT and vT

respectively, integrating in space and adding them together. The result is

(3) d dt(kuk 2 2+ α p ckvk 2 2) = w TEw| x=0.

In (3), w = [u, v]T and E = diag(−A, αp

cB). In the following, all terms are evaluated

at x = 0, unless stated otherwise.

We denote k = kA−+ k+B, as the number of positive eigenvalues of E, where kA−, kB+ are the number of positive eigenvalues of −A and B, respectively. This means that exactly k interface (or accuracy) conditions are required for a well posed problem. Let the interface conditions be described by

(4) Cpu = Dpv,

where the matrices Cpand Dphave k linearly independent rows. We will seek matrices

Cpand Dpsuch that the coupled problem (1)-(2) is well posed according to Definition

(5)

Since A and B are symmetric, we have (5) A = XΛAXT, X =[X+, X−], ΛA=  Λ+A 0 0 Λ−A  , B = Y ΛBYT, Y =[Y+, Y−], ΛB=  Λ+B 0 0 Λ−B  .

In (5), X+, X− and Y+, Y− are column matrices containing the eigenvectors related

to the positive and negative eigenvalues of A and B, respectively. The diagonal block matrices Λ+A, Λ−A and Λ+B, Λ−B contain the positive and negative eigenvalues of A and B, respectively. For simplicity (and without restriction) we assume that there are no zero eigenvalues.

By using (5), the relation (3) can be reformulated as

(6) d dt(kuk 2 2+ α p ckvk 2 2) = − (X Tu)TΛ A(XTu) + αpc(Y Tv)TΛ B(YTv) = −(XT +u) TΛ+ A(X T +u) + (X T −u) TΛ− A(X T −u)  + αpc(YT +v) TΛ+ B(Y T +v) + (Y T −v)TΛ−B(Y T −v).

We consider interface conditions of the general forms

(7) (Λ − AX T −− RAX+T)u =(TAY−T+ SAY+T)v, (Λ+BY+T − RBY−T)v =(TBX+T+ SBX−T)u.

The matrices RA, TA, SA, RB, TBand SBare of sizes kA−× k + A, k − A× k − B, k − A× k + B, k + B×

k−B, k+B× k+A and kB+× kA−, respectively, where k+A and k−B are the number of negative eigenvalues of −A and B, respectively.

2.1.1. Strongly imposed interface conditions. By inserting (7) into (6) we obtain (8) d dt(kuk 2 2+ α p ckvk 2 2) = −     XT +u XT −u YT +v YT −v     T    m11 m12 m13 m14 (m12)T m22 m23 m24 (m13)T (m23)T m33 m34 (m14)T (m24)T (m34)T m44     | {z } M     XT +u XT −u YT +v YT −v     , where m11=Λ+A+ R T A(Λ − A) −1R A− αpcT T B(Λ + B) −1T B, m23=0, m12=αpcT T B(ΛB)−1SB m24=αpcS T B(Λ + B)RB, m13=RTA(Λ−A) −1S A, m33=SAT(Λ−A) −1S A, m14=RTA(Λ − A) −1T A− αpcT T B(Λ + B) −1R B, m34=SAT(Λ − A) −1T A, m22= − αpcS T B(Λ + B) −1S B, m44= − αpc(Λ − B+ R T B(Λ + B) −1R B) +TAT(Λ−A) −1T A.

To obtain a well posed problem, the matrices RA, TA, SA, RB, TB, SBand αpc must be

(6)

Proposition 2. If SA 6= 0 or SB 6= 0, then the coupled problem (1)-(2) is not

well-posed.

Proof. If SA6= 0 or SB 6= 0, then the diagonal elements m22 or m33 are negative

and consequently the matrix M in (8) cannot be positive semi-definite.

From now on we consider SA= SB = 0, which leads to the interface conditions

(9) (Λ−AXT − RAX+T)u =TAY−Tv, (Λ + BY

T

+ − RBY−T)v = TBX+Tu,

where we specify the ingoing characteristic variables in terms of the outgoing ones and incoming data.

Remark 1. The form (9) automatically gives the correct number k = kA−+ k + B of

interface conditions.

For future reference we write the final coupled problem as

(10) ut+ Aux= 0, −1 < x < 0, t > 0, u(x, 0) = f (x), −1 < x < 0, t = 0, vt+ Bvx= 0, 0 < x < 1, t > 0, v(x, 0) = g(x), 0 < x < 1, t = 0, Cpu(0, t) = Dpv(0, t), x = 0, t > 0, where Cp=  TBX+T Λ−AXT −− RAX+T  , Dp=  Λ+BYT + − RBY−T TAY−T  .

Now, (8) can be rewritten as

(11) d dt(kuk 2 2+ α p ckvk 2 2) = − " XT +u YTv #T" mp 11 m p 12 (mp12)T mp22 # | {z } Mp " XT +u YTv # , where (12) mp11= m11, m p 12= m14, m p 22= m44.

In the following proposition, we will find matrices RA, RB, TA and TB such that the

coupled problem (10) leads to an energy estimate.

Proposition 3. If the matrices RA, RB, TA, TB and the parameter αpc are chosen

such that (13) Λ + A+ R T A(Λ−A) −1R A>0, TBT(Λ + B) −1T B≤(Λ+A+ R T A(Λ−A) −1R A)/αpc, Λ−B+ RBT(Λ+B)−1RB <0, TAT(Λ − A) −1T A≥αpc(Λ − B+ R T B(Λ + B) −1R B), and (14) RTA(Λ−A)−1TA= αpcT T B(Λ + B) −1R B,

then, the matrix Mp in (11) is positive semi-definite and an energy estimate is

(7)

Proof. Since (14) is satisfied, mp12= 0 and the matrix Mp has the form " mp 11 0 0 mp22 # ,

which is positive semi-definite due to the conditions in (13). Remark 2. Note that if mp11= 0 or m

p

22= 0 and m p

126= 0, the coupling can not

be done.

2.1.2. Weakly imposed interface conditions. In order to prepare for the numerical approximation, we impose the interface condition (4) weakly. The result is

(15) d dt(kuk 2 2+α p ckvk 2 2) = −u TAu + αp cv TBv +uTΣL(Cpu − Dpv) + (uTΣ p L(C pu − Dpv))T +αpcvTΣ R(Dpv − Cpu) + (vTΣR(Dpv − Cpu))T ,

where ΣL and ΣR are penalty matrices of sizes m × k and n × k, respectively. The

relation (15) can be rewritten in matrix form as

(16) d dt(kuk 2 2+ α p ckvk 2 2) = − " u v #T" np11 np12 (np12)T np 22 # | {z } Np " u v # , where (17) n p 11=A − ΣLCp− (ΣLCp)T, np12= ΣLDp+ αpc(ΣRCp)T, np22= − αpc B + ΣRDp+ (ΣRDp)T.

Next, we must find the matrices ΣLand ΣRsuch that the matrix Npin (16) is positive

semi-definite. Let ˜ Np=  I Ep 0 I T Np  I Ep 0 I  ,

and note that if ˜Np is positive semi-definite, so is Np. The choice Ep= −(np 11)−1n p 12, gives ˜ Np=  np11 0 0 −(np12)T(np 11)−1n p 12+ n p 22  . By choosing ΣL and ΣRsuch that

(18) np11> 0, −(np12)T(np11)−1np12+ np22≥ 0, then, the right-hand side of (16) is bounded.

The above procedure is formalized as

Proposition 4. By choosing the penalty matrices ΣL and ΣR such that (18) is

(8)

In the following proposition, we specify special choices of matrices ΣL and ΣR, which

lead to well-posedness. In Section4, it will be shown that the discrete approximation of (10) is dual consistent for these specific choices.

Proposition 5. If the matrices RA, TA, RB, TB and αpc satisfy (13) and (14),

then the matrices

(19) ΣL= X  0 0 0 IA−  , ΣR= −Y  IB+ 0 0 0  ,

guarantee that (18) is satisfied. In (19), IA− and IB+ are identity matrices of size k−A and k+B, respectively.

Proof. By inserting the matrices ΣL and ΣR into (17) we obtain

Np=     Λ+A RT A −α p cTBT 0 RA −Λ−A 0 TA −αp cTB 0 αpcΛ + B −αpcRB 0 TT A −αpcRTB −αpcΛ − B     .

The matrix Np can be split into

Np= ˜Mp+ ˜Np, M˜p=     mp11 0 0 mp12 0 0 0 0 0 0 0 0 (mp12)T 0 0 mp 22     ,

where mp11, mp12, mp22 are given in (12). The matrix ˜Mp is positive semi-definite and

corresponds to the result obtained for the strong interface conditions. Furthermore, we have ˜ Np=     −RT A(Λ − A)−1RA RTA 0 −RAT(Λ − A)−1TA RA −Λ−A 0 TA 0 0 0 0 −TT A(Λ − A)−1RA TAT 0 −T T A(Λ − A)−1TA     +αpc     +TBT(Λ+B)−1TB 0 −TBT +T T B(Λ + B)−1RB 0 0 0 0 −TB 0 Λ+B −RB +RT B(Λ + B) −1T B 0 −RTB +RTB(Λ + B) −1R B     ,

which is due to the weak imposition of interface conditions. The matrix ˜Np can be

rewritten as ˜ Np= LA(CA⊗ (Λ−A) −1)LT A+ α p cLB(CB⊗ (Λ+B) −1)LT B, where LA=     RA 0 0 0 0 Λ−A 0 0 0 0 I 0 0 0 0 TA     , CA=     −1 1 0 −1 1 −1 0 1 0 0 0 0 −1 1 0 −1     , LB =     TB 0 0 0 0 I 0 0 0 0 Λ+B 0 0 0 0 RB     , CB=     −1 0 1 −1 0 0 0 0 1 0 −1 1 −1 0 1 −1     ,

(9)

and ⊗ denotes the Kronecker product [10]. The eigenvalues of the matrices CA and

CB are {−3, 0, 0, 0}, which implies that the matrix ˜Np is positive semi-definite.

The difference between the right-hand side in the strong estimate (11) and the right-hand side in the weak estimate (16) is

R = −     XT +u XT −u YT +v YT −v     T ˜ Np     XT +u XT −u YT +v YT −v     .

We can expand the term R by using

CA=XAΛACX T A, XA= 1 √ 3     1 −1 0 1 −1 0 0 1 0 0 1 0 1 1 0 0     , ΛAC= diag([−3, 0, 0, 0]), CB =XBΛBCX T B, XB = 1 √ 3     1 0 −1 1 0 1 0 0 −1 0 0 1 1 0 1 0     , ΛBC = ΛAC,

and find that

R = −     RAX+Tu Λ−AXT −u Y+Tv TAY−Tv     T (XAΛACX T A⊗ (Λ − A) −1)     RAX+Tu Λ−AXT −u Y+Tv TAY−Tv     −αp c     TBX+Tu XT −u Λ+BYT +v RBY−Tv     T (XBΛBCX T B⊗ (Λ + B) −1)     TBX+Tu XT −u Λ+BYT +v RBY−Tv     = + (Λ−AXT− RAX+T)u − TAY−Tv T (Λ−A)−1 (Λ−1A XT− RAX+T)u − TAY−Tv  −αp c h (Λ+BY+T − RBY−T)v − TBX+Tu T (Λ+B)−1 (Λ+BY+T − RBY−T)v − TBX+Tu i ≤ 0.

Remark 3. The additional seemingly dissipative term R in the weak energy rate is proportional to the interface condition (9) squared and is obviously zero. A non-zero truly dissipative term of the same form will appear in the discrete approximation.

2.2. The semi-discrete primal problem. We now consider finite difference approximations of (10) on SBP form [23,6,16]. The interface conditions are imple-mented using simultaneous approximation terms (SAT) as described in [7,24,25,20]. The semi-discrete SBP-SAT formulation of (10) is,

(20) ut+ (Du⊗ A)u = P

−1

u eN ⊗ ΣΣΣL(CpuN− Dpv0),

(10)

where the outer boundary conditions are ignored as in the continuous case. The discrete solutions are arranged as

u = (u01, .., u0m, .., uN 1, .., uN m)T, v = (v01, .., v0n, .., vM 1, .., vM n)T,

where uN = [uN 1, . . . , uN m]T, v0 = [v01, . . . , v0n]T. In (20), Du,v = Pu,v−1Qu,v are

difference operators approximates the first derivative and the subscripts u, v denote the direction they operate along. The difference operators satisfy Pu,v = Pu,vT > 0

and Qu,v+ QTu,v = diag([−1, 0, . . . , 0, 1]). eN = (0, . . . , 1)T and e0 = (1, . . . , 0)T are

N + 1 and M + 1 unit vectors, respectively. The penalty matrices ΣΣΣL and ΣΣΣR have

the same dimensions as in the continuous case, and will be chosen for stability. For a comprehensive review on the SBP-SAT technique, see [26].

2.2.1. Stability conditions at the interface. The discrete energy method is applied to (20) by multiplying the two equations with uT(P

u⊗ I), vT(Pv ⊗ I)

respectively, and adding the result. By defining the discrete norms kuk2 Pu⊗I =

uT(Pu⊗ I)u, kvk2Pv⊗I = v

T(P

v⊗ I)v, using the symmetry properties of A, B and

the summation-by-parts property of Qu,v, we obtain,

(21) d dt(kuk 2 Pu⊗I+α p dkvk 2 Pv⊗I) = −u T NAuN + αpdvT0Bv0 +uTNΣΣΣL(CpuN − Dpv0) + (uTNΣΣΣL(CpuN − Dpv0))T +αpdvT 0ΣΣΣR(Dpv0− CpuN) + (v0TΣΣΣR(Dpv0− CpuN))T ,

which is similar to the continuous (weak) energy rate in (15). The relation (21) can be rewritten in the matrix form

(22) d dt(kuk 2 Pu⊗I+ α p dkvk 2 Pv⊗I) = −  uN v0 T Np  uN v0  , where we find that Np= Np in (16) by letting αp

d= α p c.

We immediately arrive at the following proposition.

Proposition 6. By choosing the penalty matrices ΣΣΣL and ΣΣΣR such that (18) is

satisfied, the semi-discrete approximation (20) of the coupled problem (10) is stable. Similar to the continuous case, we will specify a set of penalty matrices which leads to stability.

Proposition 7. If the matrices RA, TA, RB, TB and αpd satisfy (13) and (14),

then the matrices

(23) ΣΣΣL= X  0 0 0 IA−  , ΣΣΣR= −Y  IB+ 0 0 0  , satisfy (18).

Proof. See proof of Proposition5.

Remark 4. The derivation in the discrete case is analogous to the continuous one above due to the mimicking properties of the SBP-SAT technique. In fact, the interface conditions and penalty matrices are already derived in the continuous setting, see [17] for more details on this technique.

The discrete energy estimate is similar to the continuous one with the additional term

R =(Λ−AX−T− RAX+T)uN−TAY−Tv0 T (Λ−A)−1(Λ−1A X−T− RAX+T)uN−TAY−Tv0  −αpd   (Λ+BY T + − RBY−T)v0−TBX+TuN T (Λ+B) −1 (Λ+BY T + − RBY−T)v0−TBX+TuN  .

(11)

This term was identically zero in the continuous case, but now adds a small amount of dissipation, that vanishes with mesh refinement.

3. The dual problem. Consider the linear functional J (w) = (u, h)L+ (v, l)R,

where h and l are weight functions and (u, h)L=

R0 −1u Th dx, (v, l) R= R1 0 v Tl dx. To

derive the dual problem, we add on forcing functions fL, fR to the right-hand sides

of (1)-(2), and seek functions φ and ψ such that Z T 0 J (w)dt = Z T 0 (φ, fL)L+ (ψ, fR)R dt.

A formal computation gives Z T 0 J (w) dt = Z T 0 J (w) dt − Z T 0 (φ, ut+ Aux− fL)L dt − Z T 0 (ψ, vt+ Bvx− fR)R dt = Z T 0 (φ, fL)L+ (ψ, fR)R dt (24) + Z T 0 (u, φt+ Aφx+ h)Ldt + Z T 0 (v, ψt+ Bψx+ l)R dt + Z 0 −1 [φTu]T0 dx + Z 1 0 [ψTv]T0 dx − Z T 0 ([φTAu]0−1+ [ψTBv]10) dt. The dual boundary and interface conditions are the minimal number of conditions such that

(25)

Z T 0

([φTAu]0−1+ [ψTBv]10)dt = 0.

In the next subsection, explicit interface conditions will be derived. By choosing homogeneous initial and final conditions

u(x, 0) = v(x, 0) = φ(x, T ) = ψ(x, T ) = 0, for the primal and dual problem, the terms R0

−1[φ Tu]T 0 dx and R1 0[ψ Tv]T 0 dx in (24)

vanish, and the following dual equations are obtained (26) −φt− Aφx= h, −1 < x < 0,

−ψt− Bψx= l, 0 < x < 1.

3.1. The dual boundary and interface conditions. The interface conditions (9) grouped together with the interface terms in (25) yields

0 = (φTAu − ψTBv) =(φTX+Λ+A+ φ TX −RA− ψTY+TB)X+Tu +(φTX−TA− ψTY+RB− ψTY−Λ−B)Y T −v,

and arrive at the dual interface conditions

(12)

Remark 5. Note that the dual interface conditions are given by the primal ones. See [2,3,4] for similar effects regarding boundary conditions.

The dual interface condition (27) can be written more compact as

(28) Caφ = Daψ, where Ca=  Λ+AXT ++ RTAX T − TT AX T −  , Da=  TT BY+T Λ−BYT − + RTBY T +  .

For completeness, the dual boundary conditions must also be determined such that

(29) φTAu|x=−1= 0, ψTBv|x=+1= 0.

To choose dual boundary conditions, the boundary conditions for the primal problem must exist. We consider the following general homogeneous boundary conditions for the primal problem

(Λ+AX+T − RlX−T)u(−1, t) = 0, (Λ−BY T

− − RrY+T)v(+1, t) = 0,

where the matrices Rland Rr are such that a well posed primal problem is obtained.

This yields φTAu|x=−1= (φTX+Rl+ φTX−Λ−A)X T −u|x=−1, ψTBv|x=+1= (ψTY+Λ+B+ ψ TY −Rr)Y+Tv|x=+1,

and we choose the dual boundary conditions

(Λ−AXT+ RTlX+T)φ(−1, t) = 0, (Λ+BY+T + RTrYT)ψ(+1, t) = 0, such that all boundary terms vanish and (29) is satisfied.

The initial conditions for the dual problem are given at time t = T. The time transformation τ = t − T inserted into (26) results in

(30) φt− Aφx= h, −1 < x < 0, τ > 0, φ(x, 0) = 0, −1 < x < 0, τ = 0, ψt− Bψx= l, 0 < x < 1, τ > 0, ψ(x, 0) = 0, 0 < x < 1, τ = 0, Caφ(0, τ ) = Daψ(0, τ ), x = 0, τ > 0,

which is the final form of the coupled dual problem. Note that the dual interface conditions Caφ = Daψ, given in (28) specify the ingoing characteristic variables in terms of outgoing ones and incoming data, just as in the primal problem.

3.1.1. Strongly imposed dual interface conditions. We apply the energy method to (30), and ignore the boundary terms as for the primal problem. By using

(13)

the dual interface conditions (27), we get (31) d dτ(kφk 2 2+ α a ckψk 2 2) =φ TAφ| x=0− αacψ TBψ| x=0 = " XT −φ YT +ψ #T" ma 11 ma12 (ma 12)T ma22 # | {z } Ma " XT −φ YT +ψ # , where ma11=Λ−A+ RA(Λ+A) −1RT A− α a cTA(Λ−B) −1TT A, ma12= − RA(Λ+A) −1TT B + α a cTA(Λ−B) −1RT B, ma22= − αacΛ+B− αa cRB(Λ−B) −1RT B+ TB(Λ+A) −1TT B, and αa c is a positive weight.

In the following proposition, we will find the matrices RA, RB, TA and TB such

that the coupled problem (30) leads to an energy estimate.

Proposition 8. If the matrices RA, RB, TA, TB and the parameter αac are chosen

such that (32) Λ − A+ RA(Λ+A) −1RT A<0, TA(Λ−B) −1TT A ≥(Λ − A+ RA(Λ+A) −1RT A)/α a c, Λ+B+ RB(Λ−B) −1RT B >0, TB(Λ+A) −1TT A ≤α a c(Λ + B+ RB(Λ−B) −1RT B), and (33) RA(Λ+A) −1TT B = α a cTA(Λ−B) −1RT B,

then, the matrix Ma in (31) is negative semi-definite and an energy estimate is sat-isfied for (30).

Proof. The proof proceeds as in Proposition3.

3.1.2. Weakly imposed dual interface conditions. As in the primal prob-lem, we impose the interface condition (28) weakly. The result is

(34) d dτ(kφk 2 2+α a ckψk 2 2) = φ TAφ − αa cψ T +φTΣL(Caφ − Daψ) + (φTΣL(Caφ − Daψ))T +αacψTΣ R(Daψ − Caφ) + (ψTΣR(Daψ − Caφ))T ,

where the penalty matrices ΣLand ΣRare of size m × (k+A+ k − B) and n × (k + A+ k − B), respectively.

Remark 6. The derivation procedure below for the dual problem is the same as for the primal problem.

The relation (34) can be rewritten in the matrix form

(35) d dτ(kφk 2 2+ α a ckψk 2 2) = " φ ψ #T" na 11 na12 (na12)T na22 # | {z } Na " φ ψ # ,

(14)

where

na11=A + ΣLCa+ (ΣLCa)T, na12= −ΣLDa− αac(ΣRCa)T,

na22= − α a

c B − ΣRDa− (ΣRDa)T.

By using the rotation technique, as in the primal problem, we arrive at Proposition 9. The choice of penalty matrices ΣL and ΣR such that

(36) na11< 0, −(na 12) T(na 11) −1na 12+ n a 22≤ 0,

holds, leads to well-posedness of coupled problem (30).

As in the primal problem, we specify special choices of matrices ΣL and ΣR, which

lead to well-posedness.

Proposition 10. If the matrices RA, TA, RB, TB and αac satisfy (32) and (33),

then the matrices

(37) ΣL= −X  IA+ 0 0 0  , ΣR= Y  0 0 0 IB−  ,

guarantee that (36) holds. In (37), IA+ and IB− are identity matrices of size kA+ and k−B, respectively.

Proof. See the proof of Proposition4.

3.2. The semi-discrete approximation of the dual problem. The corre-sponding semi-discrete SBP-SAT formulation of (30) is,

(38) φφφτ− (Du⊗ A)φφφ = P −1 u eN⊗ ΣΣΣL(CaφN− Daψ0), ψ ψ ψτ− (Dv⊗ B)ψψψ = Pv−1e0⊗ ΣΣΣR(Daψ0− CaφN),

where the outer boundary conditions are ignored as in the primal case. The vectors φN and ψ0are arranged as φN = [φN 1, . . . , φN m]T, ψ0= [ψ01, . . . , ψ0n]T. The penalty

matrices ΣΣΣL and ΣΣΣR have the same dimensions as in the continuous case.

3.2.1. Stability conditions at the interface. Applying the discrete energy method to (38) leads to (39) d dτ(kφφφk 2 Pu⊗I+α a dkψψψk 2 Pv⊗I) = φ T NAφN − αadψ T 0Bψ0 +φTNΣΣΣL(CaφN − Daψ0) + (φTNΣΣΣL(CaφN − Daψ0))T +αadψ0TΣΣΣR(Daψ0− CaφN) + (ψT0ΣΣΣR(Daψ0− CaφN))T ,

which is similar to the continuous case (34). The relation (39) can be rewritten in the matrix form d dt(kuk 2 Pu⊗I+ α p dkvk 2 Pv⊗I) =  φN ψ0 T Na  φN ψ0  , where we find that Na= Na in (35) by letting αp

a= αac.

We immediately arrive at the following proposition.

Proposition 11. By choosing the penalty matrices ΣΣΣL and ΣΣΣRsuch that (36) is

satisfied, the semi-discrete approximation (38) of the coupled problem (30) is stable. Similar to the continuous case, we will specify a set of penalty matrices which leads to stability.

(15)

Proposition 12. If the matrices RA, TA, RB, TB and αad satisfy (32) and (33),

then the matrices

(40) ΣΣΣL= −X  IA+ 0 0 0  , ΣΣΣR= Y  0 0 0 IB−  , guarantee that (36) holds and that (38) is stable.

Proof. See proof of Proposition5.

Remark 7. Just as in the primal problem, the penalty matrices in the discrete case are the same as in the continuous one, due to the similarity in analysis of the continuous and discrete problem.

4. Dual consistency. To investigate dual consistency, we rewrite (20) and (38) into the following form

(41) P wt+Lpw = P F, P θθθτ+ Laθθθ = P H, where Lp= +  Qu⊗ A 0 0 Qv⊗ B  +           0 . .. −ΣΣΣLCp ΣΣΣLDp Σ Σ ΣRCp −ΣΣΣRDp . .. 0           , La = −  Qu⊗ A 0 0 Qv⊗ B  +           0 . .. −ΣΣΣLCa ΣΣΣLDa Σ Σ ΣRCa −ΣΣΣRDa . .. 0           , and P−1 =  Pu−1⊗ Im 0 0 P−1 v ⊗ In  . In (41), w = [u, v]T, θθθ = [φφφ, ψψψ]T and F = [F L, FR], H = [HL, HR], where FL, FR, HL

and HR are the discrete form of fL, fR, hL and hR, respectively.

Remark 8. The similarity between the second matrix in Lp and La is due to the similarity of right-hand side in (20) and (38).

The schemes in (41) are dual consistent [9,2] if

(42) (Lp)T = La.

Since Qu,v+ QTu,v= EN − E0, the requirement (42) leads to the condition

(43)  A − (Cp)TΣΣΣT L (Cp)TΣΣΣTR (Dp)TΣΣΣT L −(Dp)TΣΣΣTR− B  =−ΣΣΣLC a ΣΣΣ LDa Σ Σ ΣRCa −ΣΣΣRDa  . We find

(16)

Proposition 13. If the penalty matrices ΣΣΣL,R and ΣΣΣL,R are chosen as in (23)

and (40), respectively, then the SBP-SAT discretization of the coupled problem (10) is dual consistent

Proof. A direct insertion of (23) and (40) into (43) yields the result.

Remark 9. By choosing the penalty matrices such that (18) is satisfied but (43) is not, the semi-discrete approximation (20) is stable but dual inconsistent.

Remark 10. We summarize what has been done so far below:

• Well posed interface conditions for the primal problem have been derived. The interface conditions are imposed and penalty matrices are obtained such that the continuous problem is well posed. The penalty matrices for the continuous problem lead directly to a stable discrete primal problem.

• By using the same strategy as in the primal problem, the interface conditions and penalty matrices for the continuous coupled dual problem are derived. These matrices lead to a stable discrete dual problem. The dual interface condition is determined by the primal interface condition.

• A specific set of penalty matrices for the primal and dual problem can be cho-sen such that the discrete problems are stable and the primal discrete problem is dual consistent.

5. An example. As an example, we consider the linearized symmetrized Euler equations

(44) ut+ Aux= 0, −1 < x < 0,

coupled to the elastic equation

(45) vtt− c2vxx= 0, 0 < x < 1. In (44), A =     ¯ w c/¯ √γ 0 ¯ c/√γ w¯ c¯qγ−1γ 0 c¯qγ−1γ w¯     , (46) where u = " ¯ cρ √ γ ¯ρ, w, −¯cρ ¯ ρpγ(γ − 1) + ppγ/(γ − 1) ¯ ρ¯c #T ,

and ρ, w and p are respectively the density, the velocity and the pressure perturbation. The solution we have linearized around is denoted by overbars and γ is the ratio of specific heats [1]. In (45), v is the displacement and c = pE/ρs, where E is the

elastic modulus and ρsis the density of solid.

The second order equation (45) can be rewritten as a first order system given by

Ut+ BUx= 0, B =  0 c c 0  , (47)

where U = [vt, −cvx]T and vtand vxare the velocity and the stress, respectively.

We will investigate if the mathematical theory derived earlier will provide us with the physical coupling conditions. The physical interface conditions come from mechan-ical principles: (1) Continuity and (2) Force balance (Newton’s II law). According

(17)

to the first principle the fluid and solid velocities must match at the interface, other-wise the fluid will be detached from the solid. This means that one of the interface conditions is

(48) w = vt.

The fluid pushes on the solid with a traction p and the solid pushes back with an equal and opposite traction, which is called σ, and we get σ = −p. For a linear elastic material, according to Hooke’s law, we have σ = Evx, which implies that the second

interface condition is

(49) p = −Evx.

Now, we apply the mathematical approach. The matrices A and B in (46) and (47) can be written as A = X+Λ+AX T ++ X−ΛAX−T, B = Y+Λ+BY T + + Y−ΛBY−T, where X+= 1 √ 2   q2(γ−1) γ 0 −p2/γ p1/γ 1 q (γ−1) γ   T , X− = 1 √ 2[ p 1/γ, −1,p(γ − 1)/γ]T, Y+= 1 √ 2[1, 1] T, Y − = 1 √ 2[1, −1] T, and Λ+A=  ¯ w 0 0 w + ¯¯ c  , Λ−A= ¯w − ¯c, Λ+B = c, Λ−B = −c.

According to Proposition3, the unknown matrices RA, RB, TAand TBmust be chosen

such that the conditions (13) and (14) are satisfied. There are different choices for these unknown matrices which lead to well-posedness.

One of the choices is the characteristic interface conditions (50) RA= [0, 0] , TA= p c(¯c − ¯w), RB = 0, TB = h 0,pc( ¯w + ¯c)i, and αp

c = 1. But the choices in (50) do not lead to the physical conditions (48) and

(49). Hence, these conditions lack accuracy. By inspection we find that the matrices

(51) RA=  0,(¯c − ¯w)( ¯ρ¯c − E/c) ¯ ρ¯c + E/c  , RB = − c( ¯ρ¯c − E/c) ¯ ρ¯c + E/c , TA= 2E/c(¯c − ¯w) ¯ ρ¯c + E/c , TB =  0, 2c ¯ρ¯c ¯ ρ¯c + E/c  , leads to w + p = vt− Evx, ρ¯¯cw + p = ¯ρ¯cvt− Evx,

which is equivalent to the the physical conditions (48) and (49) and hence accurate. The choice (51), and αp

c = ρs(¯c − ¯w)/ ¯ρ¯c implies that the conditions (13) and (14) are

(18)

Finally, we apply also boundary conditions to the model problem. Since two eigenvalues of matrix A are positive, two boundary conditions are needed at x = −1. The matrix B has one negative eigenvalue, which means that one boundary condition is needed at x = 1. We choose the characteristic boundary conditions

X+Tu(−1, t) = gl(t), Y−Tv(1, t) = gr(t).

5.1. Numerical results. In the numerical calculations, we will use the manu-factured solutions

u1(x, t) =u2(x, t) = u3(x, t) = et/λsin(2π(x)),

v1(x, t) =v2(x, t) = et/λsin(3π(x − t)),

with λ = 0.05 and the parameters ¯w = 0.5, ¯c = 1, ¯ρ = 1, γ = 1.4 and c = 1. The manufactured solution will provide data for the forcing function, boundary conditions and initial function. The rate of convergence is calculated as

qu= ln kuN1− uk Pu⊗I kuN2− ukP u⊗I / ln N1 N2 , qv = ln kvN1− vk Pv⊗I kvN2− vkP v⊗I / ln N1 N2 , where u and v are the analytical solutions and uNi and vNi are the corresponding

numerical solutions with Ni grid points.

The time-integration in this section is done using the classical fourth order explicit Runge-Kutta scheme [5] with CF L number = 0.01. We choose the final time T = 1.0 with time step ∆t = CF L × ∆x, which makes the time-error negligible.

In Tables 1 and 2, we show the convergence rates for different SBP operators, which confirms that the scheme is design order accurate.

5.2. Superconvergence of the functional. In this section, we investigate the impact of dual consistency. A dual consistent scheme is obtained by choosing the penalty matrices as in (23). Small changes in the coefficients of the penalty matrices in (23), leads to a stable but dual inconsistent scheme. In the dual inconsistent case, we choose Σ Σ ΣL= 3 2X  0 0 0 IA−  , ΣΣΣR= − 3 2Y  IB+ 0 0 0  , which leads to stability.

In Tables3and4, the rates of convergence of the functional

J (w) = Z 0 −1 u1+ u2+ u3 dx + Z 1 0 v1+ v2 dx,

are calculated and clear superconvergence can be seen. The superconvergent func-tional is obtained without requiring any knowledge about the solution of the dual equations. In order to have a dual consistent scheme, we only need the dual equation and its interface conditions, which are obtained by knowing the interface conditions of the primal problem. This means that superconverging functionals are obtained at no extra computational cost [2].

6. Summary and conclusions. We have considered the coupling of two general hyperbolic systems. The energy method was used to derive general well posed interface

(19)

N SBP 21 SBP 42 SBP 63 SBP 84 20 - - - -40 2.009 3.826 5.135 4.189 80 2.002 4.084 5.307 4.630 120 2.000 4.165 5.379 4.690 180 2.000 4.190 5.412 4.712 Table 1

Rate of convergence qufor u = (u1, u2, u3).

N SBP 21 SBP 42 SBP 63 SBP 84 20 - - - -40 1.983 3.219 4.618 3.851 80 2.003 3.210 4.717 4.718 120 2.002 3.140 4.761 4.846 180 2.001 3.095 4.794 4.826 Table 2

Rate of convergence qvfor v = (v1, v2).

N SBP 21 SBP 42 SBP 63 SBP 84 20 - - - -40 1.999 3.247 4.384 3.355 80 1.986 3.249 4.861 4.061 120 1.993 3.172 4.945 4.824 180 1.997 3.117 4.834 5.114 Table 3

Rate of functional convergence for the dual inconsistent discretization.

N SBP 21 SBP 42 SBP 63 SBP 84 20 - - - -40 1.958 3.843 5.419 6.400 80 1.961 4.072 5.984 7.831 120 1.984 4.122 5.900 8.373 180 1.992 4.116 6.200 8.521 Table 4

Rate of functional convergence for the dual consistent discretization.

conditions. It was shown that the derived interface conditions lead to a well posed problem both for weak and strong imposition.

The equations were discretized using finite differences on SBP-SAT form. The penalty matrices were derived in the analysis of the continuous problem. Almost no additional derivations were necessary.

Next, the dual problem and its well posed interface conditions were derived using the energy method. The interface conditions were imposed weakly and strongly also

(20)

for the dual problem. The weak interface procedures lead directly to stability of the numerical approximation as in the primal problem. The numerical scheme was shown to be dual consistent for specific choices of the penalty matrices.

By considering a physical example, it was shown that the mathematical inter-face conditions contain the physically correct interinter-face conditions. The mathematical theory can also narrow the search for well posed and accurate interface conditions.

The rate of convergence was verified by the method of manufactured solution and the result was consistent with the SBP-SAT theory. Superconvergent rates were obtained for the dual consistent discretization, at no extra cost.

REFERENCES

[1] S. Abarbanel and D. Gottleib, Optimal time splitting for two- and three-dimensional Navier-Stokes equations with mixed derivatives, Journal of Computational Physics, 41 (1981), pp. 1–43.

[2] J. Berg and J. Nordstr¨om, Superconvergent functional output for time-dependent problems using finite differences on summation-by-parts form, Journal of Computational Physics, 231 (2012), pp. 6846–6860.

[3] J. Berg and J. Nordstr¨om, On the impact of boundary conditions on dual consistent finite difference discretizations, Journal of Computational Physics, 236 (2013), pp. 41–55. [4] J. Berg and J. Nordstr¨om, Duality based boundary conditions and dual consistent finite

dif-ference discretizations of the Navier-Stokes and euler equations, Journal of Computational Physics, 259 (2014), pp. 135–153.

[5] J. C. Butcher, Coefficients for the study of runge-kutta integration processes, Journal of the Australian Mathematical Society, 3 (1963), pp. 185–201.

[6] M. Carpenter, J. Nordstr¨om, and D. Gottlieb, A stable and conservative interface treat-ment of arbitrary spatial accuracy, Journal of Computational Physics, 148 (1999), pp. 341– 365.

[7] M. H. Carpenter, D. Gottlieb, and S. Abarbanel, Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: methodology and application to high-order compact schemes, Journal of Computational Physics, 111 (1994), pp. 220–236. [8] B. Flemisch, M. Kaltenbacher, and B. I. Wohlmuth, Elasto-acoustic and acoustic-acoustic

coupling on nonmatching grids, International Journal for Numerical Methods in Engineer-ing, 67 (2006), pp. 1791–1810.

[9] J. E. Hicken and D. W. Zingg, Superconvergent functional estimates from summation-by-parts finite-difference discretizations, SIAM Journal on Scientific Computing, 33 (2011), pp. 893–922.

[10] R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge University Press, 1991.

[11] X. Huan, J. E. Hicken, and D. W. Zingg, Interface and boundary schemes for high-order methods., The 19th AIAA fluid dynamics conference, AIAA paper no. 2009-3658, San Antonio, USA, (2009), pp. 22–25.

[12] A. Jameson, Aerodynamic design via control theory, Journal of Scientific Computing, 3 (1988), pp. 233–260.

[13] Y. Jarny, M. N. Ozisik, and J. P. Bardon, A general optimization method using adjoint equation for solving multidimensional inverse heat conduction, International Journal of Heat and Mass Transfer, 34 (1991), pp. 2911–2919.

[14] S. Li and L. Petzold, Adjoint sensitivity analysis for time-dependent partial differential equations with adaptive mesh refinement, Journal of Computational Physics, 198 (2004), pp. 310–325.

[15] J. Lindstr¨om and J. Nordstr¨om, A stable and high-order accurate conjugate heat transfer problem, Journal of Computational Physics, 229 (2010), pp. 5440–5456.

[16] K. Mattsson, Boundary procedures for summation-by-parts operators, Journal of Scientific Computing, 18 (2003), pp. 133–153.

[17] J. Nordstr¨om, Well posed problems and boundary conditions in computational fluid dynam-ics, 22nd AIAA Computational Fluid Dynamics Conference no. 2015-3197, Texas, USA, (2015), pp. 1–14.

[18] J. Nordstr¨om and J. Berg, Conjugate heat transfer for the unsteady compressible Navier-Stokes equations using a multi-block coupling, Computers & Fluids, 72 (2013), pp. 20–29.

(21)

[19] J. Nordstr¨om and S. Eriksson, Fluid structure interaction problems: the necessity of a well posed, stable and accurate formulation, Communications in Computational Physics, 8 (2010), pp. 1111–1138.

[20] J. Nordstr¨om, J. Gong, E. van der Weide, and M. Sv¨ard, A stable and conservative high order multi-block method for the compressible Navier-Stokes equations, Journal of Computational Physics, 228 (2009), pp. 9020–9035.

[21] J. Nordstr¨om and M. Sv¨ard, Well-posed boundary conditions for the Navier-Stokes equa-tions, SIAM Journal on Numerical Analysis, 43 (2005), pp. 1231–1255.

[22] B. Roe, R. Jaiman, A. Haselbacher, and P. H. Geubelle, Combined interface bound-ary condition method for coupled thermal simulations, International Journal for numerical methods in fluids, 57 (2008), pp. 329–354.

[23] B. Strand, Summation by parts for finite difference approximations for d/dx, Journal of Com-putational Physics, 110 (1994), pp. 47–67.

[24] M. Sv¨ard, M. Carpenter, and J. Nordstr¨om, A stable high-order finite difference scheme for the compressible Navier-Stokes equations: far-field boundary conditions, Journal of Computational Physics, 225 (2007), pp. 1020–1038.

[25] M. Sv¨ard and J. Nordstr¨om, A stable high-order finite difference scheme for the compress-ible Navier-Stokes equations: no-slip wall boundary conditions, Journal of Computational Physics, 227 (2008), pp. 4805–4824.

[26] M. Sv¨ard and J. Nordstr¨om, Review of summation-by-parts schemes for initial-boundary-value problems, Journal of Computational Physics, 268 (2014), pp. 17–38.

[27] K. Virta and K. Mattsson, Acoustic wave propagation in complicated geometries and het-erogeneous media, Journal of Scientific Computing, 61 (2014), pp. 90–118.

References

Related documents

Inhyrningen av teknikkonsulter skulle därför kunna innebära fördelar för företaget eftersom de fastanställda kan fördjupa sin kunskap och kompetens om ett

Nedan följer en kort beskrivning av tillvägagångssättet för en fallstudie utifrån Creswell (2007) och Stake (1995). Till att börja med måste forskaren fastställa om

Figure 3.19 shows the electric field magnitude and charge density along the electrode axis.. Figure 3.20 shows the electric field magnitude as a function of the

De kvinnor som hade en BRCA- mutation berättade för färre i sin omgivning om sin genetiska mutation till skillnad från de som inte bar på en BRCA-mutation.. Kvinnornas fäder och barn

In Section 8 we derived formulas to compute the overhead due to the JAMMY join procedure. Let us now consider the overhead of the join process when the centralized solution is used.

Ett negativt värde per hektar innebär att anläggningen är före- tagsekonomiskt olönsam, vilket kan tolkas som att ”kostnaden” för att tillgodogöra sig de fördelar som

Torkning och lagring, alternativ Leverans till central tork Egen tork och lagring Värde vid skördeleverans Värde vid leverans i Pool 2 Arbets- och maskinkostnad Arbets-

Genom att som sjuksköterska lyssna på patienten, förmedla kunskap om diabetessjukdomen, vara öppen och lyhörd för den enskilde patientens behov samt uppmuntra patienten till