• Nu S-Au Găsit Rezultate

View of Extending Broyden's method to interaction problems

N/A
N/A
Protected

Academic year: 2022

Share "View of Extending Broyden's method to interaction problems"

Copied!
12
0
0

Text complet

(1)

Rev. Anal. Num´er. Th´eor. Approx., vol. 37 (2008) no. 2, pp. 169–180 ictp.acad.ro/jnaat

EXTENDING BROYDEN’S METHOD TO INTERACTION PROBLEMS

ROBBY HAELTERMAN, JORIS DEGROOTE, JAN VIERENDEELS and DIRK VAN HEULE

Abstract. The solution of problems involving the interaction of different sys- tems is a domain of ongoing research, although often a good solver already exists for each system separately. In this paper we draw our ideas from one of the best known all-round quasi-Newton methods: Broyden’s rank-one update, which we extend to algorithms using 2 approximate Jacobians. A comparison is made with the iterative substructuring method and Aitken’s acceleration method. It is shown that a Broyden method using only a single approximate Jacobian per- forms best.

MSC 2000. 90C53, 15A06, 65F10, 74F10.

Keywords. Quasi-Newton method, iterative method.

1. INTRODUCTION

Simulations of coupled systems are becoming ever more important. While powerful solvers often exist for problems in a single physical domain (e.g.

structural problems, flow problems, ...), development of similar tools for multi- physics problems is still ongoing, and these can be broadly put in one of the following categories:

Monolithic or simultaneous solution: the whole problem is treated as a monolithic entity and solved with a specialized solver.

Partitioned solution: the physical components are treated as isolated entities that are solved separately.

The relative merits of these methods are very problem dependent: the ad- vantage of the monolithic approach is the enhanced stability. This comes at a cost, however, as specialized software has to be written for each type of interaction problem, which can result in very large systems. The partitioned approach allows for the use of available specialized solvers for each physical component (structure, fluid,...), on condition that the coupling effects can be treated efficiently. The latter is often no problem for loosely coupled problems,

Department of Mathematics, Royal Military Academy, Renaissancelaan 30, B-1000 Brus- sels, Belgium, e-mail: [email protected].

Department of Flow, Heat and Combustion, Ghent University, St.-Pietersnieuwstraat 41, B-9000 Ghent, Belgium, e-mail: [email protected].

(2)

which can be solved by sequentially calling the solvers for each component.

Strongly coupled problems, on the other hand, still pose a real challenge; the simplest solution method is theiterative substructuring method, but this often displays very poor stability [4].

In this paper we write the coupled problem as a root-finding problem which we solve with Broyden’s quasi-Newton method and an extension of this method in which two approximate Jacobians are used, one for each physical problem.

This results in a partitioned approach in which each of the physical problems is solved with a specialized code that we consider to be a black box and of which the Jacobian is unknown. We assume the dominant cost lies in the use of these black boxes.

This paper is organized as follows: in section 2 we describe the type of problems we are interested in; in section 3, resp. 4, we describe two well- known methods for the solution of these problems: the iterative substructuring method, resp. Aitken’s ∆2 method; in section 5 we briefly describe Broyden’s quasi-Newton method together with two proposed variations; in section 6 we test the various algorithms on two coupled problems.

2. THE GENERAL MODEL

We consider coupled problems that can be stated in the following form:

F(g) = p, (1a)

S(p) = g, (1b)

where F : DF ⊂ RM×1 → RN×1 : g 7→ F(g) and S : DS ⊂ RN×1 → RM×1 :p7→S(p). Each equation describes a physical problem that is spatially decomposed. We limit p and g to values on the interface between the two physical problems. In this way the physically decoupled nature of the problem is exploited, and the method can be regarded as a special case ofheterogeneous domain decomposition methods [2].

We want to solve (1) with a partitioned method that keeps the available solvers for each individual equation (1a), resp. (1b), of which we lack knowl- edge of the Jacobian. (1) can also be written as the root-finding problem

F(S(p))−p=H(p)p=K(p) = 0.

(2)

Depending on the algorithm we will either use the formulation of (1) or (2).

We assume thatF,S,H and K satisfy the following properties:

(1) F, resp S, H, is twice continuously differentiable in an open set DF, respDS, DH.

(2) There is a po such that the set C(po) = {p;kK(p)k ≤ kK(po)k} is compact.

(3) K(p) = 0 has a single solutionp inC(po)

(4) (K0(p))−1 exists and is continuous in an open set containing C(po).

(3)

Remark. The choice, in equation (2), betweenF(S(p)) = 0 orS(F(g)) = 0, will depend on the relative dimensions of pandg; we will choose the formula- tion such that the (square) Jacobian of K is of the lowest dimension.

3. ITERATIVE SUBSTRUCTURING METHOD

The simplest strongly coupled method to solve (1) is the iterative substruc- turing method; it proceeds as follows:

Algorithm 3.1 (Iterative substructuring method (ISM)). 1. Startup. Take an initial valuepo.

Sets= 0.

2. Loop until sufficiently converged a. gs=S(ps)

b. ps+1 =F(gs) c. s=s+ 1

Depending on the strength of the coupling between the two equations in (1) this method is not always stable [4].

4. AITKEN’S2 METHOD

Another method that is often used for the solution of (2) is Aitken’s ∆2 method [3].

Algorithm 4.1 (Aitken’s ∆2 method (A∆2)). 1. Startup:

Take a starting value po and θo; computeK(po) =H(po)−po; computep1=po+θoK(po);

set s= 1.

2. Loop until sufficiently converged:

a. Compute K(ps) =H(ps)−ps. b. θs=−θs−1K(phδKs−1)TδKs−1

s−1,δKs−1i

c. ps+1 =ps+θsK(ps) d. s=s+ 1

(δKs−1 =K(ps)−K(ps−1) and h·,·iis the classical scalar product.)

5. BROYDEN’S METHOD AND VARIATIONS

A well-known solution method for the nonlinear problem in (2) is Newton’s method. However, as we suppose we don’t have access to the Jacobians of eitherF,S,HorK, we will resort to quasi-Newton methods, which we apply either to the single equation (2) or the system (1).

(4)

5.1. Basic Broyden method. Broyden’s method [1] was developed for a sin- gle equation as in (2). The resulting algorithm can be written as follows.

Algorithm 5.1 (Broyden’s algorithm (Br)). 1. Startup:

Take a starting value po and ˆKo0; computeK(po) =H(po)−po;

computep1= (1−ω)po+ωH(po) (ω being a relaxation parameter);

set s= 1.

2. Loop until sufficiently converged:

a. Compute K(ps) =H(ps)−ps.

b. Compute the Broyden update ˆKs0 = ˆKs−10 +(δKs−1Kˆ

0

s−1δps−1)δpTs−1 hδps−1,δps−1i

c. ps+1 =ps−( ˆKs0)−1K(ps) d. s=s+ 1

(δps−1 =psps−1.)

A good initial guess for ˆKo0 is needed. As we do not have any idea what the Jacobian of F or S is, we will set them equal to zero, such that we have Kˆo0 =−I as a consequence.

5.2. Extension of Broyden’s method. We recall that K(p) =F(S(p))−p and that as such

K0(ps) = F0(S(ps))S0(ps)−I.

We therefore propose to build approximate Jacobians ˆFs0 and ˆSs0 along the lines of the Broyden update:

Fˆs0 = ˆFs−10 +(δFs−1Fˆ

0

s−1δgs−1)δgs−1T hδgs−1,δgs−1i , (3a)

Sˆs0 = ˆSs−10 +(δSs−1Sˆ

0

s−1δps−1)δpTs−1 hδps−1,δps−1i , (3b)

where δgs−1 = gsgs−1, δFs−1 = F(gs)−F(gs−1), and δSs−1 = S(ps)− S(ps−1); we set ˆKs0 = ˆFs0Sˆs0I.

To keep in line with ˆKo0 in§5.1 we set ˆFo0 = ˆS0o= 0.

It is easy to see that this construction respects the following secant equa- tions:1

Fˆs0δgs−1 = δFs−1, (4a)

Sˆs0δps−1 = δSs−1, (4b)

Kˆs0δps−1 = δKs−1. (4c)

1Broyden’s original update also respects the secant equation.

(5)

This now allows us to write the following algorithm:

Algorithm5.2 (Broyden’s method–2 Jacobians–Variant 1 (BrVar1)). 1. Startup:

Take a starting value po, ˆFo0 and ˆSo0;

computego =S(po) and p1 = (1−ω)po+ωF(go)(ω being a relaxa- tion parameter);

set s= 1.

2. Loop until sufficiently converged a. Compute gs=S(ps).

b. Make approximate Jacobian ˆSs0 according to (3b).

c. ComputeH(ps) =F(gs).

d. Make approximate Jacobian ˆFs0 according to (3a).

e. Quasi-Newton step: ps+1 =ps−( ˆFs0Sˆs0I)−1(H(ps)−ps) f. s=s+ 1

Remark. While Broyden’s original update also respects (4c), it can be shown that the approximate Jacobians differ between algorithms 5.1 and 5.2.

In [5] another approach was proposed to couple (1) with two approximate Jacobians, which we now apply to Broyden’s method. In this approach we look at each equation in (1) separately and write a Taylor expansion in which we insert the approximate Jacobians:

F(gs+1)−F(gs) = Fˆs0(gs+1gs), S(ps+1)−S(ps) = Sˆs0(ps+1ps).

We require the solution of the coupled problem to satisfy (1) at iteration s+ 1; it follows that

pFˆs0g = F(gs)−Fˆs0gs, (5a)

gSˆs0p = S(ps)−Sˆs0ps. (5b)

Solving (5) forg and p gives

p = (I−Fˆs0Sˆ0s)−1F(gs) + ˆFs0(S(ps)−Sˆs0psgs), (6a)

g = (I−Sˆs0Fˆs0)−1S(ps) + ˆSs0(F(gs)−Fˆs0gsps). (6b)

We assign the value ofpin equation (6a) to ps+1, and use it as an input for S,

ps+1 = (I−Fˆs0Sˆs0)−1F(gs) + ˆFs0(S(ps)−Sˆs0psgs).

(6)

We can then update the Jacobian ( ˆS0) in equation (6b) and assign the value of g togs+1:

gs+1 = (I−Sˆs+10 Fˆs0)−1S(ps+1) + ˆSs+10 (F(gs)−Fˆs0gsps+1). Combining this with Broyden’s method results in the following algorithm:

Algorithm 5.3 (Broyden’s method-2 Jacobians-Variant 2 (BrVar2)). 1. Startup:

Take a starting value po, ˆFo0 and ˆSo0;

computego =S(po) and p1 = (1−ω)po+ωF(go) (ω being a relaxa- tion parameter);

computeg1 =S(p1)

2. Make approximate Jacobian ˆS10 according to (3b);

set s= 1.

3. Loop until sufficiently converged a. Compute F(gs).

b. Make approximate Jacobian ˆFs0 according to (3a).

c. ps+1 = (I −Fˆs0Sˆs0)−1F(gs) + ˆFs0(S(ps)−Sˆs0psgs) d. Compute S(ps+1).

e. Make approximate Jacobian ˆSs+10 according to (3b).

f. gs+1= (I−Sˆs+10 Fˆs0)−1S(ps+1) + ˆSs+10 (F(gs)−Fˆs0gsps+1) g. s=s+ 1

Remark. It is obvious from the way we construct ˆFs0 and ˆSs0 and equations (5a) and (5b) that they respect the secant equations (4a) and (4b). For reasons of comparison with the other algorithms, we start from ˆFo0 = ˆSo0 = 0.

6. TEST-CASES

We propose two test-cases: the heat equation and a fluid-structure interac- tion problem. The convergence requirement will be defined based on a relative reduction of the residual, which is defined by kK(pkK(ps)k

o)k ≤ 10−5 for algorithms 4.1, 5.1 and 5.2, and on kFkF(g(gs)−psk

o)−pok ≤10−5 for algorithms 3.1 and 5.3, asK(ps) is not explicitly computed in these. The performance parameter is the number of times the functionF (or S) is called, as we have assumed this is the domi- nant cost of the algorithm; we call this value F C. We break off the algorithm ifF C >100 if it has not converged at that point.

(7)

6.1. Heat equation.

6.1.1. The model. The heat equation in one dimension without heat source is given by:

ρC∂T∂t∂x k∂T∂x= 0, (7)

whereT is temperature,ρ is density, C is heat capacity,kis thermal conduc- tivity. We assume we only have a linear solver to solve (7) with an imposed value of ρ, C and kthat may vary in space.

If we want to take dependency of ρ, C and k on T into account, we need to do this outside the solver for (7). They can be obtained by the following interpolation polynomials (based on values between 0 and 300C in [6]):

k = 7.0277·10−5T + 2.4388·10−2, (8a)

C = 4.3004·10−7T2+ 1.1850·10−5T+ 1.0048, (8b)

ρ = 5.3641·10−6T2−3.7809·10−3T+ 1.2781.

(8c)

We use a finite differencing scheme onN equidistant nodes for the solver of (7), with spacing ∆x, and an implicit time discretization with time step ∆t.

This gives the N linear equations (i= 1, . . . , N):

ρiCi(Tin+1Tin)−2∆x∆t2 (ki+1+ki)Ti+1n+1− (9)

−(ki+1+ 2ki+ki−1)Tin+1+ (ki+ki−1)Ti−1n+1= 0.

We will useg to denote the ensemble of discretized constantsρi, Ci and ki, which will be defined by a node-to-node relation withTi according to (8).

We only consider one time-step and solve (9) for the temperature at the next time-step. We write this asF(gs) =Ts+1. The parameters are computed for a given temperature, which is written as S(Ts) =gs.

As an initial condition we take the values at time-leveln, which is a uniform temperature ofTo= 150K, the right boundary condition is a temperature of 150K, while the left boundary condition is a function of time, defined as TL=To+T2osin πt10.

We use an initial relaxation factorω = 0.1.

6.1.2. Results and discussion. Various combinations of the ratio ∆x∆t2 were tried forN = 100 andN = 1000 (tables 1 and 2). For small values of ∆x∆t2 the three Broyden algorithms and the ISM performed equally well; Aitken’s method trailed behind. As the ratio ∆x∆t2 grows, the problem becomes harder to solve and the number of iterations becomes larger. The performance of the ISM initially degrades faster than the Broyden methods, which perform equally well. ForN = 1000 and ∆x∆t2 = 109 the ISM method catches up with the basic Broyden method, while the other three methods lag behind. For even higher

(8)

values of ∆x∆t2 convergence could not be obtained. It is thus seen that for this mildly non-linear problem there is not much difference between the methods, but the basic Broyden method is always the best. This is perhaps surprising as it only uses one Jacobian and thus less information than the Broyden variants.

Table 1. F Crequired for convergence for the heat equation. N= 100.

∆t

∆x2 Br BrVar1 BrVar2 ISM A∆2

10−10 3 3 3 3 4

10−9 3 3 3 3 4

10−8 3 3 3 3 4

10−7 3 3 3 3 5

10−6 3 3 3 4 5

10−5 3 3 3 4 5

10−4 3 3 3 4 5

10−3 3 3 3 4 5

10−2 4 4 4 5 6

10−1 4 4 4 5 6

1 5 5 5 6 7

10 6 6 6 8 8

102 6 6 6 8 8

103 6 6 6 8 9

104 6 6 6 8 9

105 6 6 6 7 9

106 6 6 6 7 9

107 6 6 6 7 9

108 6 6 6 7 9

109 6 6 6 7 9

1010 7 7 7 7 9

6.2. Fluid-structure interaction problem. In this problem both physical problems are vector-based and non-linear.

Fig. 1. One-dimensional flow in a flexible tube.

6.2.1. The model. We consider one-dimensional unsteady flow in an elastic tube (figure 1). (For a two-dimensional simulation, see [5]). The fluid is incom- pressible and inviscid and gravity is neglected. The governing equations are the conservation of mass and momentum which can be written in conservative

(9)

Table 2. F Crequired for convergence for the heat equation. N = 1000.

“div” = divergence or non-convergence: F C>100.

∆t

∆x2 Br BrVar1 BrVar2 ISM A∆2

10−10 3 3 3 3 4

10−9 3 3 3 3 4

10−8 3 3 3 3 4

10−7 3 3 3 3 5

10−6 3 3 3 4 5

10−5 3 3 3 4 5

10−4 3 3 3 4 5

10−3 3 3 3 4 5

10−2 4 4 4 5 6

10−1 4 4 4 5 6

1 5 5 5 6 7

10 6 6 6 8 8

102 6 6 6 8 8

103 6 6 6 8 9

104 6 6 6 8 9

105 6 6 6 8 9

106 6 6 6 8 9

107 6 6 6 7 9

108 6 6 6 7 9

109 8 9 9 8 9

1010 div div div div div

form as

∂g

∂t +∂(gu)

∂x = 0, (10a)

∂(gu)

∂t +∂(gu2)

∂x +1 ρ

∂(gp0)

∂xp0∂g

∂x

= 0 (10b)

to which suitable boundary conditions need to be added (velocity at the inlet and pressure at the outlet). g represents the cross-sectional area of the tube, u the velocity along the axis of the tube, ρ is the density of the fluid andp0 the static pressure.

We introduce the kinematic pressurep=p0/ρ. If the elastic wall of the tube has a constitutive law of the form g = g(p(x, t)), with its inertia neglected, then (10) can be rewritten in the following form:

∂p

∂t +u∂p

∂x+c2∂u

∂x = 0, (11a)

∂gu

∂t +∂gu2

∂x +∂gp

∂xp∂g

∂x = 0.

(11b)

The wave speed cis defined by

c2 = g

∂g

∂p

.

(10)

For the constitutive law of the structure we take the non-linear relation defined by:

g = go

p−2c2mk po−2c2o,mk

2

,

wherecmkis the Moens-Korteweg wave speedcmk =q2ρrEh

o

,EYoung’s mod- ulus,h the thickness of the tube wall and ro a transversal reference length of the tube. The length of the tube is L.

The velocity is imposed at the inlet as u = uo+uo

10sin2 πt

T

,

where uo is the initial velocity and T the period (= 1). At the outlet, a non-reflecting boundary condition is imposed:

∂u

∂t x=L

= 1 c

∂p

∂t x=L.

We use a one-dimensional mesh withN nodes of length ∆x. Fluid velocity, pressure and diameter of the tube are stored at the cell centers. Central discretization is used for all terms in (11), except for the convective term in (11b) where a first-order upwind scheme is used; pressure stabilization is added [5]. Pressure and velocity at the inlet, resp. outlet are linearly interpolated from nearby mesh-points. For the time-discretization we use backward Euler with time-step ∆t.

As in §6.1, we only consider one time-step and use those values as initial values for the iteration. We solve the discretized fluid problem for the pressure at the next time-step and a fixed geometry, which we write as F(gs) =ps+1. For the structural problem we computegfor a given pressure, which we write asS(ps) =gs.

The initial relaxation parameter is ω= 10−2, except for (κ, τ) = (10,10−2) and (10,10−3), where it is 10−4, resp. 10−5.

6.2.2. Results and discussion. Various combinations of parameters were tried forN = 100 and N = 1000. When the governing equations are non-dimensio- nalized it is seen that they are characterized by two non-dimensional param- eters κ = uco

o, being a dimensionless structural stiffness, and τ = uoL∆t, being a dimensionless time-step; co, uo, are the linearizing, i.e. reference, values of the wave speed, resp. fluid velocity. Various combinations of these were tested (tables 3 and 4). It is seen that a more flexible tube (lower value of κ) and a smaller time-step (lower value of τ) represent a more difficult problem to solve. The three Broyden methods steadily outperform the ISM and Aitken’s method but share similar performance for the ”easier” test-cases. When the problem becomes harder, the basic Broyden method appears the best, followed

(11)

by the second variant of Broyden’s method. Again this is somewhat surpris- ing as all the information of the two systems S and F is lumped into a single approximate Jacobian ˆKs0.

Table 3. F Crequired for convergence for the FSI problem. N = 100.

“div” = divergence or non-convergence: F C>100.

κ τ Br BrVar1 BrVar2 ISM A∆2

1000 10−1 3 3 3 5 4

1000 10−2 3 3 3 7 5

1000 10−3 5 5 5 14 8

100 10−1 4 4 4 7 5

100 10−2 6 6 6 17 10

100 10−3 9 9 9 div 16

10 10−1 6 6 6 28 11

10 10−2 9 9 9 div 18

10 10−3 37 38 37 div div

Table 4. F Crequired for convergence for the FSI problem. N = 1000.

“div” = divergence or non-convergence: F C>100.

κ τ Br BrVar1 BrVar2 ISM A∆2

1000 10−1 3 3 3 5 4

1000 10−2 3 3 3 7 5

1000 10−3 6 6 6 17 10

100 10−1 4 4 4 7 5

100 10−2 6 6 6 28 14

100 10−3 9 9 9 div 11

10 10−1 6 6 6 34 11

10 10−2 13 14 13 div 27

10 10−3 36 46 40 div div

7. CONCLUSION

The basic Broyden method was tested on interaction problems and com- pared with two proposed variants as well as the iterative substructuring method and the Aitken ∆2 method. For weakly non-linear problems and/or problems requiring few iterations to converge there was little to chose from between the three Broyden methods, although they were better than the remaining two. For more difficult problems not only did the three Broyden-based methods outperform the other two, but it appeared that the basic Broyden method came out best. This is surprising as, compared with the two variants, it only uses a single Jacobian to capture all the information that is available of the two constituent systems, while the other two use two approximate Jacobians, one for each system.

(12)

Acknowledgement. The first author kindly acknowledges the financial support from the “T. Popoviciu” Institute of Numerical Analysis (Romanian Academy), for attending the conference.

REFERENCES

[1] Broyden, C.G.,A class of methods for solving nonlinear simultaneous equations, Math.

Comp.,19, pp. 577–593, 1965.

[2] Deparis, S., Discacciati, M., Fourestey, G.and Quarteroni, A.,Fluid-structure algorithms based on Steklov-Poincar´e operators, Computer Methods in Applied Mechan- ics and Engineering,195/41-43, pp. 5797–5812, 2006.

[3] uttler, U.andWall, W.,Fixed-point fluid-structure interaction solvers with dynamic relaxation, Computational Mechanics 2008, DOI 10.1007/s00466-008-0255-5, to appear.

[4] Michler, C., Van Brummelen, E.H., de Borst, R., Error-amplification analysis of subiteration-preconditioned GMRES for fluid-structure interaction, Comput. Methods Appl. Mech. Engrg.,195, pp. 2124–2148, 2006.

[5] Vierendeels, J., Lanoye, L., Degroote, J. and Verdonck, P., Implicit coupling of partitioned fluid-structure interaction problems with reduced order models, Comput. &

Structures,85, pp. 970–976, 2007.

[6] http://www.engineeringtoolbox.com/air-properties-d 156.html, April 2008.

Received by the editors: May 21, 2008.

Referințe

DOCUMENTE SIMILARE

Using more precise majorizing sequences than before [1], [8], and under the same computational cost, we provide a finer semilocal convergence analysis of Newton’s method in

The elementary spline functions will be introduced here as an approximate solution of a Cauchy problem regarding linear differential equations, using the method described in [3]..

In this paper we give a general method to approximate the set of all efficient solutions and the set of all weakly-efficient solutions for a multiple criteria optimization

In this paper is presented a bicubic spline collocation method for the numerical approximation of the solution of Dirichlet problem for the Poisson’s equation.. The

Abstract. Using the idea of the least squares method, a nonlinear two point boundary value problem is transformed into an optimal control problem. For solving the optimal

following the convergence of the Chebyshev method attached to this problem and we shall apply the results obtained for the approximation of an eigenvalue and of

In this paper we present a method for the numerical solu- tion of a mod.el problem õf ttre two dimentional Self-Adjoint second order.. elliptic partial

Further, the synthesized barium and nickel ferrite nanoparticles were used as adsorbent to remove heavy metal ions such as copper, cadmium and lead from simulated