Degroote, Joris

Demeester, Toon

van Brummelen, E. Harald

Ghent University Academic Bibliography

VIII International Conference on Computational Methods for Coupled Problems in Science and EngineeringCOUPLED PROBLEMS 2019E. On˜ate, M. Papadrakakis and B. Schrefler (Eds)COMBINING A LEAST-SQUARES APPROXIMATEJACOBIAN WITH AN ANALYTICAL MODEL TO COUPLEA FLOW SOLVER WITH FREE SURFACE POSITIONUPDATESToon Demeester1,∗, E. Harald van Brummelen2 and Joris Degroote1,31 Department of Flow, Heat and Combustion Mechanics, Ghent UniversitySint-Pietersnieuwstraat 41, 9000 Ghent, Belgium, toon.demeester@ugent.be2 Department of Multiscale Engineering Fluid Dynamics, Eindhoven University of TechnologyPO Box 513-5600 MB Eindhoven, the Netherlands3 Flanders Make, BelgiumKey words: Free Surface Flow, Quasi-Newton, Fitting Method, Surrogate ModelAbstract. This paper presents a new quasi-Newton method suitable for systems thatcan be solved with a black-box solver for which a cheap surrogate model is available.In order to have fast convergence, the approximate Jacobian consists of two differentcontribution: a full rank surrogate model of the system is combined with a low rankleast-squares model based on known input-output pairs of the system. It is then shownhow this method can be used to solve 2D steady free surface flows with a black-boxflow solver. The inviscid flow over a ramp is calculated for supercritical and subcriticalconditions. For both simulations the quasi-Newton iterations converge exponentially andthe results match the analytical predictions accurately.1 INTRODUCTIONQuasi-Newton methods can be used to solve coupled problems in a partitioned way,using existing solvers for each system. A common application is fluid-structure interaction.In this paper, a new approach is presented to construct the approximate Jacobian of a(non-linear) system that can be solved with a black-box for which also a cheap surrogatemodel is available. To reduce the number of quasi-Newton iterations, two approximationsare added: a full rank surrogate model of the system is combined with a low rank modelwhich is constructed with known input-output pairs of the system using a least-squarestechnique.The method was developed in [1] for numerically calculating 2D steady free surface flowof water and air as encountered in marine engineering. The free surface is represented1Toon Demeester, E. Harald van Brummelen and Joris Degrootewith a fitting approach, i.e. the computational grid is aligned with it. As the air phasecan be neglected for the envisioned flows, the free surface becomes a domain boundarywhose position depends on kinematic and dynamic boundary conditions. Both the flowfield and the free surface position are unknown, so that this free-boundary problem mustbe solved iteratively. This is done by reformulating the problem as a root-finding problemand solving it with the new quasi-Newton approach.In Section 2 the new method to construct the approximate Jacobian is explained fora general system. The method is applied to the steady free surface problem in Section 3:first some details are given about the surrogate model for the flow solver, then the free-boundary problem is given in more detail and a solution method derived based on quasi-Newton iterations. In Section 4 the performance of the method is demonstrated bycalculating the inviscid free surface flow over a ramp at different Froude numbers.2 QUASI-NEWTON METHOD WITH LEAST-SQUARES AND SURRO-GATEA general black-box system S(x) = y is considered, with x,y ∈ Rn×1 respectivelythe input and output of S. The input vector x for which S(x) = 0 (i.e. the root ofS) can be found using quasi-Newton iterations. With superscript m the iteration index,∆xm = xm+1 − xm the change in input, and J ∈ Rn×n an approximate Jacobian of S, anew value xm+1 can be calculated using the well-known formulaJ∆xm = −ym = −S(xm). (1)The iterations start with a guess x0 and continue until ‖ym‖ drops below a small constantε. It is assumed that S is expensive to evaluate, so the number of quasi-Newton iterationsmust be as low as possible. For this purpose a Jacobian is constructed that consists of twoparts which both have advantages and limitations. By combining them, a more robustand effective Jacobian is obtained.The first approximation is a full rank surrogate model J sur. It can be based on e.g.analytical considerations or a low fidelity model of the system, depending on the problemat hand. It must be of full rank to ensure that Eq. 1 has a unique solution ∆xm. Thebetter J sur approximates the behavior of S, the faster the quasi-Newton iterations willconverge. However, if the prediction of J sur is inaccurate for certain modes, these modeswill slow down or even destabilize the iterations.The second approximation J IQN is constructed using the IQN-ILS algorithm by De-groote et al. [2], which was developed originally for doing partitioned fluid-structureinteraction with black-box solvers. J IQN is a low rank approximation which is built fromknown input-output pairs of the system. Inputs and outputs of previous iterations arecollected in the matricesV m =[∆xm−1 · · · ∆x0] , (2)Wm =[∆ym−1 · · · ∆y0] . (3)2Toon Demeester, E. Harald van Brummelen and Joris DegrooteThe Jacobian is then constructed with a least-squares technique:JmIQN = WmRm−1QmT(4)withV m = QmRm (5)the economy-size QR-decomposition. With each iteration JmIQN grows in rank and im-proves. The modes which are picked up by J IQN are those that hinder convergence of thequasi-Newton iterations, which are precisely those modes for which J sur is inaccurate.JmIQN only affects the part of x which is in range(Vm) = range(Qm). As these modescan be extracted using QmQmT , the combined Jacobian can be written asJm = JmIQN + J sur(In −QmQmT)(6)with In the identity matrix. In a nutshell: most of the work is done by the full ranksurrogate model J sur, but the low rank least-squares based J IQN is crucial to stabilizeand accelerate convergence of the iterations when J sur is inaccurate for some modes. Themethod is accordingly called Quasi-Newton method with Least-Squares and Surrogate andan overview of it is given in Algorithm 1.Algorithm 1 Quasi-Newton method with Least-Squares and Surrogate.1: m = 02: y0 = S (x0)3: while ‖ym‖ > ε do4: if m > 0 then5: update V m, Wm . Eqs. (2) and (3)6: QR-decomposition V m = QmRm7: end if8: construct Jm . Eq. (6)9: solve Jm∆xm = −ym10: xm+1 = xm + ∆xm11: m = m+ 112: ym = S (xm)13: end whileNote that if a surrogate model is available for the system’s inverse S−1, the quasi-Newton iteration of Eq. (1) can be rewritten as ∆xm = −J−1mym, avoiding the need tosolve a linear system. In addition the least-squares Jacobian J−1mIQN does not have to beconstructed explicitly, avoiding matrix-matrix products.3Toon Demeester, E. Harald van Brummelen and Joris Degroote3 APPLICATION TO STEADY FREE SURFACE FLOWBefore the free-boundary problem is solved with the new quasi-Newton technique, afree surface perturbation analysis is discussed in Section 3.1. Not only does this analysisprovide the basis for constructing the flow solver’s surrogate model, it also shows thatadditional conditions are needed to obtain a unique free surface solution for subcriticalflows. In Section 3.2 the 2D steady free surface problem is discussed in more detail,then it is expressed as a free-boundary problem and finally a solution algorithm based onquasi-Newton iterations is outlined.3.1 Surrogate model from perturbation analysisIn [3] the 2D inviscid free surface flow over a flat bottom is investigated in detail. Thisflow is characterized by its Froude number Fr = U/√gh which is based on the averageflow velocity U and flow depth h. The steady flow solution of this case is a horizontal freesurface and a uniform velocity U . It is shown in [3] that when the free surface height η(x)is perturbed with a sine wave, the free surface pressure p(x) changes proportionally. Thisanalysis leads to a relation between free surface height perturbations H(k) and pressureperturbations P (k) in the wavenumber domain:P (k) = L(k)H(k) (7)withL(k) = ρg(Fr2khtanh kh− 1). (8)Here ρ is the density, g the gravitational acceleration and k = 2pi/λ the wavenumber. Thisrelation between perturbations can be used to approximate the Jacobian of a general 2Dfree surface flow. An expression for J sur for uniform free surface grids is derived in [1]using the discrete Fourier transform. For stretched free surface grids J sur can be calculatedusing the convolution theorem as explained in [4].The behavior of free surface flow changes at Fr = 1. For supercritical flow (Fr > 1), thewave speed of surface gravity waves is smaller than the flow speed U for all wavenumbersso that disturbances can only travel downstream. For subcritical flow (Fr < 1) this is notthe case: waves below a certain wavenumber travel faster than the flow, so disturbancescan travel upstream. This change in behavior corresponds to a zero in L(k), Eq. (8):the wave with this wavenumber travels at the same speed as the flow and thus appearsstationary. This so-called steady gravity wave is a solution of the flow and therefore givesno change in pressure (L = 0). This is problematic: it means that for subcritical flow,this wave can be present with arbitrary phase and amplitude without influencing the freesurface pressure. As a consequence, two additional conditions will need to be added tofind a unique solution for the free-boundary problem.4Toon Demeester, E. Harald van Brummelen and Joris Degroote3.2 Solving the free-boundary problem with quasi-NewtonThe 2D immiscible steady flow of water and air as encountered in marine engineeringis considered as application. The free surface is represented with a fitting approach,which means that the computational grid is aligned with it. Due to the large densitydifference with water, the air phase can be neglected for these flows. The free surfacethen becomes a domain boundary whose position depends on kinematic and dynamicboundary conditions. For steady flow, the kinematic boundary condition (KBC) requiresthat the flow is tangential to the free surface. The dynamic boundary condition (DBC)requires continuity of the stresses at the free surface. This condition can be split intotangential and normal components which can be simplified. To fulfill the tangential DBCthe shear stresses at the free surface must be zero. For the normal DBC the pressuremust be constant along the free surface.The flow field and free surface position cannot be determined simultaneously by astandard CFD solver, so the free-boundary problem needs to be solved iteratively: in onestep the flow field is solved with a given free surface position, in the other step a new freesurface position is calculated. The distribution of the free surface boundary conditionsover these two steps varies between different methods found in literature [5, 6, 7].In this paper the conditions are distributed in such a way that the free-boundaryproblem becomes a root-finding problem. The black-box system of this problem is theflow solver F which uses a free-slip wall at the free surface boundary—this way the KBCand tangential DBC are satisfied in the flow solver step. The input and output of F arerespectively the free surface height η ∈ Rn×1 and the free surface pressure p ∈ Rn×1, sothat p = F(η). The normal DBC requires that p = pcst1 with 1 the all-ones vector. Thiscannot yet be formulated as a root-finding problem as pcst is unknown. It follows fromthe perturbation analysis in Section 3.1 that changing the value of pcst corresponds to achange of the average free surface height. Instead of specifying pcst it is more convenientto impose e.g. the free surface inlet height η(0) = h to fix the average free surface height.pcst is removed from the normal DBC by subtracting the average value (or constant mode)from the pressure. The flow solver F and pressure p after removing the average value aredenoted as F∅ and p∅, so that the normal DBC becomesF∅(η) = p∅ = 0. (9)This is a root-finding problem for a black-box system F∅ with input η and output p∅,which can be solved with quasi-Newton iterations. The inlet height condition must beadded to the system for a unique solution. With Jm the approximate Jacobian of F∅,the iterative scheme becomes: {Jm ∆ηm = −pm∅∆ηm(0) = h− ηm(0) (10)This system has n+ 1 equations for n unknowns and must therefore be solved with least-squares. The system is still well-determined though, because the rank of J is n − 1 due5Toon Demeester, E. Harald van Brummelen and Joris Degrooteto subtraction of the average value from p.In Section 3.1 it is explained that for subcritical flow two additional conditions need tobe added to ensure that the free surface solution is unique. A physically relevant optionis to require that the free surface is flat upstream. This is enforced by the conditionsη(i) = h with i = i1, i2 (11)where x(i1) and x(i2) are points close to the domain inlet x(0). Adding these conditions,the system that has to be solved in the quasi-Newton iterations becomes:{Jm ∆ηm = −pm∅∆ηm(i) = h− ηm(i) with i = 0, i1, i2(12)The Free Surface Quasi-Newton method with Least-Squares and Surrogate (FreQ-LeSS) issummarized in Algorithm 2.Algorithm 2 Free Surface Quasi-Newton method with Least-Squares and Surrogate.1: m = 02: p0∅ = F∅ (η0)3: while∥∥pm∅ ∥∥ > ε do4: if m > 0 then5: update V m, Wm . Eqs. (2) and (3)6: QR-decomposition V m = QmRm7: end if8: construct Jm . Eq. (6)9: solve system for ∆ηm . Fr > 1: Eq. (10) | Fr < 1: Eq. (12)10: ηm+1 = ηm + ∆ηm11: m = m+ 112: pm∅ = F∅ (ηm)13: end while4 INVISCID FREE SURFACE FLOW OVER A RAMPThe free surface method is demonstrated by calculating the inviscid free surface flowover a ramp as shown in Fig. 1, using the same geometry and flow parameters as Muzaferijaand Peric´ [5]. The total mechanical energy E (in J/kg) is conserved between inlet andoutlet of the domain. Assuming a uniform flow velocity in each section, it can be shownthatEi =U2i2+ ghi + gyb,i = cst (13)with Ui the flow velocity, hi the depth of the flow and yb,i the vertical position of thebottom. Combining Eq. (13) with the continuity equation gives analytical expressions6Toon Demeester, E. Harald van Brummelen and Joris DegrooteFigure 1: Definition of free surface flow over a ramp [5].−5.0 −2.5 0.0 2.5 5.0x [m]0.00.10.20.3y[m]Fr = 1.92−2 0 2 4x [m]−0.03−0.02−0.010.00y[m]Fr = 0.32Figure 2: Free surface for supercritical and subcritical flow calculations. The values are plotted withrespect to the inlet height. The ramp is located between x = 0 m and x = 1 m.for the outlet water depth h2 and outlet velocity U2 for known inlet conditions and rampgeometry.Two cases are simulated: a supercritical and a subcritical flow. For both cases theramp is 0.2 m high, h1 = 1 m, density ρ = 1 kg/m3 and gravity g = 9.81 m/s2. Thesupercritical flow has U1 = 1 m/s (Fr = 1.92), the subcritical flow U1 = 6 m/s (Fr = 0.32).In the flow solver the bottom and free surface are modeled as free-slip walls. At theoutlet a hydrostatic pressure is imposed, at the inlet a uniform velocity U1. For bothsimulations, the initial free surface is horizontal and corresponds to the inlet height.Structured quadrilateral grids are used, with 100 cells in the y-direction and 50 cells permeter in the x-direction. As the free surface grid is uniform, J sur is constructed based onEq. (8) and the discrete Fourier transform, as described in [1].The free surface for both cases is shown in Fig 2. The supercritical free surface is asmooth curve, while the subcritical free surface has the characteristic wave train down-stream of the obstacle. To avoid wave reflections at the outlet, a wave damping zone is7Toon Demeester, E. Harald van Brummelen and Joris DegrooteTable 1: Analytical result for outlet depth h2 and error in numerical result.Fr h2 (analytical) h2 (numerical)1.92 1.089730 m −0.008 %0.32 0.763544 m +0.127 %0 10 20 30 40coupling iteration m0.00010.0010.010.11∥ ∥ pm ∅∥ ∥ /∥ ∥p0 ∅∥ ∥Fr = 1.92Fr = 0.32Figure 3: Normalized pressure residuals of the quasi-Newton iterations.implemented near the outlet using momentum source terms [8]. In Table 1 the outletdepth of the flow h2 is given: the analytical result based on Eq. (13) is compared tothe result of the simulations, giving excellent results. The subcritical result is a bit lessaccurate, most likely because of the wave train and the damping zone.Fig. 3 shows the residuals for both cases: they converge exponentially and in a lownumber of iterations. The subcritical case takes longer: this is due to the presence of thesteady gravity wave, which makes the steady free surface problem –and the root-findingproblem– much harder to solve.5 CONCLUSIONSThis paper presents a new quasi-Newton method suitable for systems that can be solvedwith a black-box solver for which a cheap surrogate model is available. In order to havefast convergence, the approximate Jacobian consists of two different contribution: a fullrank surrogate model of the system is combined with a low rank least-squares model basedon known input-output pairs of the system. It is then shown how this method can be usedto solve 2D steady free surface flows with a black-box flow solver. The inviscid flow overa ramp is calculated for supercritical and subcritical conditions. For both simulationsthe quasi-Newton iterations converge exponentially and the results match the analyticalpredictions accurately.8Toon Demeester, E. Harald van Brummelen and Joris DegrooteREFERENCES[1] T. Demeester, E. H. van Brummelen, and J. Degroote, “An efficient quasi-Newtonmethod for 2D steady free surface flow.” Manuscript submitted for publication.[2] J. Degroote, K.-J. Bathe, and J. Vierendeels, “Performance of a new partitionedprocedure versus a monolithic procedure in fluid–structure interaction,” Computers& Structures, vol. 87, no. 11, pp. 793–801, 2009.[3] T. Demeester, J. Degroote, and J. Vierendeels, “Stability analysis of a partitionediterative method for steady free surface flow,” Journal of Computational Physics,vol. 354, pp. 387–392, 2018.[4] T. Demeester, E. H. van Brummelen, and J. Degroote, “Extension of a fast methodfor 2D steady free surface flow to stretched surface grids,” in Proceedings of the 8thInternational Conference on Computational Methods in Marine Engineering, 2019.[5] S. Muzaferija and M. Peric´, “Computation of free-surface flows using the finite-volumemethod and moving grids,” Numerical Heat Transfer, vol. 32, no. 4, pp. 369–384, 1997.[6] G. Tzabiras, “A numerical investigation of 2D, steady free surface flows,” InternationalJournal for Numerical Methods in Fluids, vol. 25, no. 5, pp. 567–598, 1997.[7] E. H. van Brummelen, H. C. Raven, and B. Koren, “Efficient numerical solution ofsteady free-surface Navier–Stokes flow,” Journal of Computational Physics, vol. 174,no. 1, pp. 120–137, 2001.[8] R. Peric´ and M. Abdel-Maksoud, “Analytical prediction of reflection coefficients forwave absorbing layers in flow simulations of regular free-surface waves,” Ocean Engi-neering, vol. 147, pp. 132–147, 2018.9

Combining a least-squares approximate Jacobian with an analytical model to couple a flow solver with free surface position updates

https://biblio.ugent.be/publication/8624068/file/8624069.pdf

Combining a least-squares approximate Jacobian with an analytical model to couple a flow solver with free surface position updates

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography