Integration of Fluid Dynamics and Solid Mechanics Models for FSI - PDF document

Transactions of the Korean Nuclear Society Virtual Spring Meeting July 9-10, 2020 Integration of Fluid Dynamics and Solid Mechanics Models for FSI Simulation using GPU- based SPH Framework Tae Hoon Lee a , So Hyun Park a , Eung Soo Kim a  a Department of Nuclear Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, South Korea * Corresponding author: kes7741@snu.ac.kr 1. Introduction 2.2 SPH Approximation Mathematically, random functions can be expressed in In many engineering fields, pressure from fluid flow integral form using delta functions. The SPH method is causes major deformation in the structure. This represented in integral form using the kernel function, a interaction is Fluid Structure Interaction (FSI). Among continuous function with a smoothing length h instead of the many ways to describe the FSI are the Arbitrary a delta function. Discretizing this integral form is like (1). Lagrangian Eulerian (ALE) typical. However, ALE has disadvantages in phenomena such as large deformation 𝑔 𝑗 (𝑠) = ∑ 𝑔 𝑘 𝑋(𝑠 𝑗𝑘 , ℎ)𝑊 (1) of the structure or rapid flow of fluid [1]. On the other 𝑘 𝑘 hand, the Smoothed Particle Hydrodynamics (SPH) method, which is a Fully Lagrangian method, has Where i stands for the central particle and j stands for advantage to this interpretation. the surrounding particle. 𝑿(𝒔 𝒋𝒌 , 𝒊) represents the kernel FSI is a method used in many areas of nuclear function, 𝒔 𝒋𝒌 = 𝒔 𝒋 − 𝒔 𝒌 , 𝑾 𝒌 is volume of adjacent engineering, such as the behavior of the fuel assembly particle and h is the smoothing length indicating the flowing axially in the direction of coolant, pressure loads range of nearby particles to be included in the on the vessel internal structures in PWR during a LOCA, approximation process. The kernel function must be able blowdown, Flow-Induced Vibration (FIV), flow-induced to approximate the delta function, so it has a very large fluid-elastic vibrations, sloshing of pressurizer on a value at the center and the farther away from the center, nuclear ship, rupture or swelling of fuel rods and In- the more convergent it is to zero. It also satisfies all the Vessel Retention (IVR) failed due to broken vessel. This mathematical properties of the delta function. is because accurately simulating the interaction of fluid There are several ways to obtain a gradient of function, and structure can help design the plant and cope with depending on the method of deriving. In this study, a accidents. gradient of function was calculated in the following In this study, the interaction of fluid with structure was manner added to the SOPHIA code to implement the FSI. The SOPHIA code is an SPH-based parallelization multi- 𝑔 i + 𝑔 j ∇ 𝑔 𝑗 (𝑠) = 𝜍 𝑗 ∑ m j ( ρ i ρ j ) ∇ 𝑋(𝑠 𝑗𝑘 , ℎ) (2) physics code developed by Seoul National University [2]. j To validation this, Benchmark experiment and numerical simulation were compared with simulation through 3. Fluid Structure Interaction (FSI) SOPHIA code with FSI added. 3.1 Fluid Dynamics 2. SPH Method In the SPH method, there are two methods for solving 2.1 Concept of SPH method mass conservation equations: mass summation method and solving continuity equation method. The mass Smoothed Particle Hydrodynamics was first summation method was used in this study. developed in astrophysics as a meshless CFD method of the Lagrangian-based [3]. The fluid is represented as a 𝜍 𝑗 (𝑠) = ∑ 𝑛 𝑘 𝑋(𝑠 𝑗𝑘 , ℎ) (2) 𝑘 collection of finite particles in a way that tracks and interprets fluid motion, not based on space and grid. The momentum conservation equation is described in Particles that are considered to be a collection of fluid Lagrangian form as (3). molecules move along with physical quantities (mass, velocity, temperature, etc.) that are determined by the ⃗⃗⃗⃗⃗ 𝑒𝑣 𝑒𝑢 = −∇𝑞 + 𝜈∇ 2 𝑣 ρ ⃗ + 𝜍𝑕 (3) type of fluid and the spacing of particles. Track the movement of particles by setting up initial conditions for them and interpreting interactions with the central Where 𝛓 is density, 𝒗 ⃗ is speed, 𝒒 is pressure, 𝒉 ⃗ ⃗ ⃗ is particle and its surrounding particles. The SPH method gravity acceleration and 𝝂 is dynamic viscosity. The first has the advantage of dealing with undetermined areas of term on the right hand side represents the pressure force interpretation or highly variable flows, thanks to the and the second term is the viscous force, which is (4) and nature of the Lagrangian-based analysis method. (5) if they are to be discretized.

Transactions of the Korean Nuclear Society Virtual Spring Meeting July 9-10, 2020 0 ) = 𝑺 𝑗 ∇𝑋 0 ). ⃗⃗⃗⃗⃗ ∗ (𝚿 𝑗𝑘 𝑞 𝑗 +𝑞 𝑘 ∇𝑋 𝑗 (𝚿 𝑗𝑘 𝑒𝑣 ( 𝑒𝑢 ) 𝑔𝑞,𝑗 = − ∑ 𝑛 𝑘 ( 𝜍 𝑗 𝜍 𝑘 ) ∇𝑋 (4) 𝑗 𝑘 𝑗𝑘 ⃗⃗⃗⃗⃗ 4𝑛 𝑘 𝜈 𝑘 𝑠 𝑗𝑘 ⃗⃗⃗⃗⃗ ∙∇𝑋 𝑗𝑘 𝑒𝑣 𝑔𝑤,𝑗 = ∑ Calculate the strain according to the linear elasticity ( 𝑒𝑢 ) ⃗⃗⃗⃗⃗ | 2 +𝜃 2 ) (𝑣 𝑗 ⃗⃗⃗ − 𝑣 𝑘 ⃗⃗⃗ ) (5) 𝑘 (𝜍 𝑗 +𝜍 𝑘 )(|𝑠 𝑗𝑘 model from the 𝐆 ∗ . The pressure is calculated as Equation of State (EOS) ∗ + 𝐆 𝑗 1 ∗𝑈 ) − 𝐉. (10) 𝜁 𝑗 = 2 (𝐆 𝑗 using the Tait's equation [4]. This equation is given as shown in (6) in a manner that assumes weakly Calculate the stress through the Piola-Kirchhoff stress compressible. tensor using strain. G and K stand for Shear and Bulk modulus, respectively. 𝛿 2 𝜍 0 𝑑 0 𝜍 p = 𝛿 [( 𝜍 0 ) − 1] (6) 2 𝑸 𝑗 = 2𝐻𝜁 𝑗 + (𝐿 − 3 𝐻) 𝑢𝑠(𝜁 𝑗 )𝐉. (11) Where 𝝇 𝟏 is the reference density of the fluid, 𝒅 𝟏 is speed of sound and 𝛿 = 7 . Finally, calculate the elastic force through the divergence of the stress sensor [5]. 3.2 Structure Dynamics 0 (𝐐 𝑗 ∇𝑋 0 ) − 𝐐 0 0 𝑊 ∗ (𝚿 𝑗𝑘 ∗ (𝚿 𝑘𝑗 0 )) 𝐠 𝑗 = ∑ 𝑊 𝑘 ∇𝑋 . (12) Elastic force is calculated using the divergence of the 𝑘 𝑗 𝑘 𝑗 𝑘 stress tensor. In order to calculate these stress tensors, we ∇𝑋 𝑘 = −∇𝑋 𝑗 allows the expression to be described as must first calculate the F (Deformation gradient). The change of the solid consists of three types: translation, (4), but in case of (12), the kernel function multiplied by deformation, and rotation. Elastic force depends solely the R of each particle, the expression is described on deformation, so translation and rotation should be separately as (12). extracted from the change of the solid. F (Deformation gradient) indicates how much of the current position has 3.3 Fluid Structure Coupling been changed from the initial position. At this time, F excludes only the effects of translation from the three Fluid and structure are calculated only between the changes in the solid. This is as shown in equation (7) same kinds of particles, so additional forces acting on when SPH approximation is performed. In contrast to A. each other should be calculated. In this study pressure Peer et al [5], F was calculated without using a kernel force was used as the interaction force between Fluid and gradient correction. Structure. Based on the points calculated using pressure force in 0 ) 0 0 𝚿 𝑘𝑗 ⨂∇𝑋 𝐆 i = ∑ 𝑊 𝑗 (𝚿 𝑗𝑘 . (7) the relationship between the Fluid particle and the 𝑘 𝑘 Boundary particle, simulation code was configured to push each other with pressure force even in the Where 𝐘 𝑘𝑗 = 𝐘 𝑘 − 𝐘 𝑗 , 𝐘 𝑗 is position vector and the relationship between the fluid and structure. superscript 0 is initial value. First, use EOS to calculate the pressure from the R (Rotation matrix) must be calculated first to exclude structure particle, as with the boundary particle. Then, for rotation from F. Instead of calculating R directly for each one central structure particle, the pressure force ( 𝑔 𝑄 ) is step, calculate R through iteration using R of the previous calculated if the adjacent particle is fluid particle and the step. Iteration is performed until ω is less th an 1e-10, elastic force ( 𝑔 𝑓 ) if the adjacent particle is structure where ω is the vector where ‖F − 𝑆‖ 2 is the minimum. particle. It is described in Fig.1. R of the next step is calculated by rotating R of the previous step in the same direction as ω in the angle of ω [6]. ∑ 𝑠 𝑗 ×𝑏 𝑗 𝑗 𝐒 ← exp(𝝏)𝐒, 𝐒 ← exp ( |+𝜁 ) 𝐒. (8) |∑ 𝑠 𝑗 ∙𝑏 𝑗 𝑗 𝒈 𝑸 Calculate the 𝐆 ∗ (Corotated deformation gradient), 𝒈 𝒇 which removes the change by rotation from F, to take into account only the effects of the deformation. 𝐆 ∗ is calculated as the difference between the change of the current position and the change of the initial position multiplied by the rotation. Finally, only the effect of the Fig. 1 Fluid Structure Interaction Force deformation will remain in the structure [5]. ∗ = 𝐉 + ∑ 𝑊 0 )⨂∇𝑋 0 ). 4. SPH Simulation on Deformation of Elastic Plate 0 0 (𝚿 𝑘𝑗 − 𝑺 𝑗 𝚿 𝑘𝑗 ∗ (𝚿 𝑗𝑘 𝐆 𝑗 (9) 𝑘 𝑘 𝑗