A Principled Approach to Design Using High Fidelity Fluid-Structure Interaction Simulations
AA Principled Approach to Design Using High FidelityFluid-Structure Interaction Simulations
Wensi Wu a, ∗ , Christophe Bonneville a, ∗ , Christopher Earls a,b a School of Civil & Environmental Engineering, Cornell University, Ithaca, NY 14850,United States b Center for Applied Mathematics, Cornell University, Ithaca, NY 14850, United States
Abstract
A high fidelity fluid-structure interaction simulation may require many days torun, on hundreds of cores. This poses a serious burden, both in terms of time andeconomic considerations, when repetitions of such simulations may be required( e.g. for the purpose of design optimization). In this paper we present strategiesbased on (constrained)
Bayesian optimization (BO) to alleviate this burden. BOis a numerical optimization technique based on Gaussian processes (GP) thatis able to efficiently (with minimal calls to the expensive FSI models) convergetowards some globally optimal design, as gauged using a black box objectivefunction. In this study we present a principled design evolution that movesfrom FSI model verification, through a series of
Bridge Simulations (bringingthe verification case incrementally closer to the application), in order that wemay identify material properties for an underwater, unmanned, autonomousvehicle (UUAV) sail plane. We are able to achieve fast convergence towards anoptimal design, using a small number of FSI simulations (a dozen at most), evenwhen selecting over several design parameters, and while respecting optimizationconstraints.
Keywords:
Bayesian Optimization, Fluid-Structure Interaction, GaussianProcess, Machine Learning, Design Optimization ∗ Equal contribution
Preprint submitted to Elsevier August 25, 2020 a r X i v : . [ c s . C E ] A ug . Introduction Fluid-structure interaction (FSI) analyses have been successfully employed inincreasing numbers of engineering application areas over the past decades. Afew examples of such efforts include FSI modeling of heart valves [1], predictionof aerodynamics flutter around an AGARD 445.6 wing [2], parachute inflationsimulation in turbulent supersonic flows [3], and slamming-induced structuralresponse on a large container ship [4]. While the idea of effectively employingFSI analyses in design is attractive because of the enhanced insight it offers( i.e. enabling study of the interplay of a solid deforming due to forcing froma fluid medium), correctly capturing the physical response of a multi-physicssystem, where each sub-system is described using a distinct set of governingequations, is a nontrivial endeavor. Particularly, in a partitioned FSI approach( i.e. where the analysis context is such that the fluid and structural systems areconsidered using separate numerical descriptions, with external coupling strate-gies imposed to ensure satisfaction of salient transmission conditions along thefluid-structure interface) [5, 6, 7], a very fine temporal and spatial discretizationis usually required to properly resolve the flow fields, and to correctly repro-duce the physical mechanisms along the fluid-structure interface. In addition,in order to mitigate instabilities arising from the presence of the so-called addedmass effects ( i.e. the two-way influence of the adjacent fluid mass contactingthe structure, which is especially pronounced in the context of incompressiblefluids) [8, 9, 10], implicit coupling strategies are required to properly satisfy theconservation equations along the FSI interface. Hence, performing FSI analysesis rather computationally expensive: a serious problem when such simulationsare required within the design cycle of some engineering artifact.This computational cost often arises as an obstacle for optimizing structuraldesigns subject to heavy fluid loads. Indeed, optimizing design parametersthrough random search, with the need to evaluate the design candidate withFSI analysis, at each new design point, would be very costly, both in time and2n computational resources. Furthermore, trying random guesses offers very lit-tle hope for effectively finding the optimal candidate, especially as the numberof design parameters increases. Standard numerical optimization methods ( e.g. gradient-based methods, etc. ) are typically not applicable in the present casebecause they require convexity in the design objective function, and further,usually require many objective function evaluations. Consequently in this pa-per, we present design optimization strategies based on Bayesian optimization(BO): a powerful optimization method that is well suited for optimizing blackbox objective functions that are potentially noisy. In addition, BO is typicallyable to achieve quick convergence to an admissible candidate with minimal ob-jective function evaluations. BO has been successfully applied to multiple fieldsof science and engineering, such as material design [11, 12], aerospace [13], biol-ogy [14], drug discovery [15], and hyperparameter tuning in data sciences [16].In the context of FSI, we particularly focus on BO with inequality constraints[17], which allows to handle more realistic design problems ( e.g. if we wantto reduce the drag of a hydrodynamic profile, we might want to do so whilemaintaining a minimal lift, maximum admissible structural deformation, etc. )In this paper, We are proposing a principled approach to design when high-fidelity FSI simulations are required as part of the design loop. We begin withthe identification of a verification problem that is “close” within the design spaceof our ultimate application (in our case, the Turek & Hron FSI3 benchmark [18]fits our needs). We subsequently incrementally extend the verification contextto bring it ever closer to our application, using what we call Bridge Simula-tions . We then apply our, now trusted, FSI simulation capability to treat thecase of composite material design for use in the sail plane of a UUAV, as ademonstration. If, when extending the proposed approach to other problems,the application in question does not have a suitable benchmark verification casein the literature, then one must resort to experimental validation of a simplifiedversion of the design context: this may be advisable in all cases, if feasible.3he outline of this manuscript is as follows: Section 2 provides details con-cerning the FSI solution framework used in the present work, along with someelements of Bayesian optimization theory; FSI validation results, as well asbridging simulation results, are presented in Section 3 and Section 4; Section 5provides a demonstration on how the proposed BO/FSI method can be appliedto the composite design for a UUAV sail plane; and Section 6 concludes the work.
2. Methods
A depiction of the computational domain for a representative fluid-structure in-teraction problem is presented in Figure 1. Within the computational domain,a structural cantilever member, denoted as Ω s , is enclosed in a fluid domain,denoted as Ω f . The fluid domain (driven by an inflow defined by u f ( t )) interactswith the structural domain through the moving boundary Γ. Mesh deformationwithin the FSI system is governed by a mapping function χ , which maps themesh motion from its reference configuration Σ to its deformed configurationΣ( t ). Figure 1: A representative computational domain within an FSI system
The present work adopts a partitioned approach to FSI: relying to two opensource analysis tools when treating the fluid and structural domains, respec-tively. The open source computational structural dynamics (CSD) finite elementsoftware used in this work is CU-BENs [19], and the open source computational4uid dynamics (CFD) finite volume software, Open Field Operation and Manip-ulation (OpenFOAM 1.6-ext) [20, 21, 22, 23], is used to treat the fluid domain.Information transfer between the CFD and CSD solvers is managed using acoupling library, denoted herein as the “
Coupler ” [24].
CU-BENs [19] is an open source structural modeling finite element solver, writ-ten in C and developed at Cornell University. CU-BENs exploits high orderstructural elements for use in fully nonlinear ( i.e. both material and geometriccontexts) transient dynamic analysis. The CU-BENs discrete Kirchhoff triangu-lar (DKT) shell element is used in the present work, when modelling the struc-tural responses within the FSI system. The Newmark implicit time integrationscheme, augmented with generalized- α method [25], for numerical stability pur-poses (artificial, numerical damping can improve solution behavior), is used fortime integration within our structural analyses. OpenFOAM version 1.6-ext is a comprehensive CFD solver based on the finitevolume method. In many ways, OpenFOAM functions as an object-orientedapplication programming interface (API) for C++. In this work, we employOpenFOAM 1.6-ext, since it provides a FEM decomposition framework [20, 26](also known as “mini-element” method), along with a Laplacian-based meshmotion equation technique [20, 27, 26], to govern the mesh motion within thefluid domain. We found the finite element vertex-based mesh motion solver toproduce superior numerical results, as compared with some of the other pop-ular mesh motion methods ( e.g. the spring analogy [28, 29] and its variants[30, 31, 32, 33]), as it guarantees boundedness in the Laplacian operator (theLaplace equation governs mesh vertex motion), even when the fluid cells ap-proach a degenerate state [34, 20]. This feature is vital for FSI simulations5nvolving large structural displacements; since such structural responses couldlead to considerable deformation and skewness of the attached fluid control vol-umes.
The coupling library (the Coupler), initially developed by the Navel SurfaceWarfare Center Carderock Division (NSWCCD) [24], was modified for use inthe present work. Due the extreme importance of the Coupler in FSI analyses,more detail is offered in what follows. In terms of functionality, the Couplerensures that the salient transmission conditions are satisfied along the fluid-structure interface. Also, the Coupler manages data transfer between the CSDand CFD models. To effect communication with the Coupler, solver interfaceclasses ( i.e.
CFD interface and CSD interface) are written for OpenFOAM1.6-ext and CU-BENs. Specifically, these C++ interface classes work with therespective CSD and CFD solvers to: prepare necessary data structures; adaptcross domain variables; and translate native input/output to a form that iscompatible with the Coupler. The interface classes are subsequently translatedto “Python-ized” modules using the Python wrapper generator: SWIG. Theresulting python files are then integrated into the Coupler. This allows the in-terface classes to communicate with the Coupler using pointers to memory; thuseliminating the need for writing and reading input/output files. A schematicrepresentation of our implicit partitioned FSI solver coupling framework is pre-sented in Figure 2.Within a given FSI simulation, it is the CFD solution that is most demand-ing, both in terms of time step size (time increment) and spatial discretization.The CSD solver can accommodate any time step size (due to its implicit timeintegration) and requires much coarser meshes (as a result of the high order el-ements). The orchestration of the FSI solutions, within each time step, followsthe general scheme: 6. OpenFOAM computes the hydrodynamic loading, given the structuraldeformation from the previous time step, using an
Arbitrary LagrangianEulerian (ALE) scheme.2. The CFD interface object takes in the resulting fluid fields and preparesnecessary data for the Coupler.3. The Coupler interpolates the hydrodynamic loading, and subsequentlymaps the loading onto the structural mesh.4. CU-BENs then computes the corresponding structural deformation.5. The CSD interface class takes in the resulting FSI boundary deformationand prepares necessary data for the Coupler.6. The Coupler checks for convergence by monitoring out of balance displace-ment and/or pressures between the fluid and structural domains: compar-ing to a set of user specified tolerances. If the response along the interfaceconverges, the simulation advances to the next time step. Otherwise, thepredicted displacements are passed to OpenFOAM, to update the meshconfiguration, and the process is repeated.
Figure 2: Implicit partitioned FSI coupling framework with OpenFOAM as the fluid solver,CU-BENs as the structural solver, and a Coupler that facilitates all necessary communicationsbetween the two solvers to ensure the fluid-structure interface conforms to the transmissionconditions
Nonconforming mesh projection [35] is utilized to guide the grid-to-grid meshmapping and updating along the fluid-structure interface. This technique al-7ows for an optimal choice of mesh, for resolving the physics within the respec-tive structural and fluid domains, as it provides flexibility for load and motiontransferring on the FSI interface, when the mesh sizes in each subsystem are in-congruous. Figure 3 provides a visual representation of nonconforming mesh ina FSI analysis. The structural domain is discretized with DKT shell elements;thus thickness is implied as a shell property (and not explicitly modeled). Insuch a context, the shell elements (green) are situated at the mid plane of thephysical structure, and enclosed by the fluid boundary (grey). Orthogonal pro-jection mapping is initiated at the beginning of a FSI simulation, to identifyappropriate coupling pairs within the structural and fluid node sets; thus en-abling mesh motion tracking throughout the analysis. It is typically the case inFSI analyses that the number of CFD nodes, occurring along the FSI interface,far exceeds the number of CSD nodes; thus leading to non-conforming mesh ge-ometries that complicate data transfer in the FSI transmission conditions. TheCoupler uses a standard approach when handling this case.
Figure 3: Nonconforming mesh: grey represents the fluid boundary that interact with thestructural elements, while green represents the structural DKT shell elements
Within our Coupler, inverse distance weighting interpolation (IDW) [36] is usedto convert hydrodynamic pressure, located at the center of each fluid patch (oneface of the finite volume), along the interface, into point loads, before they areprojected onto the structural mesh. The IDW procedure used in the currentwork, employs the following relation: F ( x ) = n (cid:88) i =1 ω i ( x ) F ( x i ) (cid:80) nj =1 ω j ( x ) , ω i ( x ) = 1 d ( x, x i ) α (1)8here x i denotes the spatial coordinate of cell center i , F ( x i ) represents the fluidpressure at x i , n is the total number of nodal points used in the interpolation, ω i stands the weighting factor, d ( x, x i ) corresponds to the distance between thenode x and cell center x i , α is the distance-decay parameter, and F ( x ) is theestimated point load at location x . Typically, α = 2 is used in the standardIDW. A simple representation of the fluid interface is shown in Figure 4 to helpvisualize the IDW interpolation process. Figure 4: Inverse distance weighting interpolation for converting fluid pressure at cell center(blue) to the point load (red)
Time step advancement, and solution convergence on the fluid-structure inter-face, are controlled by the Coupler using an iterative method: the IQN-ILSiteration scheme [37, 7]. The IQN-ILS is a quasi-Newton root finding method,within which Newton-Raphson iterations are used, within a given time step,to solve for the interface displacement R ( d ) = S ( F ( d )) = 0. In the equationmentioned, F () represents a CFD solver for a given displacement d , while, S ()is a CSD solver for a given loading resulting from the CFD solver. The Newton-Raphson iteration formulation within each time step is shown in Equation 2 and3 R k + R (cid:48) k ∆ d k = 0 (2) d k +1 = d k + ∆ d k (3)Combining Equation 2 and 3, a new expression for d k + is derived d k +1 = d k − ( δ R (cid:48) k δ d ) − R k , (4)9here R (cid:48) k represents the Jacobian of the interface displacement residual oper-ator. Moreover, the residual vector R k , which is evaluated at d k , is describedin Equation 5 R k = S ( F ( d k )) − d k = ˜d k +1 − d k . (5)Herein, ˜d k + represents the predicted FSI interface displacement at iteration k + 1, whereas d k + is the corrected FSI interface displacement at iteration k + 1. The Newton-Raphson iteration step is repeated until R k satisfies a userspecified convergence criterion (cid:13)(cid:13) R k (cid:13)(cid:13) ≤ (cid:15) , where (cid:15) is a convergence tolerance(standard Euclidean norm, implied).The procedure of IQN-ILS may seem easy enough when the exact Jacobianof R is known. However, in many cases, the Jacobians of the fluid and struc-tural domains are not readily available; as in the case of coupling black boxCFD/CSD pairs. In such cases, inverting R (cid:48) k , to compute ∆ d k , is no longera trivial task. Fortunately, we can circumvent the problem by constructing aset of basis vectors on R and ˜d , such that the solution to ( δ R (cid:48) k δ d ) − R k can beapproximated ( i.e. the inverse of the Jacobian does not need to be formed ex-plicitly). A sketch of this procedure follows next.When approximating the required Jacobian, a series of simple steps are used,such that their collective forms the required linearization within the IQN-ILSalgorithm. First, the difference of residual and interface displacement, fromthe previous iteration (denoted as superscript k −
1) and the current iteration(denoted as superscripted k ), are computed for iteration i = 0 , , , , ..., k − R k − = R k − R k − (6)∆ ˜d k +1 = ˜d k +1 − ˜d k (7)Each ∆ R i corresponds to a ∆ ˜d i +1 , and these vectors are stored as columns10ithin the matrices V k and W k V k = (cid:104) ∆ R k − ∆ R k − ... ∆ R ∆ R (cid:105) (8) W k = (cid:104) ∆ ˜d k ∆ ˜d k − ... ∆ ˜d ∆ ˜d (cid:105) (9)where the number of columns, determined by the number of iterations, q , islimited by the degrees-of-freedom, p , on the FSI interface, such that the sys-tem is over-determined. In the event that the number of iterations exceedsthe cardinality in degrees-of-freedom along the FSI interface, the rightmost col-umn of V k must be truncated, to preserve the over-determined characteristicof the system: this allows the system to be solved as a least square problem [38].Subsequently, a set of orthogonal basis vectors is constructed, where the valueof ∆ R k is approximated as a linear combination of the known ∆ R i , as:∆ R k ≈ k − (cid:88) i =0 α ki ∆ R i = V k α k (10)with α ki ∈ R q × . The matrix V k is then decomposed via economy QR decom-position, using Householder transformations V k = Q k R k (11)into an orthogonal matrix Q k ∈ R p × q , and an upper triangular matrix R k ∈ R q × q . Substituting Equation 11 into Equation 10 yields R k α k = ( Q k ) T ∆ R k , (12)where the coefficient vector α k in Equation 12 can be obtained by backwardsubstitution.Similarly, the corresponding ∆ ˜d k + can also be calculated as a linear com-bination of the ∆ ˜d i , from previous iterations, ∆˜d k + ≈ k − (cid:88) i =0 α ki ∆˜d i = W k α k . (13)11rom Equations 3, 5, 6 and 7, ∆ R k = ∆ ˜d k + − ∆ d k is derived, this leads tothe relation ∆d k = W k α k − ∆ R k . (14)Given that ∆d k = − ( δ R (cid:48) k δ d ) − R k and ∆ R k = − R k (assuming in an ideal case R k +1 = 0), the quasi-Newton Raphson iteration step in Equation 4 is thenre-expressed as d k +1 = d k + W k α k + R k . (15) We now discuss the second critical piece within our FSI design “puzzle”: Bayesianoptimization. Bayesian optimization (BO) is a powerful global optimizationframework that permits the minimization of complex, non-convex objectivefunctions [39, 40, 41]. Unlike gradient based optimization methods, BO doesnot require calculation of gradients to the objective function, and typically onlyrequires few evaluations of the objective function, in order to achieve satisfac-tory convergence towards the global minimum. Consequently, it is a methodwell suited for minimizing (noisy) black box functions that are expensive toevaluate ( e.g. experimental evaluations, computationally intensive evaluations, etc. ) The general optimization problem may be stated as follows:find ˆ f = min x ∈ Ω f ( x ) s.t. c k ( x ) ≤ λ k ( k ∈ [[1 , n c ]] ) (16)Where f is the objective function, ˆ f its global minimum, Ω the domain of defi-nition of f , and c k is a set of n c inequality constraints ( c k may be a black boxfunction that is evaluated at the same time as f ) [17].The main idea of BO is to approximate f with a Bayesian regression surro-gate that interpolates on a set of already known values of f . The Bayesianregression is vital, since it outputs “confidence intervals” on its prediction; thusproviding useful hints concerning the veracity of our interpolations of the ob-jective function. Additionally, this built-in uncertainty quantification on the12bjective function allows for a principled formulation for the requisite acqui-sition function ; used in guiding our search for a minimum towards the mostpromising regions of the design space. Evaluating f within the regions of designspace that are specified by the acquisition function provides new data that aresubsequently used to refine the Bayesian regression surrogate, and update theacquisition function prediction [41]. Iterating this process usually yields conver-gence with only a small number of objective function evaluations.In the following subsections, we first provide more details about Gaussian pro-cesses (GP) [42], which is the most widely used Bayesian regression methodin BO. Then, we provide more insights concerning acquisition functions andconstrained Bayesian optimization (cBO) [17]. In this brief introduction of GP regression, we employ the same notationsas in [43]. Let us consider a machine learning training dataset made of n input-output couples ( training data ): D = ( X, y ) = ( x i , y i ) i ∈ [[1 ,n ]] . x i is an input vector( e.g. the optimization/design parameters in the case of BO) and y i a collectionof output scalars ( e.g. the corresponding value of the objective function in thecase of BO). y i is assumed to be contaminated with Gaussian white noise (othernoise models can be accommodated) and modeled as follows: y i = f ( x i ) + (cid:15) i (cid:15) i ∼ N (0 , σ ) (17)In machine learning (of which GP regression is a type), the main goal is tolearn a representation of the unknown function f , given the available trainingdata, in order to make prediction for new (yet unobserved) input points x ∗ [44, 45, 46]. Assumptions must be made about the functional form of f , andin GP regression it is taken as a stochastic vector sampled from a multivariate13aussian distribution [42] : f ∼ GP (0 , k ( X, X )) (18)The covariance matrix k ( X, X ) depends on the input training data, and is com-puted using an arbitrary covariance kernel function. The kernel function, k , is adesign choice that reflects prior knowledge one may have about the nature of f (smoothness, linearity, etc. ) In BO, it is often hard to have any intuition abouthow the objective function behaves (especially in black-box optimization). Apopular choice, that we will adopt in the following, is the 5 / k ( x i , x j ; θ ) = θ Γ( ν )2 ν − (cid:18) √ νθ || x i − x j || (cid:19) ν K ν (cid:18) √ νθ || x i − x j || (cid:19) (19)where ν = 5 /
2, Γ is the standard gamma function, K ν is the modified Besselfunction of the 2 nd kind, and θ = ( θ , θ ) is a set of tunable hyperparameters.These hyperparameters can be learned directly from the data, by maximizingthe marginal likelihood: p ( y | X, θ ) = (cid:90) p ( y | f ) p ( f | X, θ ) df (20)The likelihood p ( y | f ) is an immediate consequence of Equation 17, and p ( f | X, θ )is the prior distribution defined in Equation 18. In the present case, the marginallikelihood is analytically tractable and exact Bayesian inference can be per-formed [42]. In order to make predictions at any location, a new input point x ∗ is appended to the training data X , and the noise variance is incorporatedwithin the upper left block of the covariance matrix. Equation 18 thereforebecomes: yf ∗ ∼ N (cid:32) , k ( X, X ) + σ I n k ( X, x ∗ ) k ( x ∗ , X ) k ( x ∗ , x ∗ ) (cid:33) (21)The joint distribution above can then be combined with the marginal likelihoodto find the predictive distribution, conditioned on the training data: p ( f ∗ | X, y, x ∗ ) = p ( y, f ∗ | X, x ∗ ) p ( y | X ) = N ( f ∗ | µ ∗ ( x ∗ ) , Σ ∗ ( x ∗ )) (22)14n the present case, the distribution over the predictive output, f ∗ , is Gaussian,and analytical expressions can be found for its mean value µ ∗ and variance Σ ∗ .More details on GPs and practical implementation considerations can be foundin the reference book from Rasmussen & Williams, Gaussian Processes for Ma-chine Learning [42]
Gaussian processes are a powerful tool for building a surrogate of the objectivefunction that we want to minimize. As a Bayesian method, it typically avoidsthe caveat of overfitting [47], and the predictive distribution offers insight onhow much is truly known about the objective function. To find out where theobjective function should be evaluated next, acquisition functions are used. Theacquisition function embodies a principled means to quantify, for any x ∈ Ω,how likely it is that evaluating f ( x ) will yield improvement towards the globalminimum of our objective function. Many types of acquisition functions can befound in the literature, but in the current work we adopt Expected Improvement(EI) : a popular method that is both analytically tractable and can easily handleconstrained BO [41, 17].A Bayesian optimization iteration loop must be initialized using a few eval-uations of the objective function (typically at random input locations). Amongall of the evaluations, the (currently known) minimum is denoted as ˜ f . For any x ∈ Ω, the potential improvement toward the global minimum, given ˜ f is: I ( x ) = max(0 , ˜ f − f ( x )) (23)where f is unknown, but approximated with a GP, so the expected improve-ment can be derived and expressed in term of the GP predictive mean and15ariance [41].EI( x ) = E p ( f ∗ | X,y,x ) [ I ( x )]= (cid:90) + ∞−∞ max(0 , ˜ f − f ∗ ) N ( f ∗ | µ ∗ ( x ) , Σ ∗ ( x )) df ∗ = ( ˜ f − µ ∗ ( x ))Φ( ˜ f | µ ∗ ( x ) , Σ ∗ ( x )) + Σ ∗ ( x ) N ( ˜ f | µ ∗ ( x ) , Σ ∗ ( x )) (24)Φ refers to the univariate normal cumulative distribution function (CDF). Theexpected improvement can now be maximized using standard optimization meth-ods ( e.g. BFGS, gradient descent, etc. ) [48], and the objective function can beevaluated at the new input x next = arg max EI( x ). Finally, the evaluated value f ( x next ) can be used to update the GP predictions, and if a better global min-imum candidate value is obtained, the value of ˜ f , as well. This process can beiterated until convergence (as gauged by a user specified tolerance) or until thebudget for evaluating the objective function has been exhausted. Note that eachcomplete BO iteration comprises two optimization sub-iterations: The first onefor maximizing the marginal likelihood when fitting the GP to the data, andthe second one for maximizing the acquisition function. As a result, each BOiteration is fairly computationally intensive (but also worth it, when the cost ofevaluating the objective function is itself significantly more expensive).In order to deal with constrained optimization context, expressed in Equation16, the standard expected improvement acquisition function in Equation 24must be modified [17]. To account for the fact that some inputs x ∈ Ω are sim-ply not feasible with respect to constraint c k ( x ) ≤ λ k , the modification resultsin: cI ( x ) = ∆( x ) I ( x ) (25)∆( x ) = c k ( x ) ≥ λ k c k ( x ) ≤ λ k (26)16ince c k may be a black box function, evaluated at the same time as the objectivefunction, it may be approximated with a GP as well. As a result, ∆( x ) isa random variable following a Bernoulli distribution of parameter ρ ( x ) [17].Furthermore: ρ ( x ) = p ( c k ( x ) ≤ λ k ) = (cid:90) λ k −∞ p ( c k ∗ | X, y, x ) dc k ∗ (27)Where p ( c k ∗ | X, y, x ) is the Gaussian predictive distribution of c k ∗ in Equation22 and conveniently ρ is a univariate Gaussian cumulative distribution function[17]. Finally, the expected constrained improvement becomes:cEI( x ) = E[∆( x ) I ( x )]= E B (1 ,ρ ( x )) [ ρ ( x )]E p ( f ∗ | X,y,x ) [ I ( x )]= ρ ( x )EI( x ) (28) To implement our Bayesian optimization codes, we use botorch [49], a dedicatedpython library. botorch relies on gpytorch [50] for performing all the relevantGP inferences, and both of these libraries are built on top of the deep learningframework pytorch [51].
3. FSI Validation
As part of our design approach, we carefully searched the literature for a canon-ical FSI verification problem that was similar enough to our ultimate designapplication ( i.e. fiber reinforced plastic UUAV sail plane), to be useful is build-ing or confidence in the veracity of our FSI modeling approach. To that end, theTurek & Hron’s benchmark cases [18] are considered here as we judge the accu-racy of our implicit partitioned FSI approach using CU-BENs and OpenFOAMsolvers. The general problem setup for our verification context is shown in Fig-ure 5. Within the structural domain, an elastic cantilever beam is attached toa fixed cylinder whose radius is 5 cm; the center of the cylinder is positioned at170.2 m, 0.2 m). The cantilever beam is 0.35 m long, 0.02 m thick, and 0.01 mwide.
Figure 5: Turek and Hron FSI3 case
No-slip boundary conditions are imposed at the top and bottom walls of thechannel, as well as the interface along the cylinder and the cantilever. Thepressure outflow condition along the right edge of the channel is prescribed toatmospheric pressure. The channel inflow condition is defined along the leftedge of the domain as a parabolic velocity profile, described in Equation 29,where U is the mean inflow velocity. Equation 30 describes the evolution of theinflow velocity as a function of time. u f (0 , y ) = 1 . U . y (0 . − y ) (29) u f ( t, , y ) = u f (0 , y )[1 − cos( πt )]2 if t < . u f (0 , y ) otherwise (30)The verification procedure sets out to verify three test cases ( i.e. the CSM3test, the CFD3 test, and the FSI3 test) detailed in [18], each test case eval-uates the behavior of the corresponding component response within our CU-BENs/OpenFOAM implicit FSI solver, so that comparisons can be made withthe ground truth data reported in [18]. Mesh refinement studies are undertakenat this point, so that we may know what spatiotemporal descretizations aresuitable as we move towards our ultimate design case.18 .1. CSM3 Test An elastic cantilever beam subjected to gravitational force, only, with g =2 m/s is carried out in the CSM3 test. The elastic beam has a Young’s modu-lus of 1.4 MPa, Poisson ratio of 0.4, and mass density of 1000 kg/m . ImplicitNewmark time integration scheme, coupled with Newton Raphson sub-iterationwithin each time step, is used to capture the geometric nonlinearity effectswithin the CU-BENs structural model. The displacements at the free end ofthe cantilever are recorded. Figure 6 and Table 1 presents the x– and y– tipdisplacements over the time span t = [8 ,
10] seconds. As shown, that the tipdisplacements tend to the reference values provided by Turek & Hron, as thenumber of shell elements increases; negligible differences among the displace-ments are observed.
20 shell elements40 shell elements80 shell elements
20 shell elements40 shell elements80 shell elements
Figure 6: x– and y– displacements at the free end of the cantilever at time t = [8 ,
10] seconds: δt = 0 .
005 sec
Number of shells x– displacement [centerline ± amplitude] frequency [x– disp] y– displacement [centerline ± amplitude] frequency [y– disp]20 -15.602 ± ± ± ± ± ± ± ± Table 1: Peak x– and y– displacement at the free end of the cantilever at time t = [8 , δt = 0 .
005 sec .2. CFD3 Test Next, the fluid dynamic aspects of the coupled system are verified using theCFD3 test case, within which the cantilever is modeled as a rigid beam. Aparabolic inflow described in Equation 30, with 2 m/s mean inflow velocity,is prescribed along the left boundary edge of the channel. The fluid flow issimulated with OpenFOAM’s incompressible solver: PimpleDyMFoam. Thetotal drag and lift forces along the cylinder, as well as along the rigid beam,are recorded. The drag and lift over the time span t = [9 ,
10] seconds can befound in Figure 7 and Table 2. We see that the drag and lift forces computedusing PimpleDyMFoam are converging to Turek & Hron’s values with increasingnumber of control volumes. At the finest level of refinement, with 606,208control volumes, the percentage differences for the peak drag and lift values arewell below 2 %.
Figure 7: Drag and lift forces along the cylinder and the cantilever at time t = [9 ,
10] seconds: δt = 0 .
005 sec
Number of control volumes drag force [centerline ± amplitude] frequency [drag] lift force [centerline ± amplitude] frequency [lift]37,888 470.29 ± ± ± ± ± ± ± ± Table 2: Peak drag and lift forces along the cylinder and the cantilever at time t = [9 , δt = 0 .
005 sec .3. FSI3 Test The complete FSI modeling context is treated in FSI3 test case, where an elasticstructure (Young’s modulus of 5.6 MPa, Poisson ratio of 0.4, and mass densityof 1000 kg/m ) is submerged in a fluid channel that is 2.5 m long, 0.41 mtall, and 0.01 m wide; a summary of the material properties for the fluid andstructural domains is provided in Table 3. The channel inlet, on the left, is pre-scribed with a parabolic inflow speed with mean inflow velocity of 2 m/s. Thefluid flow is modeled with OpenFOAM’s PimpleDyMFoam solver, while, thestructural deformation is estimated using the Generalized- α integration methodwithin CU-BENs (a spectral radius of 0.8 is applied to ensure numerical sta-bility.) Deformation along the fluid and structural interface is governed by theIQN-ILS scheme. Displacements along the moving boundary are iterated untilthe residual vector ( R k ) satisfies the convergence criterion (cid:13)(cid:13) R k (cid:13)(cid:13) ≤ (cid:15) , where theconvergence tolerance is (cid:15) = 10 − . Table 3: Engineering properties
Fluid Properties Structural PropertiesDensity, ρ f kg/m Density, ρ s kg/m Viscosity, µ f − m /s Poisson ratio, µ s δt = 0 . t = [19 . ,
20] seconds as well as thepeak x– and y– tip displacements are provided in Figure 8 and Table 4. Theanalysis is repeated using three levels of mesh refinement in the fluid domain,meanwhile the discretization in the structural domain is fixed at 80 DKT shellelements. As shown, the implicit FSI solver is able to predict the structural de-formation reasonably well compared to the reference values provided by Turek& Hron. The peak displacements are slightly higher than the benchmark values21ith the peak x– displacement converging towards − . ± .
75 mm and thepeak y– displacement converging towards 1 . ± .
97 mm. Nonetheless, thediscrepancies among the peak displacements are within 1 mm in x– directionand 2 mm in y– direction, approximately 12 .
5% difference for x– displacementand 4 .
1% difference for y– displacement.
80 CSD shells & 9,472 CFD CVs80 CSD shells & 26,160 CFD CVs80 CSD shells & 37,888 CFD CVs
80 CSD shells & 9,472 CFD CVs80 CSD shells & 26,160 CFD CVs80 CSD shells & 37,888 CFD CVs
Figure 8: x– and y– displacements at the free end of the cantilever during t = [19 . , δt = 0 . Number of shells Number of control volumes x– displacement [centerline ± amplitude] frequency [x– disp] y– displacement [centerline ± amplitude] frequency [y– disp]80 9,472 -2.46 ± ± ± ± ± ± ± ± Table 4: Peak x– and y– displacements at the free end of the cantilever during time t =[19 . ,
20] seconds: δt = 0 . The drag and lift responses over the time span t = [19 . ,
20] seconds are pre-sented in Figure 9 and Table 5. Minor nonphysical oscillations can be observeddue to the different spatial discretization in the fluid and structural systems aswell as nonconforming mesh along the fluid-structure interface. The implicitFSI solver estimates the peak drag converge towards 489 . ± .
44 N, and thepeak lift converge towards 11 . ± .
21 N. The discrepancies among the peakdrag and lift are within 40 N and 170 N respectively, with percentage differenceof 8% for drag and 26 .
5% for lift. The discrepancies are due to the difference influid solver used in the present work ( i.e. different spatial temporal discretiza-tion scheme, different approach to solving the pressure-velocity equations, etc. )22s well as the different FSI approach compared to Turek & Hron. As specificallynoted in [52], with different discretization, solver, and coupling mechanism, dif-ferences in different approach could lead to discrepancies of up to 50% for thedrag and lift values and 10% difference for the displacement values. With thisin mind, we are confident in our FSI predictions, as we are well below thesethresholds.
80 CSD shells & 9,472 CFD CVs80 CSD shells & 26,160 CFD CVs80 CSD shells & 37,888 CFD CVs
80 CSD shells & 9,472 CFD CVs80 CSD shells & 26,160 CFD CVs80 CSD shells & 37,888 CFD CVs
Figure 9: Drag and lift forces along the the cantilever during time t = [19 . ,
20] seconds: δt = 0 . Number of shells Number of control volumes drag force [centerline ± amplitude] frequency [drag] lift force [centerline ± amplitude] frequency [lift]80 9,472 488.9 ± ± ± ± ± ± ± ± Table 5: Peak drag and lift forces during time t = [19 . ,
20] seconds: δt = 0 .
4. Bridging Simulations
In the previous section, we identified canonical verification models to furnishground truth results for assessing the performance of our proposed FSI simu-lation framework. These versification analyses allowed us to identify suitablespatiotemporal discretization for our application space, along with gauging theefficacy of our selected numerical schemes. As it is that we need to extend ouranalyses to be more realistic with respect to our ultimate design application(fiber reinforced polymer UUAV sail plane), we now incrementally increase our23odeling complexity using what refer to as bridging simulations . We performBayesian optimization (both constrained and unconstrained) using our bridgingsimulations. The geometric configuration and boundary conditions of the FSI3benchmark, outlined in Section 3, are leveraged in the bridging simulations, asdetailed in Section 4.1 to 4.3. As it is that our ultimate design application in-volves the material property design for a composite fiber reinforced structuralcomponent of fixed geometry, in our bridging simulations the elastic modulusof our structure is specified as a function of its cantilevered geometry withinour bridging simulations: as linear, uniform, and box functions, respectively.A schematic representation of the BO-FSI procedure is presented in Figure 10,where x is the optimization parameter (also referred to as input), f ( x ) is theobjective function, c k ( x ) is the set of constraint functions, and cEI ( x ) is theacquisition function. The BO process goes as follows:1. We initialize the optimizer by defining a set of initialization inputs (chosenrandomly except in section 4.3 and 5).2. We pass the initialization inputs to the FSI solver, which outputs corre-sponding data on the objective function and the constraints.3. The BO solver computes the GP surrogates, maximizes the acquisitionfunction, and provide a new input point.4. The new input point is passed to the FSI solver, and the process it repeateduntil the desired design is uncovered. Figure 10: A schematic representation of the interplay between the FSI and Bayesian opti-mization solvers a priori stopping crite-rion; though there may be a hard limit on the number of feasible iterations ( e.g. available computational budget for evaluating the objective function): knowingwhen to stop may be difficult. There is a trade-off to balance between objec-tive function evaluation cost and the likelihood of finding a better (the global)minimum. The optimizer is typically stopped when no significant progress hasbeen made for the last few iterations, or if the uncovered design is simply “goodenough” for engineering purposes.
In this first example, we aim to minimize the vertical tip displacement of thecantilever beam, used in the FSI3 verification, subject to an outlet velocity con-straint. In this section, and those subsequent, the tip displacement correspondsto the maximum displacement attained at any time point within the stabilizedFSI run. The beam elasticity modulus, E , is assumed to vary linearly along thelength of the cantilever (with zero coordinate corresponding to the fixed end, atthe cylinder): E ( x ) = Ax + B x ∈ [0 , l ] (cid:90) l E ( x ) dx = 1 .
96 MPa (31)The integral of the elastic modulus function, taken along the beam, remainsconstant for all time; this results in B being dependant on A . In addition, toensure E ( x ) > A is bounded in the interval [ − . , . A , such that the vertical beam tip displacement, δ , is minimized, while also making sure that we do not exceed our prescribedmaximum outlet fluid velocity, v c = 1 . A δ s.t. v x ≤ v c (32)The Bayesian optimizer is initialized using four different FSI runs; instantiatedwith random values for A. Table 6 shows the successive input values proposed bythe optimizer, as well as the corresponding tip displacement and outlet velocity.Figure 11 shows plots of the surrogate GPs for δ ( A ) and v x ( A ), along withthe constrained expected improvement at the initialization, and after each BOiteration. Ultimately, the tip displacement, and the outlet velocity, both appearto have a monotonic trend; thus the optimizer pushes us towards the lowerbound of A ( i.e. where the stiffness is concentrated at the beam’s fixed end).Once ( A min = − . , δ min = 0 . A (MPa) δ (m) v x (m/s)Init − .
565 0 . . − .
957 0 .
035 1 . .
033 0 .
048 1 . .
764 0 . . − . . . − . . . − . . . − . . . − . . . Table 6: Successive BO iterations - tip displacement with outlet velocity constraint T i p D i s p l a c e m e n t Predictive Mean * Predictive Variance *
30 20 10 0 10 20 301.21.41.6 O u t l e t V e l o c i t y Unfeasible Regions
30 20 10 0 10 20 30A0.00000.00020.0004 A c q u i s i t i o n F un c t i o n x next (a) Initialization
30 20 10 0 10 20 300.030.040.05 T i p D i s p l a c e m e n t
30 20 10 0 10 20 301.21.41.6 O u t l e t V e l o c i t y
30 20 10 0 10 20 30A0.00000.00050.0010 A c q u i s i t i o n F un c t i o n (b) Iteration 1
30 20 10 0 10 20 300.030.040.05 T i p D i s p l a c e m e n t
30 20 10 0 10 20 301.01.21.41.6 O u t l e t V e l o c i t y
30 20 10 0 10 20 30A0.000000.000250.000500.00075 A c q u i s i t i o n F un c t i o n (c) Iteration 2
30 20 10 0 10 20 300.030.040.05 T i p D i s p l a c e m e n t
30 20 10 0 10 20 301.21.41.6 O u t l e t V e l o c i t y
30 20 10 0 10 20 30A0.00000.00020.00040.0006 A c q u i s i t i o n F un c t i o n (d) Iteration 3
30 20 10 0 10 20 300.030.040.05 T i p D i s p l a c e m e n t
30 20 10 0 10 20 301.01.21.41.6 O u t l e t V e l o c i t y
30 20 10 0 10 20 30A0246 A c q u i s i t i o n F un c t i o n
1e 8 (e) Iteration 4
30 20 10 0 10 20 300.030.040.05 T i p D i s p l a c e m e n t
30 20 10 0 10 20 301.01.5 O u t l e t V e l o c i t y
30 20 10 0 10 20 30A0.00000.00020.0004 A c q u i s i t i o n F un c t i o n (f) Iteration 5Figure 11: Successive BO iterations - tip displacement with outlet velocity constraint. Thesolid blue and orange lines represent the predictive mean of the two surrogate GPs ( δ and v x respectively) with their 95% confidence interval (shaded areas), given the already known tipdisplacements and outlet velocity values. The acquisition function indicates what value of A is the most likely to yield a better global minimum ( i.e. the value that should be use in thenext FSI run). T i p D i s p l a c e m e n t Figure 12: Best minimum recorded at each BO iteration - tip displacement with outlet velocityconstraint .2. Example 2: Unconstrained Uniform Stiffness In this second example, we now aim to minimize the vertical tip displacement ofthe cantilever with respect to B in Equation 31, assuming that the beam has anuniform stiffness ( A = 0) ( i.e. integral condition on E ( x ) is no longer enforced). B is bounded in the interval [0 , B δ (33)This optimization problem may seem overly simple, because one might intu-itively expect that larger stiffness values always yield a smaller tip displacement.However, the initialization data shows the opposite, as B = 2 leads to a smaller δ than B > . B = 5 . B = 2, are compared. As shown in Figure 13, underthe same conditions, the cantilever with lower stiffness exhibits a higher modeof oscillation, under the action of the complex vortex shedding within the vonK´arm´an vortex street . As a result of this more complex motion, the less stiffbeam exhibits a smaller tip deflection than the stiffer case. This highlights someinteresting, and counter intuitive, results that may arise within FSI contexts.Figure 14a shows the GP surrogate after 11 BO iterations and indicates that δ ( B ) follows a highly non-monotonic and non-convex behavior, with a significantdrop in tip displacement for B ≈
2. This drop is sharp, which makes the searchfor the global minimum difficult, but after 6 iterations, ( B min = 1 . , δ min =0 . a) simulation time = 6 .
25 sec;B= 5.6 (b) simulation time = 6 .
25 sec;B= 2(c) simulation time = 12 . . B = 5 . B = 2 at time t= 6 .
25 sec and 12 . T i p D i s p l a c e m e n t A c q u i s i t i o n F un c t i o n (a) Tip displacement with uniformstiffness (11 th iteration) T i p D i s p l a c e m e n t (b) Best minimum recorded at eachBO iterationFigure 14: Successive BO iterations - unconstrained uniform stiffness B (MPa) δ (m)Init 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 7: Successive BO iterations - unconstrained uniform stiffness
In this third example, we consider the stiffness function along the beam to bedescribed by a box-function as illustrated in Figure 15. Formally, the stiffnessis defined by: E ( x ) = . x ∈ [ x b − l/ , x b + l/ . x b . Figure 16a shows the surrogate GP after 11iterations. The tip displacement tends to decrease, overall, when the zone ofconcentrated material stiffness ( i.e the box ) is positioned close to the beam’sfixed end. Interestingly, the minimum tip deflection value, or δ min = 0 . xll /3 l /6 l /6 x b Figure 15: Stiffness along the beam given the box mid-point, x b not obtained for a box fully shifted to the far left ( i.e. corresponding to the fixedend), but rather for a box located at about a third of the beam length ( x b, min =0 . x b value proposed by the optimizer is x b = 0 . T i p D i s p l a c e m e n t x b A c q u i s i t i o n F un c t i o n (a) Tip displacement with boxfunction stiffness (11 th iteration) T i p D i s p l a c e m e n t (b) Best minimum recorded at eachBO iterationFigure 16: Successive BO iterations - non-uniform stiffness x b (m) δ (m)Init 0 . . . . . . . . . . . . . . . . . . . . . . . .
10 0 . .
11 0 . . Table 8: Successive BO iterations - non-uniform stiffness
5. Composite Material Properties Selection for UUAV Sail Plane
Fortified with our experience with verification and bridging simulations, we nowturn our attention towards a realistic application: a UUAV sail plane made fromFRP (fiber reinforced polymer), as shown in Figure 17. The sail plane geometryis fixed: 0 . .
01 m wide, with some adjustable angle of attack, θ (with respect to the horizontal). The leading edge of the sail plane is 0 .
038 mthick; reducing linearly with length, to assume a value of 0 .
004 m thickness atthe trailing edge. The sail plane has a density of 2000 kg/m . A uniform inflowvelocity of 5 m/s is prescribed along the left boundary. We wish to find theoptimal (uniform) stiffness, E , and angle of attack θ , that minimizes the dragover time, while imposing conditions on the lift, the vertical tip displacement,as well as the outlet pressure differential. Unlike previous cases, this is a 2-32arameters optimization problem: Figure 17: UUAV geometry configuration
Find min θ,E F D s.t. F L ≥ F L,c δ ≤ δ c ∆ p ≤ ∆ p c (35)The drag F D and lift F L values are respectively taken as the sum of the totalpressure acting on the sail plane in the tangential and normal direction, withrespect to the orientation of the sail plane, over the time span [0 . , . δ is the plane vertical tip displacement. The pressure differential, ∆ p , is takenas ∆ p = p max − p min , where p max and p min are the pressures recorded at theoutlet, at a height of 5 / H , during the time interval [0 . , . F L,c = 12000 N, δ c = 2 . p c = 120 kPa. The BayesianOptimizer is initialized with eight data points (since the input space is now intwo dimensions, more data are required for fitting the surrogate GPs properly).Table 9 shows the successive input parameters proposed by the optimizer, alongwith the corresponding FSI analysis results. Figure 19 shows the predictivemean of the drag force predicted by the surrogate GP at, each BO iteration, withthe feasible/unfeasible inputs regions (given the constraint set for the surrogateGPs). Figure 18 shows the best minimum found after each iteration. Thedrag appears to follow a very smooth and fairly monotonic behavior, and theoptimizer converges quickly toward a minimum. As seen in Figure 19, thefeasible regions provide significant freedom for E to vary, whereas the feasibleregions for θ are much more narrow. This is intuitive since one might expectthe angle of the sail plane to influence the drag more than the stiffness. Yet,33he optimal stiffness value is not trivial. Given the optimizer progress after 8iterations (Figure 18), it seems unlikely that a significantly better minimumdrag value can be achieved, and the iteration loop is stop there. The globalminimum is taken as ( θ min = 3 . , E min = 35 . , F D,min = 2303).Iteration θ ( ◦ ) E (MPa) F D (N) F L (N) δ (mm) ∆ p (kPa)Init 0 .
00 40 .
00 350 3540 0 .
565 9 . .
50 40 .
00 1482 9214 1 .
206 39 . .
00 40 .
00 3419 13736 2 .
146 89 . .
50 40 .
00 5042 16183 2 .
593 141 . .
00 40 .
00 7148 18333 3 .
399 165 . .
50 30 .
00 1452 10051 1 .
813 42 . .
50 50 .
00 4767 16336 1 .
884 139 . .
00 30 .
00 7026 17716 4 .
840 173 .
801 3 .
85 39 .
93 2454 12025 1 .
739 74 .
532 3 .
66 30 .
00 2256 11891 2 .
355 61 .
833 3 .
75 30 .
00 2247 11320 2 .
183 71 .
424 4 .
02 35 .
77 2466 12467 1 .
749 89 . .
74 35 .
39 2303 12184 1 .
956 50 . .
63 35 .
12 2172 11961 1 .
939 77 .
977 3 .
65 35 .
10 2285 11858 1 .
941 76 .
858 3 .
69 35 .
11 2164 11547 1 .
934 78 . Table 9: Successive BO iterations - UUAV sail plane D r a g Figure 18: Best minimum recorded at each BO iteration - UUAV sail plane E (a) Initialization E (b) Iteration 1 E (c) Iteration 2 E (d) Iteration 3 E (e) Iteration 4 E (f) Iteration 5 E (g) Iteration 6 E (h) Iteration 7 E (i) Iteration 8Figure 19: Successive BO iterations - UUAV sail plane. The color surface plots represent theGP predictive mean for the drag, and the shaded dark areas represent the unfeasible inputregions imposed on the surrogate GPs for the lift, tip displacement, and pressure. . Conclusion This paper has demonstrated a principled approach to design using high fidelityFSI simulations. Our proposed approach consists of three stages. Firstly, verifyand/or validate the FSI computational tool using canonical cases and/or exper-iment results from the literature with design configurations that are close to theultimate application. This helps identify the optimal spatiotemporal discretiza-tion, as well as numerical schemes, for the design application. Subsequently,incrementally and systematically increase the complexity of the verification caseto bring it closer to the design context, using what we term bridge simulations .This step provides information to gauge whether adjustments need to be madeto further improve the computational model. Finally, apply the verified FSIframework, thus now trusted, to the design problem of interest.Our experience with BO has shown it to be able to handle complicated de-sign scenarios that have multiple design parameters, and multiple constraints.Furthermore, it is able to effectively optimize design cases that exhibit bothsmooth and non-smooth objective functions, even with counter-intuitive opti-mal parameters directions.
7. Acknowledgement
The support of the Office of Naval Research (ONR), under the grant N00014-19-1-2034 and contract N6833518C0217, is gratefully acknowledged.
ReferencesReferences [1] Sarah C. Vigmostad, Holavanahalli S. Udaykumar, Jia Lu, and Krishnan B.Chandran. Fluid–structure interaction methods in biological flows with spe-cial emphasis on heart valve dynamics.
International Journal for NumericalMethods in Biomedical Engineering , 26(34):435–470, 2010.362] Ramji Kamakoti and Wei Shyy. Fluid–structure interaction for aeroelasticapplications.
Progress in Aerospace Sciences , 40(8):535 – 558, 2004.[3] Zhengyu Huang, Philip Avery, Charbel Farhat, Jason Rabinovitch, ArmenDerkevorkian, and Lee D. Peterson.
Simulation of Parachute Inflation Dy-namics Using an Eulerian Computational Framework for Fluid-StructureInterfaces Evolving in High-Speed Turbulent Flows .[4] J T Tuitman and ˇS Malenica. Fully coupled seakeeping, slamming, andwhipping calculations.
Proceedings of the Institution of Mechanical En-gineers, Part M: Journal of Engineering for the Maritime Environment ,223(3):439–456, 2009.[5] Hermann G. Matthies and Jan Steindorf. Partitioned strong coupling algo-rithms for fluid–structure interaction.
Computers and Structures , 81(8):805– 812, 2003. K.J Bathe 60th Anniversary Issue.[6] Hermann G. Matthies, Rainer Niekamp, and Jan Steindorf. Algorithms forstrong coupling procedures.
Computer Methods in Applied Mechanics andEngineering , 195(17):2028 – 2049, 2006. Fluid-Structure Interaction.[7] Joris Degroote, Robby Haelterman, Sebastiaan Annerel, Peter Brugge-man, and Jan Vierendeels. Performance of partitioned procedures in fluid–structure interaction.
Computers and Structures , 88(7):446 – 457, 2010.[8] P. Causin, J.F. Gerbeau, and F. Nobile. Added-mass effect in the designof partitioned algorithms for fluid–structure problems.
Computer Methodsin Applied Mechanics and Engineering , 194(42):4506 – 4527, 2005.[9] Christiane F¨orster, Wolfgang A. Wall, and Ekkehard Ramm. The artifi-cial added mass effect in sequential staggered fluid–structure interactionalgorithms. TU Delft, The Netherlands, 2006. European Conference onComputational Fluid Dynamics.[10] Christiane F¨orster, Wolfgang A. Wall, and Ekkehard Ramm. Artificialadded mass instabilities in sequential staggered coupling of nonlinear struc-37ures and incompressible viscous flows.
Computer Methods in Applied Me-chanics and Engineering , 196(7):1278 – 1293, 2007.[11] Yichi Zhang, Daniel Apley, and Wei Chen. Bayesian optimization for ma-terials design with mixed quantitative and qualitative variables, 2019.[12] Peter I. Frazier and Jialei Wang. Bayesian optimization for materials de-sign.
Springer Series in Materials Science , pages 45–75, Dec 2015.[13] Remi Lam, Matthias Poloczek, Peter Frazier, and Karen Willcox. Advancesin bayesian optimization with applications in aerospace engineering. 012018.[14] Doniyor Ulmasov, Caroline Baroukh, Benoit Chachuat, Marc Peter Deisen-roth, and Ruth Misener. Bayesian optimization with dimension scheduling:Application to biological systems, 2015.[15] E. O. Pyzer-Knapp. Bayesian optimization for accelerated drug discovery.
IBM Journal of Research and Development , 62(6):2:1–2:7, 2018.[16] V. Nguyen. Bayesian optimization for accelerating hyper-parameter tuning.In , pages 302–305, 2019.[17] Jacob R. Gardner, Matt J. Kusner, Zhixiang Xu, Kilian Q. Weinberger, andJohn P. Cunningham. Bayesian optimization with inequality constraints.In
Proceedings of the 31st International Conference on International Con-ference on Machine Learning - Volume 32 , ICML’14, pages II–937–II–945.JMLR.org, 2014.[18] Stefan Turek and Jaroslav Hron. Proposal for numerical benchmarkingof fluid-structure interaction between an elastic object and laminar in-compressible flow. In Hans-Joachim Bungartz and Michael Sch¨afer, edi-tors,
Fluid-Structure Interaction , pages 371–385, Berlin, Heidelberg, 2006.Springer Berlin Heidelberg. 3819] Wensi Wu, Justyna Kosianka, Heather Reed, Christopher Stull, andChristopher Earls. Cu-bens: A structural modeling finite element library.
SoftwareX , 11:100485, 2020.[20] Zeljko Tukovic Hrvoje Jasak, Aleksandar Jemcov, editor.
OpenFOAM:A C++ Library for Complex Physics Simulations , Dubrovnik, Croatia,September 2007. IUC, International Workshop on Coupled Methods in Nu-merical Dynamics.[21] H.G. Weller and G. Tabor. A tensorial approach to computational contin-uum mechanics using object-oriented techniques.
Computers in Physics ,12(6):620, 1998.[22] Hrvoje Jasak.
Error analysis and estimation for the finite volume methodwith applications to fluid flows . PhD thesis, Imperial College of Science,Technology and Medicine, 1996.[23] Bernhard Gschaider, Hakan Nilsson, Henrik Rusche, Hrvoje Jasak, MartinBeaudoin, and Vanja Skuric. OpenFOAM-1.6-ext.[24] Ronald Miller, Sung-Eun Kim, and Sung Won Lee. A computational frame-work for fluid–structure interaction on flexible propellers involving largedeformation. Pasadena, California, September 2010. 28th Symposium onNaval Hydrodynamics.[25] J. Chung and G.M. Hulbert. A time integration algorithm for structuraldynamics with improved numerical dissipation: the generalized- α method. Journal of applied mechanics , 60(2):371–375, 1993.[26] Hrvoje Jasak.
Dynamic Mesh Handling in OpenFOAM .[27] Christophe Kassiotis. Which strategy to move the mesh in the computa-tional fluid dynamic code openfoam. 2008.[28] JOHN BATINA.
Unsteady Euler airfoil solutions using unstructured dy-namic meshes . 3929] Frederic J. Blom. Considerations of the spring analogy. volume 32, pages647–668. International Journal for Numerical Methods in Fluids, 2000.[30] C. Farhat, C. Degand, B. Koobus, and M. Lesoinne. Torsional springs fortwo-dimensional dynamic unstructured fluid meshes.
Computer Methods inApplied Mechanics and Engineering , 163(1):231 – 245, 1998.[31] Y. Zhao, J. Tai, and F. Ahmed. Simulation of micro flows with mov-ing boundaries using high-order upwind fv method on unstructured grids.
Computational Mechanics , 28(1):66–75, 2002.[32] Blair Perot and Ramesh Nallapati. A moving unstructured staggered meshmethod for the simulation of incompressible free-surface flows.
Journal ofComputational Physics , 184(1):192 – 214, 2003.[33] Christoph Degand and Charbel Farhat. A three-dimensional torsionalspring analogy method for unstructured dynamic meshes.
Computers andStructures , 80(3):305 – 316, 2002.[34] Hrvoje Jasak and Zeljko Tukovi´c. Automatic mesh motion for the unstruc-tured finite volume method, 2004.[35] C. Farhat, M. Lesoinne, and P. Le Tallec. Load and motion transfer algo-rithms for fluid/structure interaction problems with non-matching discreteinterfaces: Momentum and energy conservation, optimal discretization andapplication to aeroelasticity.
Computer Methods in Applied Mechanics andEngineering , 157(1):95 – 114, 1998.[36] Donald Shepard. A two-dimensional interpolation function for irregularly-spaced data. In
Proceedings of the 1968 23rd ACM National Conference ,ACM ’68, pages 517–524, New York, NY, USA, 1968. Association for Com-puting Machinery.[37] Joris Degroote, Klaus-J¨urgen Bathe, and Jan Vierendeels. Performance of anew partitioned procedure versus a monolithic procedure in fluid–structure40nteraction.
Computers & Structures , 87(11):793 – 801, 2009. Fifth MITConference on Computational Fluid and Solid Mechanics.[38] Gene H Golub and Charles F. Van Loan.
Matrix Computations . JohnsHopkins University Press, 4th edition, 2013.[39] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesianoptimization of machine learning algorithms. In F. Pereira, C. J. C. Burges,L. Bottou, and K. Q. Weinberger, editors,
Advances in Neural InformationProcessing Systems 25 , pages 2951–2959. Curran Associates, Inc., 2012.[40] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas. Takingthe human out of the loop: A review of bayesian optimization.
Proceedingsof the IEEE , 104(1):148–175, 2016.[41] Peter I. Frazier. A tutorial on bayesian optimization, 2018.[42] Carl Edward Rasmussen and Christopher K.I. Williams.
Gaussian Pro-cesses for Machine Learning . MIT Press, Cambridge, Mass. [u.a.], 2006.[43] Christophe Bonneville, Maxwell Jenquin, Juan Lundono, Alex Kelly, Jef-frey Cipolla, and Christopher Earls. Gaussian processes for shock testemulation.
In Review , 2020.[44] Christopher M. Bishop.
Pattern Recognition and Machine Learning (Infor-mation Science and Statistics) . Springer-Verlag, Berlin, Heidelberg, 2006.[45] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.
The Elements ofStatistical Learning . Springer Series in Statistics. Springer New York Inc.,New York, NY, USA, 2001.[46] Kevin P. Murphy.
Machine learning : a probabilistic perspective . MITPress, Cambridge, Mass. [u.a.], 2013.[47] Radford M. Neal.
Bayesian Learning for Neural Networks . PhD thesis,CAN, 1995. AAINN02676. 4148] Jungtaek Kim and Seungjin Choi. On local optimizers of acquisition func-tions in bayesian optimization, 2019.[49] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton,Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoTorch:Programmable Bayesian Optimization in PyTorch. arXiv e-prints , pagearXiv:1910.06403, October 2019.[50] Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, andAndrew Gordon Wilson. Gpytorch: Blackbox matrix-matrix gaussian pro-cess inference with gpu acceleration. In
Advances in Neural InformationProcessing Systems , 2018.[51] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, EdwardYang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, andAdam Lerer. Automatic differentiation in pytorch. 2017.[52] S. Turek, J. Hron, M. Razzaq, H. Wobker, and M. Sch¨afer. Numericalbenchmarking of fluid-structure interaction: A comparison of different dis-cretization and solution approaches. In Hans-Joachim Bungartz, MiriamMehl, and Michael Sch¨afer, editors,