Adjoint-Based Optimization of Time-Dependent Fluid-Structure Systems using a High-Order Discontinuous Galerkin Discretization Matthew J. Zahr † and Per-Olof Persson ECCOMAS VII International Conference on Coupled Problems in Science and Engineering Rhodes Island, Greece, June 12-14, 2017 † Luis W. Alvarez Postdoctoral Fellow Department of Mathematics Lawrence Berkeley National Laboratory University of California, Berkeley 1 / 23
Optimization of unsteady fluid-structure systems • Discover energetically optimal flapping motions • understand biological systems, design Micro Aerial Vehicles (MAVs) • Design minimal weight structures and vehicles that will not flutter • Design energy harvesting mechanisms θ ( µ , t ) h ( u s ) c 2 / 23
Unsteady single physics optimization formulation Goal : Find the solution of the unsteady PDE-constrained optimization problem minimize J ( U , µ ) U , µ subject to C ( U , µ ) ≤ 0 ∂ U ∂t + ∇ · F ( U , ∇ U ) = 0 in v ( µ , t ) where • U ( x , t ) PDE solution • µ design/control parameters � T f � • J ( U , µ ) = j ( U , µ , t ) dS dt objective function T 0 Γ � T f � • C ( U , µ ) = c ( U , µ , t ) dS dt constraints T 0 Γ 3 / 23
High-order discretization of PDE-constrained optimization • Continuous PDE-constrained optimization problem minimize J ( U , µ ) U , µ subject to C ( U , µ ) ≤ 0 ∂ U ∂t + ∇ · F ( U , ∇ U ) = 0 in v ( µ , t ) • Fully discrete PDE-constrained optimization problem minimize J ( u 0 , . . . , u N t , k 1 , 1 , . . . , k N t ,s , µ ) u 0 , ..., u Nt ∈ R N u , k 1 , 1 , ..., k Nt,s ∈ R N u , µ ∈ R n µ subject to C ( u 0 , . . . , u N t , k 1 , 1 , . . . , k N t ,s , µ ) ≤ 0 u 0 − g ( µ ) = 0 s � u n − u n − 1 − b i k n,i = 0 i =1 Mk n,i − ∆ t n r ( u n,i , µ , t n,i ) = 0 4 / 23
Highlights of globally high-order discretization n da • Arbitrary Lagrangian-Eulerian formulation: N dA G , g , v X v Map, G ( · , µ , t ) , from physical v ( µ , t ) to reference V x 2 V x 1 X 2 � ∂ U X � + ∇ X · F X ( U X , ∇ X U X ) = 0 � X 1 ∂t � X Mapping-Based ALE • Space discretization : discontinuous Galerkin M ∂ u ∂t = r ( u , µ , t ) • Time discretization : diagonally implicit RK s � DG Discretization u n = u n − 1 + b i k n,i i =1 c 1 a 11 Mk n,i = ∆ t n r ( u n,i , µ , t n,i ) c 2 a 21 a 22 . . . ... . . . . . . • Quantity of interest : solver-consistency c s a s 1 a s 2 · · · a ss b 1 b 2 · · · b s F ( u 0 , . . . , u N t , k 1 , 1 , . . . , k N t ,s ) Butcher Tableau for DIRK 5 / 23
Adjoint method to efficiently compute gradients of QoI • Consider the fully discrete output functional F ( u n , k n,i , µ ) • Represents either the objective function or a constraint • The total derivative with respect to the parameters µ , required in the context of gradient-based optimization, takes the form N t N t s d F d µ = ∂F ∂F ∂ u n ∂F ∂ k n,i � � � ∂ µ + ∂ µ + ∂ u n ∂ k n,i ∂ µ n =0 n =1 i =1 • The sensitivities, ∂ u n ∂ µ and ∂ k n,i ∂ µ , are expensive to compute, requiring the solution of n µ linear evolution equations • Adjoint method : alternative method for computing d F d µ that require one linear evolution equation for each quantity of interest, F 6 / 23
Adjoint equation derivation: outline • Define auxiliary PDE-constrained optimization problem minimize F ( u 0 , . . . , u N t , k 1 , 1 , . . . , k N t ,s , µ ) u 0 , ..., u Nt ∈ R N u , k 1 , 1 , ..., k Nt,s ∈ R N u subject to R 0 = u 0 − g ( µ ) = 0 s � R n = u n − u n − 1 − b i k n,i = 0 i =1 R n,i = Mk n,i − ∆ t n r ( u n,i , µ , t n,i ) = 0 • Define Lagrangian N t N t s � � � T R 0 − T R n − T R n,i L ( u n , k n,i , λ n , κ n,i ) = F − λ 0 λ n κ n,i n =1 n =1 i =1 • The solution of the optimization problem is given by the Karush-Kuhn-Tucker (KKT) sytem ∂ L ∂ L ∂ L ∂ L = 0 , = 0 , = 0 , = 0 ∂ u n ∂ k n,i ∂ λ n ∂ κ n,i 7 / 23
Dissection of fully discrete adjoint equations • Linear evolution equations solved backward in time • Primal state/stage, u n,i required at each state/stage of dual problem • Heavily dependent on chosen ouput T ∂F λ N t = ∂ u N t s T ∂F ∂ r � ∂ u ( u n,i , µ , t n − 1 + c i ∆ t n ) T κ n,i λ n − 1 = λ n + + ∆ t n ∂ u n − 1 i =1 s T ∂F ∂ r � ∂ u ( u n,j , µ , t n − 1 + c j ∆ t n ) T κ n,j M T κ n,i = + b i λ n + a ji ∆ t n ∂ u N t j = i • Gradient reconstruction via dual variables N t s d F d µ = ∂F T ∂ g T ∂ r � � ∂ µ + λ 0 ∂ µ ( µ ) + ∆ t n κ n,i ∂ µ ( u n,i , µ , t n,i ) n =1 i =1 [Zahr and Persson, 2016] 8 / 23
Energetically optimal flapping vs. required thrust Energy = 0.21935 Energy = 3.00404 Energy = 6.2869 Thrust = 0.0000 Thrust = 1.5000 Thrust = 2.5000 Optimal Optimal Optimal T x = 0 T x = 1 . 5 T x = 2 . 5 9 / 23
Structure: semi-discretization, first-order form M s ∂ u s ∂t = r s ( u s ; t ) = r ss ( u s ) + r sf · t • Semidiscretization (CG-FEM) of continuum (hyperelasticity) ∂ p in Ω 0 ∂t − ∇ · P ( G ) = b P ( G ) · N = t on Γ N x = x D on Γ D • Force balance on rigid body M ∂ 2 q ∂t 2 + C ∂ q ∂t + Kq = t θ ( µ , t ) h ( u s ) c 10 / 23
Coupled fluid-structure formulation • Write discretized fluid and structure equations as ODEs M f ˙ u f = r f ( u f ; x ) M s ˙ u s = r s ( u s ; t ) = r ss ( u s ) + r sf · t in the fluid u f and structure u s variables • Apply couplings • Structure-to-fluid: deform fluid domain x = x ( u s ) • Fluid-to-structure: apply boundary traction t = t ( u f ) • Write coupled system as M ˙ u = r ( u ) � � � � � � u f r f ( u f ; x ( u s )) M f u = r ( u ) = M = M s u s r s ( u s ; t ( u f )) 11 / 23
High-order partitioned FSI solver: IMEX Runge-Kutta 1 • Exploit linear dependence of structure residual ( r s ) on traction ( t ) � � � � � � r f ( u f ; x ( u s ) r f ( u f ; x ( u s )) r ( u ) = = + r sf · ( t ( u f ) − ˜ r s ( u s ; ˜ r s ( u s ; t ( u f )) t ) t ) � �� � � �� � f ( u ) g ( u ) • Apply high-order implicit-explicit Runge-Kutta scheme to discretize M ∂ u ∂t = r ( u ) = f ( u ) + g ( u ) � �� � ���� explicit implicit c, ˆ A, ˆ • Explicit Runge-Kutta scheme ˆ b for f ( u ) • Diagonally implicit scheme c, A, b for g ( u ) s s � ˆ b i ˆ � u n = u n − 1 + k n,i + b i k n,i i =1 i =1 � i − 1 i � � a ij ˆ � Mk n,i = ∆ t n g u n − 1 + ˆ k n,j + a ij k n,j j =1 j =1 � i − 1 i � M ˆ � a ij ˆ � k n,i = ∆ t n f u n − 1 ˆ k n,j + a ij k n,j j =1 j =1 1 [van Zuijlen and Bijl, 2005, Froehle and Persson, 2014] 12 / 23
High-order partitioned FSI solver: IMEX Runge-Kutta 1 • Exploit linear dependence of structure residual ( r s ) on traction ( t ) � � � � � � r f ( u f ; x ( u s ) r f ( u f ; x ( u s )) r ( u ) = = + r sf · ( t ( u f ) − ˜ r s ( u s ; ˜ r s ( u s ; t ( u f )) t ) t ) � �� � � �� � f ( u ) g ( u ) • Solve: (1) implicit structure , (2) implicit fluid , (3) explicit structure • Due to choice of IMEX partition: no explicit fluid stages 1 [van Zuijlen and Bijl, 2005, Froehle and Persson, 2014] 12 / 23
High-order partitioned FSI solver: IMEX Runge-Kutta • Define stage solutions i i − 1 � � s a ij k s a ij ˆ u s n,i = u s n − 1 + n,j + ˆ k n,j j =1 j =1 i � u f n,i = u f a ij k f n − 1 + n,j j =1 • Define traction predictor as true traction at previous stage ˜ t n,i = t ( u n,i − 1 ) • Solve for stage velocities ( i = 1 , . . . , s ) M s k s n,i = ∆ t n r s ( u s n,i ; ˜ t n,i ) M f k f n,i = ∆ t n r f ( u f n,i ; x ( u s n,i )) s M s ˆ n,i = ∆ t n r sf ( t ( u f n,i ) − ˜ k t n,i ) • Update state solution at new time s s s � � � s n = u f b j k f ˆ b j ˆ u f u s n = u s b j k s n − 1 + n,j , n − 1 + n,j + k n,j j =1 j =1 j =1 13 / 23
Validation: benchmark pitching airfoil system • Simple FSI benchmark problem for studying the high-order accuracy of the IMEX scheme • Rigid pitching/heaving NACA 0012 airfoil, torsional spring • Smooth heaving step y ( t ) prescribed, angle θ ( t ) measured Setup Mach number 14 / 23
Validation: benchmark pitching airfoil system • Up to 5th order of convergence in time. • Similar accuracy as solving fully coupled system 0.06 nofluid 2e-1 1e-1 10 0 4e-2 0.05 2e-2 2e-3 2e-4 0.04 10 -1 theta 0.03 0.02 10 -2 1 0.01 1 Relative error in θ ( t ) 0.00 0.0 0.5 1.0 1.5 2.0 T 10 -3 Angle θ ( t ) vs time t 10 -4 3 1 10 -5 5 Weak Coupling 4 10 -6 ARK3 1 1 FC-ARK3 ARK4 10 -7 FC-ARK4 ARK5 FC-ARK5 10 -8 10 -2 10 -1 10 0 Entropy Time step ∆ t 15 / 23
Validation: cantilever system • Standard FSI benchmark problem. • Elastic cantilever behind a square bluff body in incompressible flow. • Cantilever: ρ s = 100 kg / m 3 , ν s = 0 . 35 , E = 2 . 5 × 10 5 Pa. • Fluid & Flow: 12.0 1.0 0.06 4.0 ρ f = 1 . 18 kg / m 3 , 1.0 ν f = 1 . 54 × 10 − 5 m 2 / s, v f = 0 . 513 m / s, Re = 333 , Ma = 0 . 2 . 5.5 14.0 • Vortex shedding frequency: ∼ 6 . 3 Hz • Cantilever first mode: 3 . 03 Hz 16 / 23
Recommend
More recommend
