Acceleration of fluid-structure interaction procedures by anticipatory coupling

Acceleration of fluid-structure interaction procedures by anticipatory coupling VII International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2017 M. Papadrakakis, E. O˜nate and B. Schrefler (Eds) ACCELERATION OF FLUID-STRUCTURE INTERACTION PROCEDURES BY ANTICIPATORY COUPLING J.H. Seubers∗and A.E.P. Veldman† ∗†Computational Mechanics and Numerical Mathematics University of Groningen Nijenborgh 9, 9747 AG Groningen, The Netherlands web page: http://www.rug.nl/fmns-research/cmnm e-mail: ∗h.seub[email protected], †a.e.p.v[email protected] Key words: Fluid-solid Interaction, Numerical Method, Strong Coupling, Added Mass Abstract. Simulating the hydrodynamics of floating structures using a two-way partitioned coupling poses a major challenge when the coupling between the fluid and the structure is strong. The incompressibility of the fluid plays an important role, and leads to strong coupling when the ratio of so-called added mass to structural mass is considerate. Existing fluid-structure interaction procedures become less efficient in such cases, and can even become unstable. This paper proposes a coupling method that deals with the added-mass effect by anticipation, and remains stable and efficient at all times. 1 INTRODUCTION Traditionally, multi-physics problems are classified as ‘strongly’ or ‘weakly’ interacting problems. From a physical perspective, the interaction is called weak if one subsystem dominates the behaviour of the coupled problem, and it is called strong if more than one of the subsystems ‘equally contribute to the interaction’ [1] or ‘have an equal say’ [2]. So the physical interaction strength is a scale running from a one-way hierarchy between systems to a two-way complementary interaction. An example of a hierarchy is the case of a very light particle (e.g. a ping-pong ball) in a large water wave: the motion of the ball follows completely from the motion of the wave. The wave is not affected by the presence of the ball. The exact opposite hierarchy occurs for heavy objects (e.g. a mammoth tanker) in quiet water: the flow of the water is completely determined by the motion of the ship. Of course, many real situations are somewhere in between these asymptotic cases, with two-way interaction between the subsystems: some feedback occurs from the water to the ship or from the particles to the water. The more feedback, the stronger the interaction. In hydrodynamic applications with moving structures, a major factor affecting the interaction strength is the ratio of the added mass of fluid to the structural mass. In 1 558 J.H. Seubers and A.E.P. Veldman the traditional formulation, where the fluid loads are imposed on the structure and the structural motions imposed on the fluid, higher added mass ratios increase the interaction strength. This effect makes simulation by traditional coupling methods of slender structures in large waves computationally expensive. The objective of the proposed method is to reduce the computation time for such applications. Section 2 provides the motivation and physical background of fluid-structure coupling problems in marine hydrodynamics. The mathematical model of this problem is explained in section 3, which provides the necessary ingredients to analyze the coupling method. The new coupling method is introduced in section 4, where it is compared to existing methods. The properties of the new method are analyzed in section 5, and the results of some numerical experiments are discussed in section 6, leading to the conclusions in section 7. 2 PHYSICAL MODEL Interactive simulations are important for predicting the behaviour of moored or freefloating ships or platforms in different operating conditions. On deck various operations may be performed that affect the load or inertia distribution of the vessel. The vessel responds not only to incoming waves but to the flow caused by its own motions as well: wave slamming, launching, green water events. The inertia of the water mass involved in these interacting flows is important for predicting the forces on the vessel. In other words, the inertia is an important feedback mechanism that leads to a strong coupling between the flow and the vessel motion. The ship or platform, which will be referred to as the structure, is modelled as a rigid body with elastic mooring lines. The structure can have an arbitrary shape and can perform large but finite translations and rotations in three dimensions. It cannot deform or change in volume. The water is modelled as an incompressible, viscous fluid with a free surface. Although the inviscid flow behaviour dominates the coupling, the vorticity and viscosity are included in order to show that the story remains essentially the same. The air flow is not modelled, a vacuum takes its place instead. The interaction is modelled by conservation of momentum and geometric compatibility between the structure and fluid surface. Note that the topology and position of the fluidstructure interface can change in time. 3 MATHEMATICAL MODEL Because of its flexibility, the partitioned approach will be adopted in this work. The partitioning cuts the system into two parts, a fluid subsystem (subscript f) and a structure subsystem (subscript s). The two subsystems with appropriate boundary conditions are represented as dynamical systems in a state-space representation, governed by the mass- 2 559 J.H. Seubers and A.E.P. Veldman spring and Navier-Stokes equations respectively. M¨x +Kx=B T sfsρΩ˙u +Gp+ρC(u)u−µLu=B T fff(1a) ys=B s¨x y f=B f˙u (1b) where [x,u,p]Tare the internal states of the fluid-structure system, [yf,ys]Tare the motions at the component boundaries, and [ff,fs]Tare the distributed loads at these boundaries. The domain where these variables live may deform over the time interval [t,t + ∆t], see fig. 1. t x y internal fluid state u,p external interface interface states yf,ys and forces ff,fs internal structure state x Figure 1: Model interface in space-time. This mathematical model is not yet complete since the forces are not given. These are determined implicitly by two coupling criteria. The kinematic criterion requires that the motions on both sides of the fluid-structure interface are the same, δy:= yf(t)−ys(t)=0.(2a) The dynamic criterion expresses the balance of forces over the fluid-structure interface, f:= ff(t)+fs(t)=0.(2b) Since the interaction is concerned with the variables that live on the interface, the internal states [x,u,p]Tare eliminated from the system by linearizing and substituting (1a) into (1b). This will lead to two operators that give the motions yin terms of the loads f, the so-called Dirichlet-to-Neumann (DtN) operators Afand As. yf(t)=y0 f(t)+A f(t)∗ff(t),(3a) ys(t)=y0 s(t)+A s(t)∗fs(t).(3b) 3 560 J.H. Seubers and A.E.P. Veldman Together with the unloaded motions y0 f,s, these DtN operators completely describe the response of both subsystems to any load. Therefore, the difficulties of the interaction can be found by studying the properties of the DtN. It is easier to derive the DtN in the Laplace domain, where the time derivatives can be manipulated algebraically. To show how this is done, some simplified models are considered first. 3.1 Response of simplified models Consider a cylinder on a spring, moving horizontally in a quiescent potential flow (fig. 2). The cylinders response is found by considering the equations of motion in the Laplace domain: Figure 2: Spring-fixed cylinder moving in infinite potential flow s2mˆx+ kˆx=ˆ fs+ m(sx0+˙x0), ˆys=s2ˆx−(sx0+˙x0). where m is the cylinder mass and k is the spring stiffness. This yields the response ˆys=s2m s2m+k−1(sx0+˙x0)+s2 s2m+kˆ fs.(4) The DtN operator is recognized as ˆ As=s2 s2m+k, which represents the acceleration of the cylinder due to an impulsive force at t= 0. The acceleration due to any other force can be found by convolution. In particular, the instantaneous acceleration due to a step force f0can be found from the initial value theorem, ys(t=0) = f0lim s→∞ ˆ As(s)=f0 m.(5) The fluid response is simply given by the added mass force. In summary, this simple interaction problem is governed by the two DtN operators ˆ As=s2 s2m+k,ˆ Af=1 ma .(6) 4 561 J.H. Seubers and A.E.P. Veldman 3.2 Interactive response Now recall the physical description of the ping-pong ball and the mammoth tanker. Since the feedback from the ping-pong ball on the water is small, we create an asymptotic expansion for the fluid motion ˆy f(s) starting from the unforced fluid motion ˆy 0 f(s), ˆy f=ˆy 0 f+ˆ Afˆ A−1 sˆy 0 s−ˆy 0 f+O(ˆ2 f).(7) On the other hand, the asymptotic expansion for the mammoth tanker will start with the unforced vessel motion ˆy 0s(s), since the feedback from the water is small. ˆy s=ˆy 0 s+ˆ Asˆ A−1 fˆy 0 f−ˆy 0 s+O(ˆ2 s),(8) where the feedback strengths ˆf,ˆsare measured by the disturbance of the motion ˆf(s)=ˆ Af(s)ˆ A−1 s(s),ˆs(s)=ˆ As(s)ˆ A−1 f(s).(9) At most one of these asymptotic expansions (7) or (8) will converge for a given problem, since ˆf<1 implies ˆs>1. Supposing that (8) converges, the motion will be given by ˆy s=ˆy 0 s+∞  i=1 ˆ Asˆ A−1 fiδˆy 0=ˆy 0 s+I−ˆ Asˆ A−1 f−1δˆy 0(10) In particular, this provides the interactive motion of the simplified models ˆys=ˆy0 s+∞  i=1 s2ma s2m+ki δˆy0=ˆy0 s+s2m+k s2(m −ma)+kδˆy0.(11) Only for m >ma, the roots of the denominator are in the left half plane, hence an oscillatory solution bounded by the initial disturbance δy0exists. In more complex cases, it could happen that neither expansion converges. In that case, both subsystems contribute equally: the physical interaction is strong. Therefore it makes sense to define the physical interaction strength as a product of the feedback strengths: Definition 1. The physical interaction strength κof the closed system (eqs. (1a), (1b), (2a) and (2b)) is a number between one and infinity, given by the initial value of the product of the feedback strengths, κ= lim s→∞ ˆf(s)ˆs(s) This definition can be seen as the sensitivity of the responses for t→0. In the simplified scalar model (section 3.1) these sensitivities are the added mass ratio and its reciprocal, lim s→∞ ˆs=ma m,lim s→∞ ˆf=m ma(12) 5 562 J.H. Seubers and A.E.P. Veldman Generalizing this to systems, the sensitivities are the maximal and minimal eigenvalues of the matrix AsA−1 f. Note that κis also the condition number of this matrix. Its eigenvalues can be interpreted physically as ‘directional’ added mass ratios, i.e. depending on the direction of the motion vector. When the interaction strength equals one, the coupled motion is simply a linear combination of the unforced motions yf= (1 −α)y0 f+αy0 s. In general however, this rarely occurs and interaction strengths may be higher. For rigid bodies floating in incompressible flow, it will be shown in section 5 that the interaction strength is still a function of mass ratios. But first, the performance of coupling algorithms will be directly related to the interaction strength in section 4. 4 NUMERICAL COUPLING METHODS The basic coupling methods are related to the asymptotic expansions in eqs. (7) and (8). In marine hydrodynamics, the expansion (8) dominated by the ship motion is most natural, and it converges provided that the added mass is smaller than the ship mass, ˆ As(s)ˆ A−1 f(s)<1 i.e.  mas2 ms2+k<1.(13) This coupling approach works well for weakly coupled problems, provided that the ship is indeed the dominant subsystem. If the ship is less dominant, it may be required to mix the approaches of the mammoth tanker and the ping-pong ball, using combinations of eqs. (7) and (8). Indeed modern domain decomposition approaches like FETI [3] are based on the difference between eqs. (7) and (8), δy=yf−ys. This difference is iteratively reduced by splitting it over the two domains and then enforcing the dynamic criterion eq. (2b), δyi+1 =δyi+Af+A sαA−1 f+ (1 −α)A−1 sδyi,(14) where iis the iteration index. An early precursor to this approach is the semi-inverse method by Le Balleur (1978), δyi+1 =δyi−αAf+A sδyi.(15) In both methods, a new force is estimated from δyand fed identically (but with opposite sign) to both subsystems to produce the new motions. These forces and motions are notated here as simple vectors, not as functions of time since eqs. (14) and (15) are steady state methods. An extension of the semi-inverse method that operates on time series of forces and motions is known as waveform relaxation [4]. In each iteration, both subsystems are integrated in time based on an estimated force series. The difference in the resulting motion series are multiplied by the relaxation parameter αto produce a new force estimate, exactly as in (15). FETI, semi-inverse and waveform relaxation methods can all be seen as local preconditioners for the kinematic criterion eq. (2a), see fig. 3. For a suitable range of α, such methods can deal with strongly coupled problems. However, their convergence (δy→0) slows down as the coupling strength increases. 6 563 J.H. Seubers and A.E.P. Veldman Figure 3: Preconditioning schemes: semi-inverse, FETI, waveform relaxation Figure 4: Extrapolation schemes: Gauss-Seidel, IQN, IBQN, manifold mapping Another family of methods for strongly coupled problems is based on the Gauss-Seidel approach (fig. 4). Instead of applying the coupling criteria (2) directly, the expansion (8) is modified by introducing extrapolation steps. The simplest type of extrapolation is under-relaxation, yi+1 =yi+αAsA−1 fδyi.(16) This method will again converge for suitable α, albeit slower when the interaction becomes stronger. In fact, even for the optimal choice in α, the performance of this method deteriorates linearly with the coupling strength. 4.1 Anticipatory coupling The quasi-simultaneous method [5] however, avoids the need for extrapolation or relaxation by reformulating the original problem to an equivalent one with reduced coupling strength. This is achieved by replacing the system of coupling conditions (2) by an equivalent system, yf+D iff=ys−Difs,(17a) ff+fs=0.(17b) Note that eq. (17a) is formed by eq. (2a) plus an arbitrary operator Ditimes eq. (2b). Since eq. (2a) was previously used as a boundary condition for the fluid, this simply amounts to a more general boundary condition. A suitable choice of the operator Di would contain some approximate physics of the structure. In the anticipatory coupling method, we choose Dias the instantaneous approximation of the structure: Di= lim s→∞ ˆ As.(18) This choice is motivated as follows. The approximation (18) is 7 564 J.H. Seubers and A.E.P. Veldman •easy to obtain, only the inertia properties such as the mass of the structure are needed. •exact at t= 0, so zero-stability of the subsystems implies zero-stability of the coupling. •physically consistent, no artificial physics are introduced into the problem. The difference between Diand Aswill produce a finite error in the solution at finite timesteps, but this error can be controlled by varying the timestep size, or resorting to any of the above families of coupling methods. The anticipatory algorithm thus takes the following steps: 1. Initialize structural motion yold sand loads fold sfrom previous timestep 2. Predict fluid velocities ˜u due to convection and diffusion 3. Move the geometry based on the structural motion yold s 4. Compute new fluid force fnew fin Poisson equation with anticipative condition (17a) ˜y f+D ifnew f=yold s−Difold s 5. Compute new structural response ynew swith dynamic condition (2b) fnew s=fnew f 6. Enforce the kinematic condition (2a) ynew f=ynew s(discarding ˜y f) 7. Correct fluid velocities uand update the free-surface position 8. Go to the next timestep To see the effect of the anticipative condition, it is illustrative to look at the simplified model from section 3.1 again. In this example, the modified Dirichlet-to-Neumann maps become ˆ Af(s)+D i=1 ma +d i,(19a) ˆ As(s)−Di=s2 s2m+k−di.(19b) Although these are only a scalar equation, it is clear that this affects the asymptotic expansion in general as ˆy s=ˆy 0 s+(ˆ As−Di)(ˆ Af+D i)−1δˆy 0+O(ˆ2 s),(20) hence the feedback strength becomes ˆs(s)= s2ma s2m+k 1−dim 1+d ima +O(s−2).(21) Therefore any choice 0 <d−1 i<2mmakes ˆs(∞)<1, and then the expansion (8) converges, at least for small enough time intervals. 8 565 J.H. Seubers and A.E.P. Veldman 5 ANALYSIS OF THE METHOD To extend the convergence result for the simple models to the Navier-Stokes and massspring model, we will take the following steps. •Obtain the DtN operators for the Navier-Stokes and mass-spring models •Choose Dibased on the instantaneous response from the DtN operators •Show that the anticipative scheme has a feedback of order s−1∼∆t 5.1 Linearized DtN operators The discrete mass-spring model in (1a) for the structure is transformed into the Laplace domain, (s2M + K)ˆx =B T sˆ fs+ M(sx0+˙x 0) (22a) ˆy s=s2Bsˆx −Bs(sx0+˙x 0).(22b) By eliminating the internal unknown ˆx by the same procedure used in eq. (4) the response of the structure is found ˆy s=ˆy 0 s+B sM+s−2K−1BT sˆ fs.(23) The treatment of the Navier-Stokes model in (1a) contains an additional step, since both ˆp and ˆu must be eliminated. Starting from the basic equations linearized around u0in the Laplace domain −GTˆu =0(24a) (sρΩ+ρC(u0)−µL)ˆu =B T fˆ ff−Gˆp +ρΩu0(24b) ˆy f=sBfˆu −Bfu0,(24c) the pressure is eliminated first by using the continuity equation (24a). Denoting the action of convection and diffusion by ˆ T−1=I+(sΩ)−1(C(u0)−νL), an analog for the pressure Poisson equation is found GTˆ TΩ−1Gˆp =G Tˆ TΩ−1BT fˆ ff+G Tˆ Tρu0.(25) The momentum (24b) and pressure (25) equations are then used to eliminate ˆu and ˆp from (24c). What remains is the response of the fluid to the imposed loads, ˆy f=1 ρBfˆ TΩ−1BT fˆ ff−Gˆp +B f(ˆ Tu0−u0) =1 ρBfI+ˆ TΩ−1∆ˆ TΩ−1BT fˆ ff+ˆy 0 f(26) 9 566