Variable resolution smoothed particle hydrodynamics schemes for 2-D and 3-D viscous flows
Abstract
Smoothed Particle Hydrodynamics (SPH) is a Lagrangian particle-based method for the numerical solution of the partial differential equations that govern the motion of fluids. The main aim of this thesis work is to better enable the applicability of SPH to problems involving multi-scale fluid dynamics. In the first part of the thesis, the capability of the SPH method to simulate three-dimensional isotropic turbulence is investigated with a detailed comparison of Lagrangian and Eulerian SPH formulations. The main reason for this first investigation is to provide an assessment of the error introduced by the particle disorder on the SPH discrete operators when being purely Lagrangian. When the free decay of isotropic turbulence in a triple periodic box is studied, the Eulerian SPH formulation achieves a very good agreement with other well validated reference solutions, whereas Lagrangian SPH yields an inaccurate prediction of turbulent energy spectra. When considering linearly forced isotropic turbulence, the use of a Godunov-type SPH scheme becomes essential for the achievement of a stable solution. The efficacy of the particle shifting technique applied to turbulent SPH flows is also studied in this part of the thesis and numerical findings indicate that corrective terms derived from the arbitrary Lagrangian–Eulerian theory are essential for a proper estimation of turbulence characteristics. Subsequently, numerical analyses of a decaying isotropic turbulent flow are carried out for the first time using SPH schemes based on high-order kernels. A dramatic increase in the accuracy of the results is observed when high-order SPH is employed, especially for the description of the vorticity dynamics. Motivated by the findings of the computational investigations above, the second part of this thesis focuses on the implementation and testing of a novel SPH variable-resolution algorithm. A domain-decomposition approach is adopted to partition the computational domain into regions having different particle resolutions. Each numerical sub-problem is then closed by appending buffer regions to every sub-domain, and populating these regions with particles whose physical quantities are obtained by means of interpolations over adjacent sub-domains. These interpolations are carried out using a second-order kernel correction procedure to ensure the proper consistency and accuracy of the interpolation process. The mass transfer among sub-domains is modeled by evaluating the Eulerian mass flux at the domain boundaries. Particles that belong to a specific zone are created/destroyed in the buffer regions and do not interact with fluid particles that belong to a different resolution zone. The algorithm is implemented in the DualSPhysics open-source code [62] and optimized thanks to DualSPHysics’ parallel framework. The algorithm is tested on a series of different fluid dynamics problems: a 2-D hydrostatic tank case, a flow past cylinder for different values of the Reynolds number, a flow past an oscillating cylinder in the cross-flow direction, and the propagation of regular waves across a rectangular tank. The present algorithm is able to simulate efficiently fluid dynamics problems characterized by a wide range of spatial scales, achieving a ratio between the coarsest and the finest resolution up to a factor equal to 256.The investigation is then extended to 3-D fluid dynamics problems, such as the flow past a sphere with a Reynolds Number equal to 300 and 500, for which a SPH solution using a uniform resolution is unfeasible due to the high computational cost, showing a good agreement between the results obtained with the variable-resolution algorithm herein presented and relevant numerical investigations in the literature. The work is then concluded with the simulation and validation of a 3-D dam-breaking flow impacting a cubic obstacle.
Full text
ABSTRACT VARIABLE RESOLUTION SMOOTHED PARTICLE HYDRODYNAMICS SCHEMES FOR 2-D AND 3-D VISCOUS FLOWS by Francesco Ricci Smoothed Particle Hydrodynamics (SPH) is a Lagrangian particle-based method for the numerical solution of the partial differential equations that govern the motion of fluids. The main aim of this thesis work is to better enable the applicability of SPH to problems involving multi-scale fluid dynamics. In the first part of the thesis, the capability of the SPH method to simulate three-dimensional isotropic turbulence is investigated with a detailed comparison of Lagrangian and Eulerian SPH formulations. The main reason for this first investigation is to provide an assessment of the error introduced by the particle disorder on the SPH discrete operators when being purely Lagrangian. When the free decay of isotropic turbulence in a triple periodic box is studied, the Eulerian SPH formulation achieves a very good agreement with other well validated reference solutions, whereas Lagrangian SPH yields an inaccurate prediction of turbulent energy spectra. When considering linearly forced isotropic turbulence, the use of a Godunov-type SPH scheme becomes essential for the achievement of a stable solution. The efficacy of the particle shifting technique applied to turbulent SPH flows is also studied in this part of the thesis and numerical findings indicate that corrective terms derived from the arbitrary Lagrangian–Eulerian theory are essential for a proper estimation of turbulence characteristics. Subsequently, numerical analyses of a decaying isotropic turbulent flow are carried out for the first time using SPH schemes based on high-order kernels. A dramatic increase in the accuracy of the results is observed when high-order SPH is employed, especially for the description of the vorticity dynamics.
Motivated by the findings of the computational investigations above, the second part of this thesis focuses on the implementation and testing of a novel SPH variable-resolution algorithm. A domain-decomposition approach is adopted to partition the computational domain into regions having different particle resolutions. Each numerical sub-problem is then closed by appending buffer regions to every sub-domain, and populating these regions with particles whose physical quantities are obtained by means of interpolations over adjacent sub-domains. These interpolations are carried out using a second-order kernel correction procedure to ensure the proper consistency and accuracy of the interpolation process. The mass transfer among sub-domains is modeled by evaluating the Eulerian mass flux at the domain boundaries. Particles that belong to a specific zone are created/destroyed in the buffer regions and do not interact with fluid particles that belong to a different resolution zone. The algorithm is implemented in the DualSPhysics open-source code [62] and optimized thanks to DualSPHysics’ parallel framework. The algorithm is tested on a series of different fluid dynamics problems: a 2-D hydrostatic tank case, a flow past cylinder for different values of the Reynolds number, a flow past an oscillating cylinder in the cross-flow direction, and the propagation of regular waves across a rectangular tank. The present algorithm is able to simulate efficiently fluid dynamics problems characterized by a wide range of spatial scales, achieving a ratio between the coarsest and the finest resolution up to a factor equal to 256. The investigation is then extended to 3-D fluid dynamics problems, such as the flow past a sphere with a Reynolds Number equal to 300 and 500, for which a SPH solution using a uniform resolution is unfeasible due to the high computational cost, showing a good agreement between the results obtained with the variable-resolution algorithm herein presented and relevant numerical investigations in the literature.
The work is then concluded with the simulation and validation of a 3-D dam-breaking flow impacting a cubic obstacle.
VARIABLE RESOLUTION SMOOTHED PARTICLE HYDRODYNAMICS SCHEMES FOR 2-D AND 3-D VISCOUS FLOWS by Francesco Ricci A Dissertation Submitted to the Faculty of New Jersey Institute of Technology in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Mechanical Engineering Department of Mechanical and Industrial Engineering August 2023
Copyright ©2023 by Francesco Ricci ALL RIGHTS RESERVED
APPROVAL PAGE VARIABLE RESOLUTION SMOOTHED PARTICLE HYDRODYNAMICS SCHEMES FOR 2-D AND 3-D VISCOUS FLOWS Francesco Ricci Dr. Angelantonio Tafuni, Dissertation Advisor Date Assistant Professor of Mechanical Engineering, NJIT Dr. Samaneh Farokhirad, Committee Member Date Assistant Professor of Mechanical Engineering, NJIT Dr. Samuel Lieber, Committee Member Date Assistant Professor of Mechanical Engineering Technology, NJIT Dr. Simone Marras, Committee Member Date Assistant Professor of Mechanical Engineering, NJIT Dr. Jos´e Manuel Dom´ınguez Alonso, Committee Member Date Assistant Professor of Applied Physics, Universidade de Vigo, Ourense, Spain
BIOGRAPHICAL SKETCH Author: Francesco Ricci Degree: Doctor of Philosophy Date: August 2023 Undergraduate and Graduate Education: •Doctor of Philosophy in Mechanical Engineering, New Jersey Institute of Technology, Newark, NJ, US, 2023 •Master of Science in Computational Fluid Dynamics, Cranfield University, Cranfield,UK, 2018 •Master of Science in Mechanical Engineering, Politecnico di Bari, Bari, Italy 2016 •Bachelor of Science in Mechanical Engineering, Politecnico di Bari, Bari, Italy, 2013 Major: Mechanical Engineering Presentations and Publications: F., P. A.S.F. Silva, P. Tsoutsanis, . F. Antoniadis, “Hovering rotor solutions by highorder methods on unstructured grids,,” Aerospace Science and Technology, Volume 97, 2020. F. Ricci, R. Vacondio, A. Tafuni; Direct numerical simulation of three-dimensional isotropic turbulence with smoothed particle hydrodynamics.” Physics of Fluids, 2023; 35 (6): 065148. F. Ricci, R. Vacondio and A. Tafuni, “A variable resolution SPH scheme based on independent domains coupling,” Proceedings of the 17th SPHERIC International Workshop, Rhodes, 2023. F. Ricci, R. Vacondio and . Tafuni, “High-order SPH schemes for DNS of turbulent flows,” Proceedings of the 2022 SPHERIC International Workshop, Catania, 2022. iv
To my grandfather and my nephew v
ACKNOWLEDGMENT First and foremost, I express my gratitude to my advisors, Prof. Angelo Tafuni and Prof. Renato Vacondio, for allowing me to pursue this Ph.D. program and for their constant support and guidance throughout all these years. Special thanks also to my defense committee, in the person of Prof. Farokhirad, Prof. Lieber, Prof. Marras and Prof. Dominguez for their feedback on my research. I would also like to thanks the financial support received by General Motors under Grant No. GAC3794, the National Science Foundation under Grant No. 2209793 and the department of Mechanical Engineering Technology. I’m also grateful for the help I received from all the researchers of the DualSPHysics group, in particular, the people of the SPH research group of the University of Manchester for their constructive criticism of this research during our meetings and to the ePhysLab for their technical support and for the opportunity of spending a period of stay at the University of Vigo. Most importantly, I want to mention my parents for all the sacrifices they have made for me since I was born and my beloved sister for being my confidant. Also, special thanks to my friends in Italy for showing their love and support even with an ocean between us and to the people in the US with which I shared these years far from home. vi
LIST OF FIGURES (Continued) Figure Page 4.19 Vorticity contours for ω= 1,5,10,20,30 at x=−0.5 for the (a) 2nd, (b) 4th and (c) 6th order schemes for the Taylor-Green Vortex at Re=1,600. (d) Reference solution in [253] . ...................... 79 4.20 Turbulent energy spectrum 2nd, 4th and ) 6th order schemes for the Taylor-Green Vortex at Re=1,600 with N= 2563particles. Numerical results are compared to the reference solution in [253]. ......... 80 5.1 (a) Example of two sub-domains Γ1and Γ2. (b) Buffer regions ∂Γ2 1and ∂Γ1 2with widths l∂Γ2 1= 2h1and l∂Γ1 2= 2h2are appended to their respective sub-domains. .......................... 82 5.2 Coupling procedure between sub-domains Γ1and Γ2. Buffer particles (orange) interpolate their properties over the fluid particles (light blue) in the coupled subdomain. ........................ 83 5.3 A buffer particle that moves into the fluid domain is transformed into a fluid particle (process 1). A fluid particle that enters the buffer region is transformed into a buffer particle (process 2). A buffer particle that moves outside the extended subdomain ∂Γj i∪Γiis deleted (process 3). 84 5.4 Particle insertion procedure: at each time step, the normal mass flux at the outer boundary of the sub-domain is calculated and added to the mass accumulation points (red squares). When the mass at the accumulation points reaches the reference particle mass, new particles (green) are created. ............................ 86 5.5 Sketch of the regularization procedure for buffer particles. The shifting correction is applied only in the direction tangential to the interface, while neglected in the normal direction n. In the corner region, this procedure is deactivated. ......................... 87 5.6 Call function for the DualSPHysics GPU code using a symplectic integrator 92 5.7 Call function for the main loop in the new multi-resolution algorithm. .96 5.8 Computational domain for the hydrostatic tank case ........... 98 5.9 Hydrostatic tank case: (a) density contours and (b) pressure distribution against the hydrostatic solution at t= 20s ............... 99 5.10 Computational domain for 2-D flow past a circular cylinder ....... 100 5.11 Dimensionless pressure for flow past a cylinder with (a) Re = 100 and (b) Re = 200 .................................. 101 xiii
LIST OF FIGURES (Continued) Figure Page 5.12 Dimensionless vorticity for flow past a cylinder with (a) Re = 100 and (b) Re = 200 .................................. 103 5.13 Streamlines for flow past a cylinder with (a) Re = 100 and (b) Re=200 .103 5.14 Time history of the drag and lift coefficient for a flow past cylinder with Re = 100 and Re = 200 .......................... 105 5.15 Time history of the lift coefficient for flow past a cylinder with (a) Re = 100 and (b) Re = 200. CLsindicates the lift coefficient obtained from a single resolution simulation while CLmbelongs to the multi-resolution simulation. ................................. 105 5.16 Sketch of the different SPH sub-domains created to resolve the flow past a cylinder at various Reynolds numbers ................. 107 5.17 Drag coefficient for flow past a cylinder with (a) Re = 1000, (b) Re = 3000 and (c) Re = 9500. These SPH solutions are compared against numerical results in [118]. ......................... 108 5.18 Vorticity contours for flow past a cylinder at Re = 1000 ......... 109 5.19 Streamlines for flow past a cylinder at Re = 1000 ............. 110 5.20 Vorticity contours for flow past a cylinder at Re = 3000 ......... 111 5.21 Streamlines for flow past a cylinder at Re = 3000 ............. 112 5.22 Vorticity contours for flow past a cylinder at Re = 9500 ......... 113 5.23 Streamlines for flow past a cylinder at Re = 9500 ............. 114 5.24 Time history of the lift coefficient for flow past an oscillating cylinder at Re = 100 for different amplitude Aand frequency Fratios: (a) (A, F) = (0.25,0.9), (b) (A, F) = (0.25,0.5), (a) (A, F) = (0.25,1.5), (a) (A, F) = (1.25,1.5) .......................... 117 5.25 Power Spectral Density (PSD) for the flow past an oscillating cylinder at Re=100 for different amplitude Aand frequency Fratio ........ 118 5.26 Contours of dimensionless vorticity for flow past an oscillating cylinder at Re = 100 with different amplitude Aand frequency Fratios ..... 119 5.27 Comparison of SPH aerodynamic forces against results in [194] for flow past an oscillating cylinder with A= 0.25 ................ 120 5.28 Computational domain for the propagation of regular waves case .... 120 xiv
LIST OF FIGURES (Continued) Figure Page 5.29 Comparison of free-surface elevation (a, b) and orbital velocities (c–f) between the multi-resolution and uniform resolution SPH simulations regular waves propagation ......................... 121 5.30 (a) Density and (b) velocity contours for the multi-resolution simulation for the wave propagation test case .................... 122 5.31 Time history of the mass variation in the multi-resolution simulation . . 123 6.1 Longitudinaland crosssections of the computational domain for the flow past a sphere. ............................... 125 6.2 (a) Time histories and (b) normalized power spectral densities of the aerodynamic coefficients for the flow past a sphere at Re = 300. .... 127 6.3 (a) Time history and (b) normalized power spectral density of the value of the streamwise coefficient at a x/D = 5.75 on the wake centerline for the flow past a sphere at Re = 300. ................... 128 6.4 Comparison of the average and root-mean-square values of the streamwise velocity along the wake centerline with the numerical results in [242] .128 6.5 Flow visualization by the Q criterion, colored by the velocity magnitude U, at every quarter of period from a view normal to the (x, z) plane for the flow past a sphere at Re = 300. ................... 130 6.6 Flow visualization by the Q criterion, colored by the velocity magnitude U, at every quarter of period from a view normal to the (x, y) plane for the flow past a sphere at Re = 300. ................... 131 6.7 (a) Time histories and (b) normalized power spectral densities of the aerodynamic coefficients for the flow past a sphere at Re = 500. .... 132 6.8 Normalized power spectral density for the streamwise velocity and pressure at (a)-(b) (x, y, z) = (2.5D, 0,0), (c)-(d) (x, y, z) = (2D, 0.3,0), (c)-(d) (x, y, z) = (2D, 0.0,0.3), for the flow past a sphere at Re = 500. 133 6.9 Flow visualization by the Q criterion, colored by the vorticity ωx, at every period T1= 1/St1from a view normal to the (x, y) plane for the flow past a sphere at Re = 500. ........................ 135 6.10 Configuration of the experimental setup for the three-dimensional dambreaking test case [113]. .......................... 136 6.11 Velocity magnitude contours for the three-dimensional dam break using the multi-resolution algorithm. The refinement region is highlighted in red. ..................................... 137 xv
LIST OF FIGURES (Continued) Figure Page 6.12 Comparison of the time history of the pressure at the front side measurements gauges between the present simulation and the experimental data [113]. 138 6.13 Comparison of the time history of the pressure at the top side measurements gauges between the present simulation and the experimental data [113]. 139 6.14 Comparison of the time history of the wave elevation at the measurements gauges between the present simulation and the experimental data [113]. 140 6.15 Snapshots of the density [kg/m3] contours across the center-line section in the longitudinal direction during the first impact of dam-breaking flows against the obstacle. ........................... 142 xvi
CHAPTER 1 INTRODUCTION 1.1 Background and Motivation The equations that model the motion of an incompressible viscous flow, i.e., the Navier-Stokes equations, are characterized by non-linear terms that make their mathematical treatment still an open challenge 100 years after their first derivation. Analytical solutions have been developed for specific cases by linearizing terms or reducing dimensionality, providing qualitative descriptions of flow in simple geometries. However, in order to obtain quantitative results, solving the Navier-Stokes equations numerically is crucial. Computational Fluid Dynamics (CFD) is the computational science encompassing all numerical methods used to analyze fluid motion and serves two main purposes. From a scientific perspective, it helps better understand complex phenomena like turbulence, multiphase flows, and flow-induced noise. In industry, it reduces costs associated with experimental aspects of the design process by narrowing the range of variables in experimental runs, optimizing designs, and shortening the time from conceptual design to production in various sectors such as automotive, aerospace, and energy. The earliest but still most popular numerical methods in CFD are mesh-based methods such as the Finite Difference Method (FDM) [92], the Finite Volume Method (FVM) [129], and the Finite Element Method (FEM) [284]. In mesh-based methods, the physical domain is discretized by computational nodes, which are topologically connected. Generating the computational grid is a complex and crucial part of the numerical computation workflow. It is often the most time-consuming task and demands 1
expertise to ensure the proper placement of computational nodes, considering the physics of specific phenomena while maintaining grid quality. High-aspect ratio or entangled grid elements can negatively impact the accuracy and stability of the numerical solution, particularly in the presence of intricate geometrical features [73]. Before applying the discretization specific to a particular numerical technique, the first choice concern the kinematic description of the continuum. Usually, two different formulations are employed [64]: in the Eulerian description, the computational nodes remain fixed in space and time, and convective terms are added into the governing equation, while the continuum moves and deforms with respect to the mesh, as opposed to the Lagrangian description, where each computational element is associated to a material particle thus following it during its motion. Both formulations have their advantages and disadvantages, depending on the particular physics of the problem: for fluid dynamics flows, the Eulerian description is usually preferred because it is conceptually simpler. However, it is generally unable to describe accurately moving interfaces. The Lagrangian formulation is preferred when dealing with structural mechanics because it is implicitly able to track interfaces and deal with a time-dependent stress-strain relationship. However, as opposed to the Eulerian description, when there is large material deformation, the motion of the computational nodes can result in distorted or entangled elements that can deteriorate the accuracy and the convergence of the numerical methods. Due to the successes in predicting aerodynamical flows, the CFD technique has also been extended to fluid-structure interaction (FSI) problems, where movable or deforming objects interact with surrounding fluid flows. This class of problems is characterized by strong non-linearities and a multi-physics nature, making the numerical solutions in a single mathematical framework rather challenging. The range of applications of FSI problems spans different engineering areas, including aeroelas2
ticity [104], bio-medical flows [96], biological flows [241], structural engineering [124], coastal and marine applications [278] among the others. Many different CFD methods have been proposed to tackle FSI problems. Among those, the two most popular techniques are Arbitrary Lagrangian-Eulerian (ALE) [93] and the Immerse Boundary Method (IBM)[193]. The idea behind the Arbitrary-Lagrangian-Eulerian (ALE) description is to retain the advantages of both the Eulerian and the Lagrangian formulation while mitigating their drawbacks. Usually, in ALE methods, a body-fitted grid discretizes the solid domain, while an Eulerian description is employed for the fluid far region. Instead, near the interface, an arbitrary velocity is imposed on the mesh to avoid excessive mesh distortion, and convective fluxes account for this arbitrary motion in the governing equations. However, in the presence of large displacement, the re-meshing and remapping process is unavoidable, along with the associated computational cost and the challenges in remapping the solution over the new mesh, which can deteriorate the accuracy of the numerical solution. In the Immersed Boundary Method, non-conformal meshes are used to tessellate the fluid and the solid region, which are usually described using a different kinematic description [159]. The ability to employ overlapping grids, avoid ALE methods’ computationally cumbersome re-meshing process and greatly simplify the grid generation process. The coupling between the non-conformal meshing is usually achieved by adding a local volumetric force into the governing equations of the fluid part, using a smoothing function to redistribute the effect of the interface over a range of computational nodes. However, different strategies have been proposed [87]. Another coupling approach is the Cut-Cell Finite-Volume approach [275], which doesn’t use a volumetric forcing and has better mass and momentum conservation properties. Nevertheless, two main drawbacks characterize the IBM: the first one, as discussed in [25], concerns the ”added-mass” effect when dealing with a high 3
fluid-to-solid density ratio; the second one regards the simulation of a high Reynolds number flows, for which a mesh refinement procedure is required in order to ensure an adequate resolution near the solid boundaries [254]. Another class of problems that involve moving interfaces is free-surface and multiphase flows problems. In this case, two main strategies are used in mesh-based methods, front-capturing and front-tracking approaches [243]. The Volume of Fluid (VOF) technique [94] is the most popular front-capturing method. In the VOF, the interface between two fluids is represented by a color function, which value is based on the volume fraction of each phase in a particular computational element. The method is composed of two steps: in the first one, the interface is reconstructed from the value of the color function in each cell. Early studies used a Simple Line Interface Construction (SLIC) [181, 94], representing the interface by segments aligned with the mesh. Although it is very simple, this method leads to large interface smearing. More accurate approaches, such as the Piecewise Linear Interface Construction (PLIC) [207], significantly improve the method’s accuracy. The color function is advected using the velocity field in the second part of the approach. Despite the mathematical formulation ensuring good conservation properties, the discontinuous reconstruction can lead to instabilities in the presence of high-curvature interfaces. Another approach is the Level-Set method [190], in which the interface is represented by a smooth function that moves according to an advection equation. The advantage of the Level-set method is that with respect to the VOF, the function is smooth; however, it has a worse conservation mass property with respect to the former method. In front-tracking techniques [244], instead, the boundaries between different phases are represented by a set of marker points that are connected and move along with the fluid. The drawbacks of these methods are the additional data structure to 4
describe the front and the explicit treatment that require the topology change of the interface. Opposed to grid-based methods, in meshless-based methods the approximation of the governing equations is built upon a set of computational nodes for which grid connectivity is not specified. Among them, there are the Smoothed Particle Hydrodynamics (SPH) method [80], the Meshless local Petrov-Galerkin method [12], the Diffuse Element Method [180], the Element-Free Galerkin methods [19], and the Moving-particle semi-implicit method [116]. The SPH method is arguably the most popular for simulating incompressible viscous flows. In the SPH method, the Navier-Stokes equations are approximated upon a set of discrete particles that carry the physical properties of the fluid. The interpolation is based on the convolution with a kernel smoothing function. One of the advantages of the SPH method is its robustness; in fact, as demonstrated in [24], the SPH formulation is consistent with a variational approach, ensuring conservation of all relevant physical quantities. Moreover, due to the Lagrangian formulation, the fluid properties are advected exactly and able to implicitly describe interfaces, suach as free-surface. The SPH rEsearch and Engineering International Community (SPHERIC) is an international organization that groups the community of SPH researcher and industrial practitioners. The main objective of SPHERIC is to steer the research focus on the SPH method. In fact, the steering committee has identified five different aspects of SPH that need to be addressed in order to encourage the widespread adoption of the method for CFD study. The SPH Grand Challenges are [245]: Convergence, consistency and stability, Boundary conditions, Adaptivity, Coupling to other models, Applicability to industry. One topic that has been poorly addressed by the SPH research and limits the range of applicability and the fidelity of the method concern the simulation of turbulent flows. In this work is studied the issue of the 5
inclusion of turbulence effects in the SPH method. In particular, a major focus is given to the relationship between turbulence modeling and the issue of adaptivity within the SPH method. 1.2 Aim and Objectives The main goal of the present dissertation is to extend the range of applicability of the Smoothed Particle Hydrodynamics method, with a focus on turbulent flows. To this end, the present objectives are set: 1. Gain insight into the performance of the SPH method in the computation of turbulent flows. 2. Analyze the effect of the discretization error due to the discrete operator, the density diffusion term and the particle shifting on the numerical computation of isotropic turbulence. 3. Assess the performance of high-order kernel scheme in the numerical solution of isotropic turbulence 4. Develop a novel and highly-efficient variable-resolution approach for turbulence problems where high local resolution is required. 5. Implement the new algorithm in the DualSPHysics open-source code and validating across different test cases. 6. Extension and validation of the approach to the numerical computation of threedimensional flows. 1.3 Structure of the Dissertation This dissertation is structured as follows: 1. In Chapter 1 is presented an overview of the numerical approaches for the computation of flows characterized by moving interfaces 2. Chapter 2 present a literature review over the state of the art of the SPH method with a particular focus over the turbulence and its modeling within the SPH method and the issue of adaptivity. 3. Chapter 3 present the mathematical basis of the SPH along with the numerical technique adopted in this work. 4. In Chapter 4 the numerical computation of homogenous isotropic turbulence within the SPH method is adressed. The effect of the discrete operator over the accuracy is studied, and the effect of a Density Diffusion Term and the Particle Shifting Technique is assessed. The numerical simulation results with high-order SPH schemes of decaying isotropic turbulence are presented. 6
approach where the flow is characterized by complex moving interface,e.g., free-surface or moving objects. As previously discussed, one of the primary limitations of the Smoothed Particle Hydrodynamics (SPH) method, which initially hindered its broad application in engineering, is its high computational cost. This is due to a larger computational stencil (typically on the order of 30+ and 300+ particles for 2-D and 3-D simulations, respectively) compared to Finite Volume Method (FVM) or Finite Element Method (FEM). However, with the rise of massively parallel architectures, such as those based on Graphics Processing Units (GPUs), numerous in-house or open-source ([62, 21, 34]) codes have rapidly developed. Originally, GPUs were utilized in graphics applications, such as image/video processing and video gaming. Their architecture is built around the Streaming Multiprocessor unit, composed of several Arithmetic Logic Units (also known as CUDA cores). The latest generations of GPUs have hundreds of streaming multiprocessors, enabling them to perform thousands of arithmetic operations simultaneously. The SPH method, particularly in its Weakly-Compressible formulation utilizing an explicit time-scheme, is especially suited for such architectures. This is due to the high arithmetic intensity of particle-particle interaction calculations which, as noted in Dominguez et al. [59], generally represents the most computationally expensive operation in an SPH simulation. In the SPH method, the computational stencil is computed through the creation of a neighbor list. Dominguez et al. [59] have discussed the primary techniques for creating this neighbor list, emphasizing the importance of particle reordering to ensure optimal coalesced memory access. Additionally, in Dominguez et al.[61], various optimizations are detailed to further enhance the efficiency of a GPU implementation of an SPH model. 13
2.2.2 Convergence and accuracy of SPH method The discretization error in the SPH method is composed of two contributions: the first one is due to the smoothing procedure that, for a conventional smoothing kernel function, has an order of O(h2), where his the smoothing length. The second source of error is due to the discrete approximation and is a function of both the smoothing length and the number of neighbors Nbincluded in the kernel support. Monaghan [164] conjectured that the SPH method because the discrete pressure gradient tends to arrange the particles in a glass-like configuration, presents more favorable convergence properties than Monte-Carlo methods. In Zhu et al. [283] is noted that to maintain constant the rate of convergence, one must have h→0, Nb→ ∞ and N→ ∞ simultaneously, and propose a power-law to relate hand Nbto the total number of particles Nin the computational domain for a typical 2nd order smoothing kernel function as: h∝N−1/6, Nb∝N0.5(2.1) The same conclusions were reached in [199], where the discretization error of a 1D SPH approximation was analyzed through a second Euler-McLaurin summation formula. The immediate consequence of these findings is that to preserve the second-order accuracy, the computational stencil, which is already large in comparison to meshbased methods, must grow as hshrinks, increasing the computational cost. Besides, early SPH practitioners run into the so-called ”pairing instability” over a certain threshold of neighbors. As demonstrated in Dehnen and Aly[55], this phenomenon is caused by negative values in the Fourier transform of the smoothing kernel. For this reason, Wendland kernels [265] have become the standard in the SPH method. Another widespread technique to decrease the inaccuracies due to a disorder distribution is the Particle Shifting Technique (PST) [136], which aims to restore a 14
more regular distribution by moving the position of the particles, typically modeling the displacement with a Fick’s law based on the concentration gradient. A similar approach has also been proposed within the SPH-ALE formulation [187]. The PST has soon become a cornerstone when dealing with incompressible viscous flows, although it requires careful treatment in the presence of a free surface due to the zero-th error introduced in the SPH gradient operator by the truncated support. In this case, the formulation of the PST is usually modified by eliminating the normal component to the free surface of the shifting vector [136]. Different procedures have been proposed in order to identify the fluid particles that belong to the free-surface: in the method proposed in Lee et al.[125], the free-surface particles are identified through the value of the divergence of the position vector. A more computationally costly but accurate approach has been proposed in Marrone et al. [152], where the minimum value of the eigenvalues of the renormalization matrix [202] and an additional procedure, based on scanning the ”umbrella-shaped” regions are used to determine the particles belonging to the free-surface. Recently, different works have addressed the inconsistency introduced at the free-surface from the shifting algorithm [109, 262, 145, 120]. Additionally, iterative explicit and implicit shifting algorithms [203] have also been proposed to ensure a better regularization of the particle distribution. Another aspect closely related to the convergence property is the consistency of the SPH method, in particular in the presence of boundaries that truncate the support domain of the smoothing kernel. In that case, neither the zero nor the first-order consistency is ensured. Different numerical techniques have been proposed to correct the inconsistency: approach based on the renormalization matrix [202, 101], Moving Least-Square schemes [58], the Corrective Smoothed Particle Hydrodynamics Method (CSPH) [36], the Finite Particle Method [140] (FPM), the modified Smoothed Particle Hydrodynamics method (MSPH) [16]. 15
The downside of these approaches is the computational cost associated with the solution of the linear system, which increases steeply with the dimensionality and the order of consistency required. Besides, most of these methods break the symmetricity of the particle interaction, with the loss of the conservation properties of the method, although some authors [185] suggest that this property can be relaxed. Different strategies have been proposed in the last years to achieve a higher convergence rate. Using a Riemann-SPH scheme, in Avesani et al. [13] has been proposed a WENO-reconstruction where the polynomials are obtained through an MLS interpolation. The reconstruction stencils are defined by partitioning the support domain of the kernel into different sectors. This approach has been furtherly improved in Antona et al. [7], where the FPM method is used in place of the MLS scheme. This scheme has also been extended in Avesani et al. [14] to include a high-order space-time reconstruction with an ADER-WENO-SPH scheme. The drawbacks of this approach are that the convergence rate is still limited by the SPH approximation in the particle interaction and the computational cost associated with the MLS reconstruction that requires a matrix inversion for each stencil Lind and Stansby [135] have shown that it is possible to achieve high-order convergence rates by employing high-order kernels. These kernels are obtained by relaxing the non-negative property of the smoothing kernel. The shortcoming of this approach is the kernel must be very well sampled, restricting the range of applicability to a uniform distribution of particles. Following this approach in Nasar et al. [178], high-order Dirichlet and Von Neumann boundary conditions are proposed, while in Nasar et al. [177], a new formulation based on a kernel consistency correction has been proposed to limit the smoothing error due to the discrete SPH operator when high-order kernels are employed. 16
2.2.3 Weakly-Compressible vs. incompressible SPH There are two main approaches to treating incompressible flows in SPH: the WeaklyCompressible SPH (WCSPH) and the Incompressible SPH (ISPH). The latter approach has been proposed firstly in Cummins and Rudman [52], where a Chorin’s projection method is used to enforce an incompressibility of the flow through the solution of a pressure Poisson equation that arises from the divergencefree velocity field assumption. Alternative algorithms have been also proposed in Shao and Lo [221] and Hu and Adams[97]. Notably, it has been shown in [230, 231] that the ISPH and the moving particle semi-implicit (MPS) method are equivalent. The ISPH has been applied to free-surface flows [108, 107, 228, 128], offshore application [133], and multi-phase flow [132]. The advantages of the ISPH against the WCSPH are the ability to obtain a smoother pressure field, solve the well-known problems that afflict the WCSPH, and that is able to use a larger time step. However, while the latter approach is computationally efficient due to the explicit scheme that is more easily parallelizable, especially by exploiting a parallel framework with a SIMD paradigm, such as OpenMP and CUDA, the ISPH presents bottlenecks that limit the efficiency of the method. The lack of topological connectivity between computational nodes, one of the distinctive traits of the SPH method, forces the construction of the PPE matrix at every time step. Furthermore, because of the large computational stencil, the memory requirements for storing the PPE matrix are considerably larger than the WCSPH, and also, with respect to the latter formulation that enforces the free-surface boundary condition implicitly, the ISPH must rely on free-surface detection method for correctly apply the Dirichlet boundary conditions to close the PPE. In the WCSPH the density and the pressure are coupled through a stiff equation of state, usually Tait’s Equation. To avoid an excessive restriction of the time-step due to the CFL condition, the Mach number is taken around 0.1, which bound the 17
variation of the density within the 1%. The advantage of this approach is the ability to use an explicit time-stepping scheme. One of the drawbacks of the WCSPH formulation is the presence of highfrequency oscillation that affects the smoothness of the density field. This aspect has been studied by several authors [112, 68], and its stems from the combination of two different factors: on the one hand, there is the employment of a stiff equation; on the other hand, there is the collocation nature of the SPH scheme. Various technique have been proposed during the past two decades to address this issue. Several authors [255, 192, 98] have proposed an approach based on the definition of a local Riemann problem to solve the particle interaction. The RiemannSPH method has been used with different limiters([200, 174, 99, 210, 117], however, for violent free-surface flow, this approach has revealed to introduce too much dissipation. In Ferrari et al. [72], has been proposed a diffusive term obtained applying a Rusanov Flux to the Riemann-SPH approach, but introducing less dissipation with respect to the latter approach. One of the drawbacks of this term is that it is neither consistent, so that doesn’t vanish for h→0, nor preserve the hydrostatic solution. To restore the consistency in Molteni and Colagrossi[162] is proposed a density diffusion term based on the discretization of the laplacian of the density field with the Morris formula, which has been improved in Antuono et al. [10], the so-called ”δ-SPH scheme, in order to preserve the hydrostatic solution with a consistency correction at the free-surface based on the calculation of the renormalization matrix. In Green et al. [85] is presented a diffusive term based on the application of a Roe’s approximating Riemann solver, and also it is shown that the δ-SPH can be viewed as a particular case of this model. More recently, in Fourtakas et al. [75], a new diffusive model is proposed, based on the neglection of the hydrostatic pressure in the calculation of the Laplacian, which is able to preserve the hydrostatic solution without relying on the calculation of the renormalization matrix, reducing the computational cost. 18
2.2.4 Boundary condition in SPH The imposition of boundary conditions to close the numerical problem is a challenging topic within the SPH method, and it has been listed as one of the SPH ”Grand Challenges” by the SPHERIC committee. Among the reasons are the lack of the Delta Kronecker property and the inconsistency of the SPH interpolation in the presence of boundaries that truncate the support domain of the smoothing length. For the definition of inlet/outlet boundary conditions, the most popular approaches are based on the extension of the computational domain by ”buffer regions” in which the fluid properties are imposed to enforce Dirichlet BC or extrapolated by the fluid domain for the definition of Neumann BC. They also model the inflow and outflow by generating or deleting fluid particles that enter the computational domain. Different strategies have been proposed based on this approach, most notably [122, 247, 69, 238]. Regarding solid boundary modeling, various approaches have been proposed in the literature, each with advantages and shortcomings. In Monaghan [165], the solid BC was imposed by means of solid particles that exerted a repulsive force over the fluid particles to enforce the no-penetration condition. The repulsive force was modeled based on the Lennard-Jones potential; however, this formulation was unable to simulate a smooth surface, imposing an implicit roughness with a spatial scale equal to the particle distance that results in a disordered configuration of the fluid particles close to the interface. A further modification to address this issue was proposed in [171, 170] based on the definition of the normals to the solid interface to ensure that a particle moves in the parallel direction of the solid boundary experience a constant force. The enforcement of the no-slip condition is then implicitly accounted for by including the viscous term in the calculation of the interaction forces between the fluid and solid particles. 19
A different strategy, the ”Dynamic Boundary Conditions” method, has been proposed in Dalrymple and Knio [53], where typically one layer of particles is placed to describe the solid boundaries. The advantage of this method is its computational efficiency and capability to discretize complex domain because these solid particles behave as fluid particles when calculating their density value. As studied in Crespo et al. [50], the solid particles exert a force that depends on the distance and the pressure of the incident fluid particles. This approach has been used to simulate the interaction between incident waves and coastal study. The shortcomings of this approach are, however, the unphysical gaps between fluid and solid particles and the generation of large oscillations in the density field. Besides, the no-penetration condition is not explicitly enforced. In the ”Ghost Particle Approach” [149], as in the DBC, several layers of particles are created at the beginning of the simulation to represent the solid interface. To generate these boundary particles, the solid interface is represented by a piecewise linear function, usually a spline, and discretized by a set of particles with a spacing equal to the characteristic particle size of the problem. After the normal and the tangent to this interface are calculated, the first layer of solid particles and the associated interpolation points in the fluid domain are created by a simultaneous contraction and expansion of the solid surface. This process is repeated recursively to ensure a sufficient number of layers based on the width of the support domain of the kernel smoothing function. The fluid properties for the solid particles are then retrieved, interpolating over the fluid domain at the interpolation points defined in the procedure outlined above. In Marrone et al. [149], the interpolation is carried out using an MLS technique, while in English et al. [66], where the mDBC is proposed to address the issue with the DBC, the correction technique of Liu and Liu [140] to restore particle consistency is employed. 20
One shortcoming of the Ghost Particle method is that complex geometrical features, i.e., sharp angles, must be carefully treated to define the interpolation points in the fluid domain correctly. Moreover, in the case of submerged thin elements, this approach requires placing sufficient layers on both sides of the element, which can lead to an unacceptable number of particles without a variable-resolution approach. In Adami et al. [1], the pressure is assigned based on the force balance at the solid interface. Another popular approach is the ”mirroring ghost particles”, where boundary particles are generated by mirroring with respect to the solid interface, the position of fluid particles close to the contours, that either carry the field properties or obtain their values through interpolation. However, this approach is more computationally cumbersome with respect to the ”Ghost Particles approach” because the mirroring procedure must be executed at each time step. Moreover, it is difficult to handle complex 3-D geometries. The Virtual Boundary Particle [72] uses a different strategy to ensure the consistency of the SPH interpolation at the boundaries. The solid interface is discretized by a set of boundary particles whose purpose is purely geometric. An interior fluid particle close to the solid contour generated a set of fictitious particles using a local point-symmetry instead of a plane-symmetry employed in the mirrored particle approach. This procedure has been further improved in [246, 76] to ensure that the local stencil resembles that of an interior particle with full support and to have a better definition of the local stencil near corner regions. Based on this approach, in Fourtakas et al. [75], the Local-Uniform-Stencil boundary condition has been presented: for each interior particle, at the beginning of the simulation, a uniform stencil of virtual particles is created. This local stencil moves alongside the fluid particles associated. Triangular elements then discretize the solid interface, and at each time step, using a raycasting algorithm, the particles 21
of the uniform local stencil are identified in the boundary region. These particles are then activated and are used to enforce the solid boundary condition. A uniform stencil ensures a zero-th and first-order consistency. Another alternative is based on accounting for the truncated kernel at the boundary using surface integrals. These integrals calculate a corrective factor that enters the governing equations. This approach, called ”Semi-analytical wall boundary condition,” was first proposed in Kulasegram et al. [121] and then developed in [148, 71, 155]. However, this approach is not suited for an efficient implementation on GPUs. 2.3 Turbulence and its Modeling in SPH 2.3.1 Turbulence and its statistical modeling The study of turbulence is one of the most complex aspects of fluid dynamics. Its statistical rather than mathematical modeling is also crucial due to the ubiquitousness of turbulent flows in any real-life flows of practical interest. Despite that the complex nature of turbulence makes a formal definition difficult, an attempt can be made to characterize a turbulent flow by the following properties: 1. chaotic 2. has a large and continuous spectrum of spatial and time scales 3. three dimensional 4. very sensible to change in initial and boundary conditions 5. intermittent 6. mixing and dissipation are enhanced with respect to laminar flows. It must be stressed that the words ”chaotic” and ”random” are not synonymous to highlight that turbulence arises from a non-linear dynamical system, the NavierStokes equations, which is deterministic. 22
numerical results with the log-wall of the turbulent boundary layer. In Violeau and Issa (2007) [257], approaches such as the k-ϵand the Explicit Algebraic Reynolds Stress Models (EARSM) have been investigated to address the simulation of dambreaking flows. Therein, the authors report a good qualitative agreement against experimental results, especially for the EARSM model, which is more suited to study regions with high distortion near the free surface. Future contributions led to a k-ϵmodel combined with the semi-analytical wall boundary conditions for Weakly Compressible SPH (WCSPH) in Ferrand et al. (2013) [71] and Incompressible SPH (ISPH) in Leroy et al. (2014) [128], showing good agreement with the Finite-Volume Method (FVM) for a fish-pass flow. The k-ϵmodel has also been used with the ISPH to investigate wave overtopping [219] and breaking [220] and more recently solitary[260] and periodic waves[261]. In later years, a switch to Large Eddy Simulation (LES) has been observed, trying to exploit the analogy between the LES filtering procedure and the SPH interpolation. One of the pioneering contributions to LES applied to SPH can be found in Lo and Shao (2002) [142], where turbulent solitary beach waves have been studied with a focus on turbulence during the breaking phase. This model has then been extended to WCSPH by Dalrymple and Rogers (2006)[54], providing validation against numerical experiments for cases of wave overtopping and beach waves in 2-D and dam breaking flow in 3-D. In Mayrhofer (2015) [156], a LES approach in SPH has been coupled with the semi-analytical wall boundary conditions for studying turbulent channel flow, however, results therein have shown an overprediction of the streamwise velocity. The authors pointed out that this is probably due to insufficient resolution when capturing the vortex. More recently, the LES approach has been employed in the simulation of turbulent open channel flows over and within natural porous gravel beds [105] and in the modeling of oil spill[225] using the ISPH formulation. 29
In Di Mascio et al. (2017) [57] and Antuono et al. (2021) [9], a new approach to LES modeling in SPH is presented. The key idea of these contributions is to use a time-space filtering procedure, where the SPH kernel operator acts as a spatial filter while time filtering is implicitly accounted for with additional terms in the governing equations. This strategy offers a consistent approach to LES while interpreting the δ-SPH by Molteni and Colagrossi (2009) [162] from a LES perspective, i.e., the delta coefficient as the deviatoric strain of the mean flow. The δ-LES model has also been employed for studying gravity waves [158] and dam-break flow [158]. From a Direct Numerical Simulation (DNS) perspective, Robinson and Monaghan (2012) [208] have attempted the DNS of two-dimensional decaying turbulence, although 2-D turbulence is fundamentally different from 3-D because of the inverse energy cascade phenomena [22, 119]. Notable contributions have been made in Mayrhofer et al. (2015) [156], where a wall-bounded flow is simulated, and a lower bound in terms of the number of particles per vortex is suggested. From the presented literature review, it is evident that several crucial aspects remain to be addressed concerning turbulence modeling in SPH. The initial consideration involves determining the most suitable approach for turbulence modeling. While the RANS modeling appears to be a more favorable choice, given the high computational cost associated with the SPH method, it seems to lack the necessary theoretical foundations for its application in the typical SPH domain, which includes free-surface flows with considerable non-stationarity and distortions. LES modeling appears to be the most suitable approach, but there are still problematic aspects related to the SPH method that require further investigation. The first issue concerns the method’s order of accuracy: LES modeling demands a high degree of accuracy to ensure that the numerical dissipation introduced by the numerical scheme does not compromise the fidelity of the results. However, this conflicts with the current convergence order of the SPH method, which is 30
at most second-order in optimal situations, while under normal conditions, due to the truncation error caused by the irregular distribution of particles, it ranges between 1 and 2. Although approaches for achieving a convergence order higher than second-order have been proposed, as previously discussed, they do not yet seem to be mature enough for broader application. The second aspect pertains to the computational cost of the method. As previously discussed, a turbulent flow is characterized by a wide spectrum of spatial and temporal scales. In grid-based methods, this aspect is addressed by increasing the resolution in areas where vortex development is expected, particularly near solid surfaces. The introduction of variable resolution in SPH encounters challenges arising from the method’s particle-based Lagrangian formulation. In the following section, these aspects will be discussed, and an overview of the state-of-the-art for multi-resolution in the SPH method will be provided. 2.4 Adaptivity within the SPH Method One weakness of the SPH approach is that adopting a multi-resolution formulation is more challenging than in mesh-based method. There are several reasons: •As previously illustrated, in order to preserve the convergence property, the ratio between the smoothing length and the mean particle size must be kept at least constant. This means that in order to increase the numerical accuracy, the value of the kernel smoothing length cannot be reduced without increasing the total number of particles in the simulation. •The interaction between particles with different smoothing lengths must be treated carefully. Different approaches have been proposed in the literature, and some of them are derived trying to preserve the variational consistency [24, 26]. Nevertheless, the smoothed formulation of the method prevents a sharp variation of the particle size. •The symmetricity of the smoothing kernel function, a cornerstone of the SPH formulation, implies the isotropic distribution of the computational nodes, as opposed to grid-based methods that can exploit the anisotropy of the flow by using different spacing depending on the spatial coordinate. A typical example is the inflation layers used to discretize the near-wall region, where the spatial resolution in the wall-normal direction is usually one order of magnitude lower than in the streamwise or spanwise direction. 31
Historically, early attempts for the introduction of adaptivity in the SPH method were focused on the introduction of a variable smoothing length formulation coupled with the definition of regions with different particle sizes at the beginning of the simulations, for example in Bonet and Paz [24] for the collapse of a circular dam over a surface, cylindrical water blast [26], wedge water entries [184], heaving cylinder and cone in traveling waves [188][189]. However, these approaches were limited to problems with a short time scale and were unfeasible for cases in which the motion of the computational nodes highly distorted the initial configuration of particles. Instead, later efforts aimed to dynamically increase the particle resolution by introducing a dynamic refinement process under the constraints of conserving mass, momentum, and angular velocity and minimizing the error in the estimation of the density with different splitting patterns [70, 204]. The basic concept is to minimize, in a least squares sense, the error in the SPH estimation of the density field between the original particle distribution and the refined one. A parametric study on the optimal splitting pattern, along with the optimal value of the weight of each particle size λi, their smoothing length value hi, and the distance with respect to the original particle ϵi, has been presented in Vacondio et al. [248]. The dynamic splitting procedure has been coupled in Vacondio et al. [249] to a de-refinement process based on coalescing pairs of particles with similar sizes and extended to 3-D in Vacondio et al. [250]. A similar approach is used in Yang et al. [272] for the computation of multi-phase [273] and free-surface [274] flows, where the splitting criterion is not based on the geometric definition of a fixed refinement region, but instead on the position with respect to the free-surface or interface between different phase. However, one weakness of this approach is the loss of computational efficiency with the increase of the refinement ratio. In fact, as discussed in Vacondio et al. [248], in order to minimize the error introduced by the splitting procedure, the value of 32
the smoothing length between the coarse and the fine particles remains almost equal, while the particle size decreases due to the splitting procedure (typically hexagonal or square). This means that the computational stencil increase along with the resolution, introducing an important source of inefficiency of the method, which is even more serious in three-dimensional applications. Moreover, the coalescing procedure is only performed pairwise, so less frequent with respect to the splitting procedure, which causes an unnecessary overhead in terms of the number of particles. In [175] [89], this issue is tackled by taking as kernel smoothing length, instead of the optimal value prescribed by the least square minimization, the average value in the support domain of the kernel, and results are presented for the computation of flow past bluff bodies. A different approach is proposed in Barcarolo et al. [15]: as in the splittingcoalescing procedure, a particle is split into fine particles when entering the refinement region. However, in this approach, the original particle, instead of being deleted, is retained and it is ipotetically advected. A weight function governs the transition between the two zones to avoid pressure discontinuities at the interface. The same approach has been improved in Chiron et al. [37] by resolving the interaction between coarse and refined particles with the definition of buffer zones that avoid the interaction between particles of different sizes. Using this approach, in Sun et al. [235], results are presented for flows past bodies with various shapes, coupling the multiresolution algorithm with a Tensile Instability Control (TIC) term, while in Sun et al. [236] the method is employed to study water entry of circular cylinders. In Gao et al. [78], a block-based adaptive particle refinement algorithm is proposed for dynamically changing the particle refinement domain coupled with a regularization process similar Particle Shifting Technique (PTS) [136] to obtain an isotropic distribution of the refined particle in the transition zone. However, one of the weaknesses of the Adaptive Particle Refinement approach is the decoupling of the position between coarse and fine (daughter) particles, which 33
in highly distorted flows, as detailed in Chaneac et al. [35], results in an over-creation of particles which eventually leads to unstable solutions. In the domain-decomposition method, the computational domain is subdivided into many computational sub-problem, and are advanced in time with an appropriate coupling strategy. An approach based on this formulation has been proposed in Bian et al. [20]. However, the results presented only two levels of refinement, and the density formulation wasn’t able to treat free-surface flows. In Shibata et al. [224] a similar strategy based on inlet/outlet boundary conditions is introduced for the coupling of different resolution zones in the context of the MPS method. Multiresolution approaches, mainly devoted to Fluid-structure problems in which different resolutions are defined between the fluid and the solid phase, have been proposed in Zhang et al. [281] and Khayyer et al. [110]. However, these approaces are limited only to fluid-structure problem and allow a very small variation in the particle size between the different phases. 34
CHAPTER 3 NUMERICAL FRAMEWORK This chapter presents the basis of the mathematical treatment in the SPH method alongside the computational techniques that represent the state-of-the-art of the SPH methodology and that have been used in this work. 3.1 Basis of the SPH Methodology 3.1.1 SPH continuous interpolation The mathematical treatment of the SPH interpolation method starts from the convolution of a field function and the Dirac’s Delta function δ(x−x′): f(x) = ZΩ f(x′)δ(x−x′)dx′.(3.1) An approximation of the previous identity is obtained by substituting the δfunction with a smoothing kernel function W(x−x′, h): ⟨f(x)⟩=ZΩ f(x′)W(x−x′, h)dx′,(3.2) where his the smoothing length, a parameter that defines the size of the kernel support Ω. The properties of the kernel smoothing function Winfluence the convergence, accuracy, and stability of the SPH interpolation and will be discussed in the next section. Expressing the derivative of (3.2) with respect to x′and using integration by parts, it can be derived the expression for the SPH gradient operator: ∇⟨f(x)⟩=ZΩ∇f(x′)W(x−x′, h)dx′−ZΩ f(x′)∇W(x−x′, h)dx′= =Z∂Ω f(x′)W(x−x′, h)·¯ ndS −ZΩ f(x′)∇W(x−x′, h)dx′. (3.3) 35
The first integral is obtained by applying the Gauss theorem to pass from a volume to a surface integral and, assuming that the kernel function has compact support and this is not truncated, this term is equal to zero, and the following identity is obtained: ∇⟨f(x)⟩=−ZΩ f(x′)∇W(x−x′, h)dx′.(3.4) Moreover, if the smoothing kernel W(x−x′, h) is an even function, the last expression can be rewritten as: ⟨∇f(x)⟩=ZΩ f(x′)∇W(x′−x, h)dx′.(3.5) where the differential operator ∇is referred to x. From here on, the bracket notation to identify the SPH approximation will be omitted. 3.2 Properties of the Smoothing Kernel Function There is a set of properties there are desirable for the kernel smoothing function: 1. Unity: ZΩ W(x−x′, h)dx′= 1.(3.6) 2. Even function: W(x−x′, h) = W(x′−x, h).(3.7) 3. Compactely supported: W(x−x′, h) = 0 if |x−x′|> ah. (3.8) 4. Positivity: W(x−x′, h)≥0,∀x′.(3.9) 36
5. Delta function: lim h→0W(x−x′, h) = δ(x−x′).(3.10) 6. Monotonically decreasing 7. Smoothness The first condition ensures that the SPH interpolation has a C0consistency when the support domain of the kernel function is far from the boundaries and is able to reproduce a constant function exactly . Together with the requirement that the smoothing kernel is an even function, it can be demonstrated that the SPH interpolation in Equation (3.2) is converging with a 2nd order convergence rate. In fact, expanding a field function f(x′) in Taylor series: f(x′) = f((x) + f′(x)(x′−x) + 1 2f′′(x)(x′−x)2+O(x′−x)3,(3.11) and multiplying the previous equation by Equation (3.2), it is obtained: f(x) = f(x)ZΩ W(x−x′, h)dx′−f′(x)ZΩ (x−x′)W(x−x′, h)dx′ +1 2f′′(x)ZΩ (x′−x)2W(x−x′, h)dx′+ZΩO(x′−x)3dx′. (3.12) From the previous identity it can be seen that in order to ensure that the approximation has a 2nd order of accuracy, the first two moment Mkmust be equal to 0: M0=ZΩ W(x−x′, h)dx′= 1,(3.13) M1=ZΩ (x−x′)W(x−x′, h)dx′= 0.(3.14) It can be observed that these two conditions are fulfilled if the smoothing kernel W(x−x′, h) verifies the Unity and the Even properties in Equations (3.6) and (3.7). 37
Moreover, it can also be demonstrated that all odd moment Mkare identically equal to 0. Regarding the other properties: the smoothing function is chosen to have compact support to reduce the numerical stencil and the computational cost; the positivity condition, rather than a mathematical requirement, is based on the physical admissibility of hydrodynamical states, i.e., density and energy. This condition also prevents the vanishing of the even-order moments of the smoothing function so that the SPH approximation has, at most, a 2nd-order accuracy. The Delta function in Equation (3.10) properties ensure that the smoothing function recovers the Dirac distribution as the smoothing length htends to zero. The sixth property is based on the assumption that closer value must have a bigger influence over the SPH estimation, while the last property improves the accuracy of the interpolation [199]. Another interesting property that derives from the symmetric condition, is that the kernel function can be expressed as a function of the distance between spatial points: W(x−x′, h) = W(|x−x′|, h).(3.15) Furthermore, the kernel function can be written in dimensionless form as: W(|x−x′|, h) = W(q),(3.16) where qis obtained operating the following change of variable: q=|x−x′| h.(3.17) Therefore, a general expression of kernel function is: W(q) = αd hdf(q),(3.18) 38
discussed in [198], this discretization ensure a more isotropic particle arrangement and avoid the clumping of particles. Besides, as shown in [259], the set of discrete operators chosen is consistent with a variational formulation, ensuring the proper conservation of the linear and angular momentum. Nevertheless, some authors have proposed alternative formulation [235, 185], based on using both the antisymmetric and symmetric SPH operator in the discretization of the momentum equation, depending on the sign of the pressure. However, these approaches don’t explicitly preserve momentum conservation. 3.5 Viscosity Models 3.5.1 Artificial viscosity The artificial viscosity model has been introduced in the SPH method for ensure the stability of the numerical solution in the presence of strong shocks [164]. In this model, a dissipative term Φa, similar to the Von Neumann-Richtmyer viscosity, is added into the momentum equation: Γa= −αh¯cabµab ¯ρab ∇aWab,if vab ·rab ≤0 0,otherwise (3.51) with µab: µab =vab ·rab rab (3.52) The parameter αis a tunable coefficient usually taken in the range between 0.1−0.01, and ¯cab and ¯ρab define the average values of the speed of sound and the density. The advantages of this model is that it preserve the conservation of the angular momentum. As demonstrated in [67], the value of αcan be related to the physical kinematic viscosity nu as: ν=αhc0 2(D+ 2) (3.53) 45
where Dis the dimensionality of the problem. 3.5.2 Laminar viscosity The laminar viscosity model start from the SPH discretization of the viscous term for an incompressible flow Γ = µ∇2v.(3.54) Instead of discretizing using an SPH spatial operator which would be very sensitive to the particle disorder [256], an SPH and a finite-difference first derivative operator are mixed [173]in order to obtain the following discretization of the viscous stress tensor [142]: (µ∇2v)a= 4mb νrab ·∇aWab (ρa+ρb)r2 ab vab.(3.55) With respect to the artificial viscosity model, this model derives directly from the discretization of the viscous term in the governing equations, although it doesn’t conserve the angular momentum. 3.5.3 Density diffusion terms One of the drawbacks of the Weakly-Compressible Smoothed Particle Hydrodynamics is the presence of spurious oscillation in the density field at the spatial scale of particles.This issue has been adressed by adding a diffusion term in the continuity equation. In this section, the Density Diffusion Terms used in the present works are presented. δ-SPH model [8] Molteni and Colagrossi [162] have proposed a diffusive term based on the discretization of the Laplacian of the density field in the form : Da= 2hc0δX b ϕab rab ·∇Wab r2 ab Vb,(3.56) 46
where ϕab = (ρa−ρb) and δis a tunable parameter that usually is taken in the range 0.3−0.1. This term is consistent, because it goes to 0 as the smoothing length h→0 and preserves the mass conservation. However, because of the singularity of the Morris formula in the case of truncated support, it diverges in the presence of a free surface. As remedy, in [8] the δ-SPH model has been proposed, suggesting the following expression for ϕab: ϕab = (ρa−ρb)−1 2∇ρL a+∇ρL b·rba, (3.57) with the renormalized density ∇ρL agiven by: ∇ρL a=X b (ρa−ρb)La∇aWab.(3.58) The renormalization matrix is defined La: La="X b (xb−xa)⊗Wab#−1 ,(3.59) and restore the first-order consistency also in the presence of boundaries or free-surface that truncate the support. Neverthless, it require the inversion of a matrix 2x2 in 2-D and 3x3 in 3-D for every particles in the computational domain, so the additional computational cost is not negligible. Green et al. DDT [86] A different approach to stabilize the density field has been proposed in [84] and is used in this work. The basic idea behind this model is to apply a Godunov SPH scheme to the continuity equation and employing a Roe’s approximate solver for the Riemann problem at the interface of each neighboring particle. A limiter function is then applied to values of the density obtained in this fashion, so to preserve the monotonicity of the scheme and avoid oscillations in the density field. The diffusive term can ben written as: Da= 2c0X b ϕab ||xa−xb|| ∂W ∂q Vb.(3.60) 47
ϕab is now defined as: ϕab =Bab (ρa−ρb)−1 2(ξab∇ρC a+ξba∇ρC b)(xa−xb),(3.61) where: Bab =||xa−xb|| h 2ρa ρr a+ρr b .(3.62) With respect to Equation (3.57), the tuning parameter δis now absent and the magnitude of the diffusion introduced in the governing equations is now automatically adjusted by the definition of reconstructed density values: ρr a=ρa+1 2ξ(η)∇ρC ab ·(xb−xa),(3.63) where ∇ρC ab is a corrected SPH approximation of the density gradient and ξ(η) is the limiter function, with ηdefined as: η=∇ρab ·(xa−xb) ρb−ρa .(3.64) Following the indications in [84], the limiter function chosen in this work is the Van Albada [252]: ξ(η) = η2+η η2+ 1.(3.65) It can be seen that, assuming ξab =ξba = 1.0 and noticing that Bab ≈0.5, this formulation is very similar to Equation 3.57. Fourtakas et al. DDT [75] Fourtakas et al. [75] have proposed a new density diffusion model that is able to maintain the hydrostatic solution while avoiding the computation of the renormalized density gradient. The main concept is to compute the Laplacian of the density field by neglecting the hydrostatic contribution. In fact: Da=δhc0X b ψab ·∇aWabVb(3.66) 48
where ψab is equal to: ψab = 2(ρT ba −ρH ab)xab ||xab||,(3.67) where the superscript Tand Hrefers to the total and to the hydrostatic part of the density. The hydrostatic density difference can be obtained by: ρH ab =ρ0 γ sPH ab + 1 CB−1 (3.68) where PH ab is the hydrostatic pressure difference: PH ab =ρ0gzab (3.69) and CB=c2 0ρ0 γ. 3.6 Particle Shifting An important issue of the SPH method is the formation of void regions, especially in the presence of strong vortex structures, that can affect the stability and accuracy of the numerical solution. This problem has been addressed in [134], where the Particle Shifting Technique (PST) has been proposed: the basic concept is to shift the position of the particles to ensure a more uniform distribution by modelling the shifting δxapplied to the particle position with a Fickian law based on the gradient of the concentration C, in the form: δx=−D∇C(3.70) where Dis the diffusion coefficient defined as: D=Ah2.(3.71) 49
Here, Ais a dimensionless constant that is tuned based on the particular problem, and his the smoothing length. The gradient of concentration ∇Cis calculated using an SPH discretization of the gradient given by: ∇C=X b mb ρb∇aWab.(3.72) In the presence of a free surface, the SPH discretization in Equation (3.72) is inaccurate due to the truncated support, resulting in an upward motion of the particles from the bulk of the flow. To address this issue, in [134], only the tangential component is retained at the free surface while the particle shifting is neglected in the normal direction: δx= 0∇·r≤AF SM , (¯ ¯ I−¯ n⊗¯ n)δx AF ST ≤ ∇·r≤AF SM , δx ∇·r≥AF ST , (3.73) In the previous expression, ¯ ¯ Iis the second order identity tensor and ¯ nis the normal at the free surfaces estimated as: ¯ n=−∇C ||∇C||.(3.74) For detecting the free-surface interface, the method proposed by [125], based on the particle position divergence, is used: ∇·ra=X b Vbrab ·∇aWab (3.75) with AFSM and AFST equal respectively to 1.1 and 1.7 in 2-D and 2.1 and 2.8 in 3-D. 50
ALE corrective terms The shifting formulation with the ALE corrective terms is given by [237]: dρa dt =X bρa ρb mb(vab +δvab)·∇Wab +mb ρb (ρaδva+ρbδvb)·∇Wab−Da,(3.76) dva dt =X bhmbPb+Pa ρaρb∇Wab + (va⊗δva+vb⊗δvb)·∇Wab+ +va(δva−δvb)·∇Wabi+ Γa, (3.77) dxa dt =va+δva.(3.78) 3.7 Solid Boundary Treatment In this work are used the modified Dynamic Boundary Condition (mDBC), proposed in [66], for the definition of solid boundary conditions. This approach is a modification of the Dynamic Boundary Condition (DBC) to address the issue of the unphysical gap between the dummy particles and the fluid particles. In the DBC formulation, a set of dummy particles that represent the solid interface are placed in the computational domain. Their velocity is set to zero, while their density is evolved by using the SPH continuity equation. In this way, when fluid particles come closer to the solid particles, the density and the pressure of the former increase, generating a repulsive force over the fluid particles. However, as discussed in [60], this repulsive mechanism creates a gap in the order of the smoothing length so that there is an incorrect definition of the solid interface that either must be considered prior to placing the solid boundary particles, or after when the in the post-processing of the results. In the mDBC, as in the Ghost Particle approach [149], for each boundary particle, a ghost node is created by mirroring the position of the solid boundary along the solid-fluid interface (Figure 3.2). Once the position of the ghost node is defined, the density and its gradient are computed at the ghost node adoption a corrected SPH operator [140]. This operator consist in the solution opf the following 51
Figure 3.2 Mirroring of ghost nodes (crosses) and the kernel radius around the ghost nodes for boundary particles in a flat surface (a) and a corner (c). Fluid particles (pink) included in the kernel sum around ghost nodes for boundary particles in a flat surface (b) and a corner (d) [66]. 52
linear system Af =b bm=X b fbwmVb Amn =X b wmrnVb wm=Wgb Wx gb Wy gb Wz gb rn=1xgb ygb zgb f=ρgρx gρy gρz g (3.79) where ρgand ρi gare the density and its gradient at the ghost node g. Then, the density at the boundary particle is found through: ρb=ρg+ (rb−rg)·ρgρx gρy gρz g(3.80) In case the matrix Ais ill-conditioned, due to an insufficient number of particles, instead of Equation (3.79), the density is calculated with a Shepard correction: ρg=PbρgWgbVb PbWgbVb (3.81) For the velocity ub, are enforce no-slip condition by an symmetric reflection along the solid interface: ub= 2us−ug(3.82) where usis the velocity of the solid interface, and ugis the velocity at the ghost node calculated as: ug=PbugWgbVb PbWgbVb .(3.83) 53
3.8 Time Stepping Scheme In this work presented later, a second-order symplectic predictor-corrector [127] is employed as time scheme: ρn+1 2 a=ρn a+∆t 2Mn a, vn+1 2 a=vn+∆t 2Fn a, xn+1 2 a=xn+∆t 2vn a, ρn+1 a=ρn a2−ϵ 2 + ϵ;ϵ=−∆tMn a ρa1 2 , vn+1 a=vn+∆t 2Fn+1 2 a, xn+1 a=xn+∆t 2(vn a+vn+1 a) + δx, (3.84) where δx is the particle shifting obtained with Equation (3.70). The time-step is chosen according to the following CFL condition: ∆t=CCFL min(∆tf,∆tcv),(3.85) where: ∆tf= min arh Fa ,(3.86) ∆tcv = min a h c0+ maxb|h(vb−va)·(xb−xa) (xb−xa)2|.(3.87) 54
initial particle distribution generated by a preliminary simulation of the same test case. This approach is similar to the one proposed by Colagrossi et al. [44]. This initialization, shown in Figure 4.3(b), allows to bypass the aforementioned problems and to obtain a much smoother velocity field [Figure 4.3(d),4.3(f)]. Specifically, when using a random-like initial distribution over a Cartesian one, a smaller dissipation is observed at the beginning of the simulation as shown in Figure 4.4(a). This difference between the two initializations seems to lose significance with increasing values of the kernel support. However, using a larger kernel support leads to a larger smoothing error for the SPH spatial operators, therefore this could be detrimental in the latter stages when turbulent structures become much smaller and are at risk of being smoothed out artificially. To verify this and to ultimately choose an optimal initial setup, a sensitivity analysis is performed to assess the effect of the smoothing-length-to-particle-spacing ratio, h/dp, on the numerical solution. The results are depicted in Figure 4.4(b) and indicate that values of h/dp > 1.5 produce worse results for a given dp, especially for t > 5 when the initial vortices are broken into smaller structures. For this reason, a value of h/dp = 1.5 is chosen and kept fixed in all simulations hereafter. Moreover, to satisfy the weakly compressibility assumption, the speed of sound is chosen to be c0= 10 ·vmax, limiting the density variations to 1% of the reference density. 4.3.2 Eulerian SPH Figure 4.5(a) shows the history of the kinetic energy for Re = 1600 with the Eulerian SPH model deprived of any diffusive terms in the continuity equation, i.e., Equation (3.49), and for 643, 1283, 2563and 5123computational nodes. It is remarked that for this Reynolds number, a rough estimation of the ratio between the Kolmogorov length [Equation (2.3)] scale and the reference length gives L/η ≈250, indicating that the two finest resolutions adopted in this study are able to resolve the smallest 61
(a) (b) Figure 4.4 (a) Effect of the initial particle distribution on the decay of the kinetic energy for the 3-D Taylor-Green Vortex at Re = 1600 with two values of the smoothing-length-to-particle-spacing ratio, i.e., h/dp = 1.5 and 2.0. (b) Effect of the smoothing-length-to-particle-spacing ratio on the decay of the kinetic energy for the 3-D Taylor-Green Vortex at Re = 1600 with a random-like particle initial distribution. scales of motion. Results are in good agreement with the high-order pseudo-spectral solution by Van Rees et al. (2011) [253], especially for the finest resolution, for which the grid independence is almost achieved. In Figure 4.5(b), the time history of the enstrophy is reported for the same set of simulations. It is possible to observe that the numerical solution is converging towards the reference solution with the increase of the resolution. 4.3.3 Lagrangian SPH The same test case as in the previous section is simulated here using Lagrangian SPH without any diffusive terms added to the continuity equation. Looking at the kinetic energy decay in Figure 4.6(a), the solution is also converging, however, an excessive dissipation is observed when the process of vortex roll-up starts, i.e., for 3≤t≤6. A comparison with the previous Eulerian SPH results suggests that this over-diffusion could be due to the disorder in the SPH particle distribution, which leads to low accuracy of the SPH spatial operators. This over-diffusion causes the 62
(a) (b) Figure 4.5 (a) Time history of the kinetic energy and (b) time history of the enstrophy for the 3-D Taylor-Green Vortex at Re = 1600 simulated with Eulerian SPH using 643, 1283, 2563and 5123particles. SPH results are compared to the reference solution in Van Rees et al. (2011) [253]. early breakdown of the initial vortices, and a discrepancy with the reference solution is thus observed [Figure 4.6(a)]. As can be seen in Figure 4.6(b), this also causes a poor agreement between SPH and the DNS solution in Van Rees et al. (2011) [253] for what concerns the time evolution of the enstrophy. When comparing the Eulerian and Lagrangian results from Figures 4.5(b) and 4.6(b), respectively, the latter severely underestimates the peak value, remarking the lower order of accuracy of Lagrangian SPH due to the poorer particle distribution. 4.3.4 Influence of the dissipation model Additional simulations of the same 3-D Taylor-Green Vortex case are presented here with the adoption of the density diffusion term (DDT) by Green et al. (2019) [84] in the Lagrangian SPH simulations. In Figure 4.7(a), the solutions with and without the DDT are compared for 5123particles in the computational domain, revealing a negligible effect of the diffusion term from the standpoint of the time history of the flow kinetic energy. However, when looking at the time evolution of the enstrophy [Figure 4.7(b)], the addition of the Green et al. (2019) [84] DDT in the continuity 63
(a) (b) Figure 4.6 (a) History of the kinetic energy and (b) of the enstrophy for the TaylorGreen Vortex at Re = 1600 for resolutions of 643, 1283, 2563and 5123particles with the Lagrangian SPH formulation. Numerical results are compared to the reference solution in Van Rees et al. (2011) [253]. equation is able to yield a higher peak value with respect to the standard Lagrangian SPH, likely due to the less noisy density field which leads to an improvement of the accuracy for this test case. This is confirmed when looking at the contours of the density field in Figure 4.8(a), where it can be clearly seen that when using a Lagrangian SPH model without any dissipation, the density field appears dominated by strong numerical noise. Instead, when the Green et al. (2019) [84] DDT is enabled [Figure 4.8(b)], a much smoother density field is obtained. This effect can also explain the differences between the turbulent energy spectra displayed in Figure 4.9, where for wave numbers in the range of the kernel radius size, the Lagrangian SPH without any DDT shows a spurious flattening, which is instead significantly reduced when the DDT is enabled. Moreover, the Eulerian SPH is capable of computing the spectrum correctly across the whole range of frequencies. 64
(a) (b) Figure 4.7 (a) Time history of the kinetic energy and (b) time history of the enstrophy for the 3-D Taylor-Green Vortex at Re = 1600 simulated with 5123particles and with a Lagrangian SPH formulation with and without the Green et al. (2019) [84] diffusive term. Numerical results are compared against the reference solution in Van Rees et al. (2011) [253]. (a) (b) Figure 4.8 Density field at t= 7 for 3-D Taylor-Green Vortex at Re = 1600: (a) Lagrangian SPH with no dissipation; (b) Lagrangian SPH with the Green et al. (2019) [84] diffusive term enabled. 65
Figure 4.9 Turbulent energy spectra for the 3-D Taylor-Green Vortex at Re = 1600 with 5123particles in the domain and with Eulerian SPH, and Lagrangian SPH with and without the Green et al. (2019) [84] dissipation term. Numerical results are compared against the reference solution in Van Rees et al. (2011) [253]. 66
4.4 Forced Isotropic Turbulence The performance of the above models is evaluated here for an isotropic turbulent flow in a triple periodic box with a linear forcing term added to the momentum equation: dv dt =−∇P ρ+ν∇2v+Av,(4.7) where Ais the forcing constant. This particular test case is chosen to assess the stability and robustness of the above schemes when simulating long periods of physical time, which has been possible to solve mainly thanks to the highly parallel performance of the DualSPHysics code [62]. When compared to a more traditional band-limited forcing, a linear forcing injects energy at all scales of motion but is capable of achieving a stationary state with an energy spectrum as good as band-limited forcing [212]. A solenoidal velocity field peaking at a wavenumber of k0= 2 has been chosen as the initial condition [212]. As discussed by Rosales and Meneveau (2005) [212], the initial velocity field has no influence over the stationary solution which depends only on the Aconstant of Equation (4.7) and on the size do the domain. As for the previous test case, simulations have been performed with both the Lagrangian and Eulerian SPH schemes with and without the density diffusion term of Equation (3.60). A rough estimation of the kolmogorov scales gives L/η ≈200, so a resolution with 2563is chosen. A list of the different simulations that have been run together with details about the SPH formulation, density diffusion scheme and value of the coefficient A used are summarized in Table 4.1. A resolution of 2563particles is used for all cases and the fluid kinematic viscosity is set to ν= 9.9471 ×10−5 4.4.1 Role of the density diffusion term Upon tackling the forced isotropic turbulence problem with an Eulerian SPH scheme without any density diffusion term, simulation results have shown to be unstable due to the insurgence of strong oscillations in the density field [Figure 4.10(a)]. This 67
Table 4.1 Simulation parameters for the Forced Isotropic Turbulence Problem Case SPH Formulation DDT A 1a Eulerian Green et al. (2019)[84] 0.1 2a Lagrangian Green et al. (2019)[84] 0.1 3a Eulerian Green et al. (2019)[84] 0.3 4a Lagrangian Green et al. (2019)[84] 0.3 4b Lagrangian Green et al. (2019)[84] + Shift. ALE 0.3 4c Lagrangian Green et al. (2019)[84] + Shift. 0.3 numerical noise was not detected in the previous TGV case, and it is probably due to the longer simulation times which bring to light the stability issues of a centered-collocated numerical scheme. To overcome this issue, two different dissipation models have been included in the continuity equation and investigated: the δ-SPH by Antuono et al. (2010)[8] and the Green et al. (2019) [84] models. The contours of the density field obtained with the δ-SPH model[8] with δ= 0.1 are depicted in Figure 4.10(b): a higher level of noise can be seen when compared to the solution obtained with the Green et al. (2019)[84] model [Figure 4.10(c)]. This is also highlighted by the probability distributions of density values in Figure 4.11(a), while the density energy spectra for these three solutions in Figure 4.11(b) indicate that δ-SPH is unable to reproduce the low-mid range scales. Although this issue could be partially mitigated by increasing the value of the δparameter according to the magnitude of the forcing term, the Green et al. (2019) [84] model appears to be the better choice due to its ability to automatically adjust the amount of introduced dissipation based on local flow conditions. For this reason, this model is included in the continuity equation for all Eulerian and Lagrangian SPH simulations henceforth. 68
(a) (b) (c) Figure 4.10 Contours of the density field at t= 9 for the forced isotropic turbulent cases with A= 0.1: (a) Eulerian SPH; (b) Eulerian SPH with the Antuono et al. (2010)[8] dissipation model; (c) Eulerian SPH with the Green et al. (2019)[84] dissipation model. (a) (b) Figure 4.11 Forced isotropic turbulent problem: (a) Probability distribution of the density field, and (b) energy spectra at t= 9 for Eulerian SPH, Lagrangian SPH with the Antuono et al. (2010)[8] dissipation model and Lagrangian SPH with the Green et al. (2019)[84] dissipation model. 69
(a) (b) Figure 4.12 (a) Time history of v2 rms, and (b) time history of ϵ/3v2 rms for Eulerian SPH results of the forced isotropic turbulent problem (case 1a). 4.4.2 Differences between Eulerian and Lagrangian approaches Two global quantities are defined to investigate the stationary behavior of the problem: the mean square value of the velocity fluctuations, v2 rms =⟨v·v⟩/3, and the eddy turnover time, τ= v2 rms/ϵ, where ϵ=−ν⟨v·∇2v⟩is the mean dissipation rate. As can be seen in Figure 4.12(a), v2 rms reaches a stationary value for t/τ > 5 up to 20 eddy turnover times in all the different cases. From the kinetic energy balance of a steady state problem, it is known that the ratio ϵ/3v2 rms must lie around the value of the parameter Aspecified in the governing equations. This is corroborated by the time history of ϵ/3v2 rms in Figure 4.12(b), where the correct value of A= 1 is retrieved. In Figures 4.13(a) and 4.13(b), the results for the same set of simulations are shown when Lagrangian SPH is employed. It can be seen that v2 rms is slightly less oscillatory for the Lagrangian scheme. This notwithstanding, the differences between the two approaches are not significant and the final values of the mean square value of the velocity fluctuations are almost identical when the stationary solution is reached. As with the Eulerian SPH, also Lagrangian SPH predicts the correct value of ϵ/3v2 rms well, with a stable result for t/τ > 8 Figure 4.13(b)). 70
(a) (b) (c) Figure 4.17 History of the kinetic energy for the (a) 2nd, (b) 4th and (c) 6th order schemes for the Taylor-Green Vortex at Re=1,600. Numerical results are compared to the reference solution in [253]. 77
(a) (b) (c) (d) (e) (f) Figure 4.18 Time evolution of the enstrophy and kinetic energy dissipation rate for the (a)-(b) 2nd, (c)-(d) 4th and (e)-(f) 6th order schemes for the Taylor-Green Vortex at Re=1,600. Numerical results are compared to the reference solution in [253]. 78
(a) (b) (c) (d) Figure 4.19 Vorticity contours for ω= 1,5,10,20,30 at x=−0.5 for the (a) 2nd, (b) 4th and (c) 6th order schemes for the Taylor-Green Vortex at Re=1,600. (d) Reference solution in [253] . 79
Figure 4.20 Turbulent energy spectrum 2nd, 4th and ) 6th order schemes for the Taylor-Green Vortex at Re=1,600 with N= 2563particles. Numerical results are compared to the reference solution in [253]. 80
CHAPTER 5 MULTI-RESOLUTION ALGORITHM In the present chapter, an adaptive resolution algorithm for SPH is presented. The scheme presented herein is based on the decomposition of the computational domain into different sub-domains, each with its own characteristic particle size and smoothing length. At each time step, the computational problem is solved in every sub-domain independently, and the sub-domain closure is provided by a buffer region that acts as a Dirichlet boundary condition. Coupling between the different sub-domains is obtained by interpolating the physical quantities in the buffer region using information available at the fluid particles lying in adjacent domains. 5.1 Multi-Resolution Algorithm The main idea behind the variable resolution algorithm of this study is based on a decomposition of the computational domain, Γ, into a set of Nsub-domains, Γiwith i= 1, . . . , N , such that Γ = SN i=1 Γi[Figure 5.1(a)]. Each sub-domain is characterized by its own characteristic particle size, dpi, and smoothing length, hi. To solve the computational problem from tnto tn+1, a closure Dirichlet boundary condition must be provided at the boundaries of each sub-domain. To this end, each sub-domain Γiis extended by appending a buffer region, ∂Γi, with width l∂Γi= 2hiand such that ∂Γi∈Γjand ∂Γj∈Γi. The latter condition allows establishing a bijective topological connection between two sub-domains, formalized with the notation ∂Γj i, i.e., the buffer region of sub-domain iis coupled with the sub-domain j. Every buffer region is populated by special SPH particles following a strategy similar to the one adopted in [238] for open boundary conditions. The “buffer” particles are distinct from regular fluid particles as they are not updated using the governing equations of motion. Instead, at the beginning of each time step, a buffer 81
(a) (b) Figure 5.1 (a) Example of two sub-domains Γ1and Γ2. (b) Buffer regions ∂Γ2 1 and ∂Γ1 2with widths l∂Γ2 1= 2h1and l∂Γ1 2= 2h2are appended to their respective sub-domains. particle a∈∂Γj iobtains its physical properties from a spatial interpolation over fluid particles in Γj, as shown in Figure 5.2. Buffer particles located close to the interface have a truncated support, therefore the corrected SPH interpolation proposed in [139] is employed to ensure consistency. The procedure for restoring particle consistency proposed in [139] starts from a multi-dimensional Taylor series expansion of a field function f(x) multiplied by the kernel function and its derivative up to the desired order of consistency. In this work, a second-order consistency condition is enforced, which in two dimensions results in 82
Figure 5.2 Coupling procedure between sub-domains Γ1and Γ2. Buffer particles (orange) interpolate their properties over the fluid particles (light blue) in the coupled subdomain. the following linear system for the generic buffer particle a∈∂Γj i: Af =b(5.1) bm=X b∈Γj fbwmVb(5.2) Amn =X b∈Γj wmrnVb(5.3) w=Wab Wx ab Wy ab Wxx ab Wxy ab Wyy ab (5.4) r=1xab zab x2 ab 2xabyab y2 ab 2(5.5) f=fafx afy afxx afxy afyy a(5.6) where Wis the kernel smoothing function, xab = [(xa−xb),(ya−yb)] is the interparticle distance, and Vis the volume. Once the properties are obtained through the 83
Figure 5.3 A buffer particle that moves into the fluid domain is transformed into a fluid particle (process 1). A fluid particle that enters the buffer region is transformed into a buffer particle (process 2). A buffer particle that moves outside the extended subdomain ∂Γj i∪Γiis deleted (process 3). interpolation procedure, the position of the buffer particle ais updated in time using a simple Euler time integration scheme: xn+1 a= ∆tvn a.(5.7) where xn+1 ais the position of the particles at time-step tn+1, and vn ais the velocity at tn. To handle the exchange of mass among the different sub-domains, the following conditions are checked at the beginning of each time step: 1. If a buffer particle moves into the fluid domain, it is converted to a fluid particle (Figure 5.3, process 1). 2. If a fluid particle enters the buffer region, it becomes a buffer particle (Figure 5.3, process 2). 3. If a buffer particle moves outside the extended subdomain ∂Γj i∪Γi, this particle is deleted (Figure 5.3, process 3). 84
Figure 5.4 depicts the procedure of particle insertion into a given sub-domain. The external boundary of the generic sub-domain Γiis divided into mass segments of length dpiand the Eulerian mass flux across the segments is calculated at each time step using the mid-point rule: ˙ma= max(0,−ρma(vma−vbf)·ndpidt) (5.8) where ˙mais the mass flux during the time step dt,ρmaand vmaare the density and the velocity calculated at the mass accumulation point using the same corrected interpolation employed for buffer particles, and vbfis the velocity of the interface in case of a moving subdomain. After ˙mais added to ma, if ma≥(dp)D/ρ0, where Dis the number of problem dimensions, a particle is inserted in the buffer region at a distance dp/2 in the normal direction to the interface, as shown in Figure 5.4. In presence of free-surface, this procedure has to be adjusted to prevent a non-zero mass flux at the accumulation point above the water level. The following formulation based on the free-surface detection method proposed in [125] is employed for the computation of the mass-flux at the interface: ˙ma= max(0,−ρma(vma−vbf)·ndpidt) if ∇·r≥ ∇·rth 0 if ∇·r<∇·rth; (5.9) where ∇·r=Pb mb ρb·∇Wab and ∇·rth is a threshold value equal to 1.5 in 2-D and 2.75 in 3-D. Therefore, because buffer particles move with a Lagrangian velocity which is not obtained from the governing equations directly but from an interpolation over fluid particles in adjacent sub-domains, they do not benefit from the regularization effect of the pressure gradient discretization [198]. These issues can yield an irregular particle distribution in the buffer regions, which eventually affects the correct enforcement of boundary conditions for the sub-domains. A solution to this problem is found 85
Figure 5.4 Particle insertion procedure: at each time step, the normal mass flux at the outer boundary of the sub-domain is calculated and added to the mass accumulation points (red squares). When the mass at the accumulation points reaches the reference particle mass, new particles (green) are created. by using the Particle Shifting Technique (PST) in the buffer regions, with special attention to avoid shifting the buffer particles toward the edges of the sub-domain. The latter risk is mitigated by surrounding the edges of the buffer regions with layers of “fixed particles” which have a fixed position in space and constant density ρ0. Other than these two physical variables, fixed particles do not have any other physical quantity associated and they interact with buffer particles only with the sole purpose of computing the shifting correction. In contrast to the shifting formulation adopted for fluid particles, buffer particles are shifted only in the direction tangential to the interface between the buffer and the fluid region to avoid inaccuracies when computing the mass flux at the interface. Moreover, because normal vectors are ill-defined in the corner regions, the shifting is disabled in these areas (Figure 5.5). Algorithm 1 shows a pseudo-code of the proposed multi-resolution strategy including all the salient steps of the simulation. 86
Table 5.3 Main Loop Functions and their Description Function Description Predictor step. Interaction Forces Call for particle interaction (PI). PreInteraction Forces Prepares variables and assigns memory for PI. Interaction Forces Computes particle interaction. DtVariable Computes the value of the new variable time step. ComputeSymplecticPre Computes System Update using PosInteraction Forces Memory release of arrays in GPU. RunCellDivide Generates neighbour list. CORRECTOR Corrector step. Interaction Forces Call for particle interaction. PreInteraction Forces Prepares variables and assigns memory for PI. Interaction Forces Computes particle interaction. DtVariable Computes the value of the new variable time step. ComputeSymplecticCorr Computes system update using symplectic=corrector. PosInteraction Forces Memory release of arrays in GPU. RunCellDivide Generates Neighbour List. FinishRun Shows and stores final overview of execution. of different research groups, the DualSPHysics package is constantly updated with new functionalities. To keep the range of applicability of DualSPHysics intact and possibly to extend it, the variable resolution algorithm was implemented, favoring the preservation of 93
the base structure of the code, avoiding major changes that could lead to long-term maintenance issues. There are two observations that can be made by looking at the algorithm proposed in Section 5.1 for variable resolution and to the DualSPHysics code structure outlined previously: •Each computational sub-problem is dependent on the other ones only when the physical properties of the buffer particles are interpolated •To each computational sub-problem corresponds a JSphSingle class instance. Considering these two aspects, the implementation of the new algorithm has been structured by operating at the highest level of abstraction. A new class, named JSphGpuWrapper, has been implemented: in this new object, an array of JSphGpuSingle objects is allocated, and at each object corresponds a numerical subdomain. In this way, each computational sub-problem can be managed independently and synchronized to the other when needed. •The number of sub-domains is read from the configuration file •The JSphGpuSingle array is allocated •Each JSphGpuSingle is initialized and configurated •The main loop for advancing the simulation is executed. •Subroutine for completing the simulation. The new main loop is structured in the same way as in the standard DualSPHysics implementation. Still, at the end of the time-step, the routines for updating the buffer regions are executed, following the pseudocode outlined in Algorithm 1: 1. For each subdomain, the solution is advanced in time by calling the same subroutines as in the original DualSPHysics code. 94
Table 5.4 Main loop Functions and their Description Function Description Interaction BufferExtrapFlux Call for computing the flux at the accumulation point. ComputeStep Procedure for transformation between fluid and buffer particles. BufferListCreate Definition of the point in which buffer particles must be created. BufferCreateNewPart Creation of buffer particles. RunCellDivide Generates neighbour list. Interaction BufferExtrapFlux Call for interpolating the buffer particles. 2. Each sub-domain computes the coupling part of the new multi-resolution algorithm. A visualization of the structure of the main loop for the multi-resolution implementation is shown in Figure 5.7, while in Table 5.4 are described the functions of the coupling section of the multi-resolution algorithm is divided into these main steps, which can be divided in three main steps: 1. Computation of the flux at the accumulation point. 2. Evaluation of the mass segment and particle creation and deletion procedure. 3. Particle reordering and interpolation for obtaining the physical properties of the buffer particles. The main advantage of the present implementation strategy is that the main structure of the DualSPHysics code is left unchanged as the changes introduced for the variable 95
Figure 5.7 Call function for the main loop in the new multi-resolution algorithm. 96
Table 5.5 New Files Added to the Source Code Files Description JSphGpuWrapper.cpp/.h Implements the class JSphGpuWrapper. JSphGpuSingleBuffer.cpp Define the host function for the coupling algorithm. JSph Gpu Buffer.cu Define the Cuda kernel function for the coupling algorithm. JSphBuffer.cpp/.h Implements the class JSphBuffer. JSphBufferZone.cpp/.h Implements the class JSphBufferZone. Interaction BufferExtrapFlux Call for interpolating the buffer particles. resolution algorithm are additive. In Table 5.5 are listed the new files introduced in the source code. Two new classes are introduced into the code: JSphBufferZone In this class, the geometrical definition of the buffer region is defined. JSphBuffer Define an array of JSphBufferZone objects, one for each buffer region of the sub-domain. It also implements the routine needed to retrieve the list of particles in the buffer region and label them as buffer particles. 5.3 Results and Discussion 5.3.1 Hydrostatic tank Despite being a simple test case, the hydrostatic tank is often very difficult for the SPH method, because of the inconsistency due to the presence of a free-surface and the tank wall that introduce numerical error in the solution. Moreover, is an ideal test case to verify the accuracy of the coupling between sub-domain and any inaccuracies introduced by the interpolation of buffer regions. 97
Figure 5.8 Computational domain for the hydrostatic tank case For these reasons, an hydrostatic tank is chosen as a study case to validate the proposed variable-resolution algorithm. The computational domain consists of a 2-D rectangular tank, with width L= 1m, filled with water up to a height H= 1m. The water density is set equal to ρ0=ρ∞= 1000 kg/m3and the gravity g= 9.81m/s2. Regarding the numerical parameters, an artificial viscosity model is used with α= 0.01, while the speed of sound is equal to c0= 10vmax, where vmax =√gH. The Fourtakas DDT in Equation 3.66 is used to order to preserve the hydrostatic solution. Two refinement zone are employed, the coarse one with a resolution equal to H/∆x= 50, and the finest one with H/∆x= 100. The refinement region is a square box with length h1= 0.5m, centered at the midline of the rectangular tank. In Figure 5.8 the geometrical definition of the computational domain and of the refinement regions is shown. 98
(a) (b) Figure 5.9 Hydrostatic tank case: (a) density contours and (b) pressure distribution against the hydrostatic solution at t= 20s Figure 5.9(b) reports the pressure distribution against the vertical coordinate at t= 20s: as can be seen, the present multi-resolution approach is able to preserve the hydrostatic solution, and no discontinuities are visible at y/H = 0.5, where the horizontal interface between the coarse and the fine resolution is placed. Moreover, despite a low value of the artificial viscosity coefficient α, in Figure 5.9(a) the density contours at t= 20sreveal a regular distribution of the particles, that doesn’t show any spurious motion, especially at the interface between the sub-domains. 5.3.2 Flow past a fixed, circular cylinder Flow past a circular cylinder is simulated for several Reynolds numbers (Re) as a first study case to benchmark the new multi-resolution algorithm against consolidated literature results. This problem has been investigated extensively both numerically and experimentally, thus it is suitable for demonstrating the advantages of the present multi-resolution approach with respect to using a uniform SPH particle resolution. Since the Reynolds number based on the smoothing length, Reh=U∞hν−1, must be in the order of 1 to resolve the flow near the cylinder properly, this test case is 99
Figure 5.10 Computational domain for 2-D flow past a circular cylinder quite challenging to resolve with a uniform particle resolution, even at low to moderate Reynolds numbers, due to the prohibitive number of particles required for the domain discretization. Moreover, the level of resolution needed in the far-field is substantially lower than around the cylinder and thus the proposed multi-resolution algorithm is adopted to reduce the average particle resolution in the far-field without affecting the global accuracy of the simulation. The computational domain shown in Figure 5.10 is chosen following the work in [238], where a cylinder with diameter D= 0.1 m is centered at (x, y) = (0,0) and the overall domain dimensions are set to 25D×20Dto minimize the change of any blockage effect. A no-slip solid boundary condition is applied to the cylinder, while free-slip conditions are defined at the bottom and upper walls. Inflow-outflow conditions are imposed using the formulation in [238]: at the inlet, the velocity and 100
(a) Re=100 (b) Re=200 Figure 5.11 Dimensionless pressure for flow past a cylinder with (a) Re = 100 and (b) Re = 200 the density are prescribed, whereas at the outlet, the velocity is extrapolated from the fluid to the buffer region, while the density is prescribed to the reference value. The fluid is initialized with U∞(x, y) = (1,0) m/s while the density has an initial value equal to ρ0=ρ∞= 1000 kg/m3. The Reynolds number Re=U∞Dν−1is varied by changing the value of the kinematic viscosity ν. The smoothing length is set to h= 2∆x, constant across the different sub-domains, and the DDT in Equation (3.66) is activated in order to diffuse the oscillations affecting the density field due to the centered collocated SPH scheme. Finally, PST is applied to avoid the creation of empty regions due to the vortical structures that develop in the wake. For the setup of the different resolution sub-domains, a minimum resolution corresponding to D/∆x= 25 is chosen for the far field. Then, new sub-domains are nested inside the far-field, each with a resolution that doubles as they get closer to the cylinder. The number of subdomains created depends on the flow Reynolds number, with higher Re requiring higher overall resolution, therefore more sub-domains. Table 5.6 summarizes the number of sub-domains used and their respective dimensions and particle resolutions as a function of the Reynolds number. 101
Table 5.6 Number of Sub-Domains Used and their Respective Dimensions and Particle Resolution as a Function of the Reynolds Number Case Re Number of zones D/∆xmax 1 100 3 100 2 200 3 100 3 1000 4 200 5 400 6 800 4 3000 5 400 6 800 7 1600 5 9500 6 800 7 1600 8 3200 Re = 100 and 200 For the first two Reynolds numbers Re = 100,200, only three sub-domains are utilized as shown previously, and these refinement regions are centered at the cylinder and extended downstream to better resolve the wake region. The contours of dimensionless pressure P⋆=P(x/D, y/D)ρ−1U2 ∞shown in Figure 5.11 depict the vortices’ cores, clearly visible as low-pressure regions developing in the wake. These vortical structures have higher intensity with increasing Reynolds number. No discontinuities are visible through the interface between the fine and the medium refinement regions, demonstrating the robustness of the coupling procedure between sub-domains. Furthermore, the dimensionless vorticity ζ⋆=ζ(x/D, y/D)DU−1can be observed in Figure 5.12, colored such that clockwise and counter-clockwise fluid rotations tend towards blue and red, respectively. The typical von Karman street associated with the periodic vortex shedding is captured correctly and the vortical structures cross the interface of the different resolution 102
(a) t∗= 1 (b) t∗= 2 (c) t∗= 3 (d) t∗= 4 (e) t∗= 5 (f) t∗= 6 Figure 5.18 Vorticity contours for flow past a cylinder at Re = 1000 109
(a) t∗= 1 (b) t∗= 2 (c) t∗= 3 (d) t∗= 4 (e) t∗= 5 (f) t∗= 6 Figure 5.19 Streamlines for flow past a cylinder at Re = 1000 110
(a) t∗= 1 (b) t∗= 2 (c) t∗= 3 (d) t∗= 4 (e) t∗= 5 (f) t∗= 6 Figure 5.20 Vorticity contours for flow past a cylinder at Re = 3000 111
(a) t∗= 1 (b) t∗= 2 (c) t∗= 3 (d) t∗= 4 (e) t∗= 5 (f) t∗= 6 Figure 5.21 Streamlines for flow past a cylinder at Re = 3000 112
(a) t∗= 1 (b) t∗= 2 (c) t∗= 3 (d) t∗= 4 (e) t∗= 5 (f) t∗= 6 Figure 5.22 Vorticity contours for flow past a cylinder at Re = 9500 Around t∗≈2 [Figures 5.22(b) and 5.23(b)], the primary vortex pair detaches from the cylinder, as can be seen also by looking at the time history of the drag coefficient 113
(a) t∗= 1 (b) t∗= 2 (c) t∗= 3 (d) t∗= 4 (e) t∗= 5 (f) t∗= 6 Figure 5.23 Streamlines for flow past a cylinder at Re = 9500 114
in Figure 5.17(c). Meanwhile, a new vortex pair is formed at a higher angle, and the same interplay between primary and secondary vorticity is observed, increasing the pressure drag up to t∗= 3. Interestingly, a t∗= 4, the second vortex pair merges with the previous, primary vortices after being separated from the body. The evolution at later times is dominated again by the complex interaction among vortices in the boundary layer, resulting in a highly unsteady flow. The convergence study performed for this last Reynolds number has resulted in using 6, 7, and 8 different resolution zones, reaching a maximum resolution near the cylinder equal to D/∆x= 800,1600,3200, respectively. The quantitative effect of these three different resolutions can be assessed in Figure 5.17(c), where the time history of the drag coefficient is shown. While all three resolutions capture the drag coefficient well until t∗= 4, only the simulation with D/∆x= 3200 is a good match with the reference solution for t∗>4, when the flow is characterized by a high degree of unsteadiness. Moreover, with the highest resolution adopted, the flow remains almost perfectly symmetrical as is also the case for the reference solution in [118]. Contrary to what has been shown for Re = 100 and 200 in Section 5.3.2, no comparisons between multi-resolution and uniform resolution SPH simulations are shown here for Re = 1000, 3000 and 9500. This is mainly due to the prohibitive cost of running uniform resolution simulations for the three Reynolds numbers selected when the highest resolution has to be used. For example, the multi-resolution simulation for Re = 9,500 with D/∆x= 3200 requires 8 sub-domains, resulting in a total number of SPH particles deployed of N≈10 ×106. Similarly, a uniform resolution SPH simulation with this level of particle spacing would require N≈9×109particles, which can be achieved only with sophisticated memory-distributed parallelization [63]. 115
5.3.3 Flow past an oscillating, circular cylinder The flow past a circular cylinder oscillating in a transverse direction to the streamwise direction is investigated for Re = 100 for different values of the oscillation amplitude and frequency. This test case generates complicated flows and it has been investigated with different techniques, including experiments computational studies [194] and experimental investigations [32, 33, 268, 5]. Herein, flow past an oscillating cylinder has been chosen to assess the capability of the proposed multi-resolution scheme to simulate flows with moving boundaries. The computational setup is the same as the one employed for studying the flow past a fixed cylinder case (see Section 5.3.2) and the domain has been discretized with three different resolution zones, with particle size equal to D/∆x= 25,50,100 (see Figure 5.10). A sinusoidal motion y(t) is applied to the cylinder in the cross-flow direction: y(t) = ADsin(2πFfst) (5.11) where A=ymax/D, with ymax equal to the maximum displacement and F=f0/fs, where fsis the frequency of the vortex shedding when the cylinder is fixed. The refinement regions also move according to Equation (5.11), so that the relative position of the cylinder with respect to the refinement areas does not change in time. Four different configurations of amplitude and frequency have been considered which can be found in Table 5.8. In Figure 5.24(a), the time history of the lift coefficient CLis shown for [A, F] = [0.25,0.9]. This is defined as a “locked configuration” because the vortex shedding presents a dominant frequency equal to f0, also corroborated by the Power Spectral Density (PSD) of the lift coefficient time signal shown in Figure 5.25(a), where only one peak is visible. For an unlocked case, instead, the lift signal contains two or more frequencies, and this is achieved by changing the frequency ratio to F= 0.5 and F= 1.5. In the 116
Table 5.8 Amplitude and Frequency Values Chosen for the Simulation of Flow Past an Oscillating Cylinder Case A F 1 0.25 0.9 2 0.25 0.5 3 0.25 1.5 4 1.25 1.5 (a) (A, F ) = (0.25,0.9) (b) (A, F ) = (0.25,0.5) (c) (A, F ) = (0.25,1.5) (d) (A, F ) = (1.25,1.5) Figure 5.24 Time history of the lift coefficient for flow past an oscillating cylinder at Re = 100 for different amplitude Aand frequency Fratios: (a) (A, F) = (0.25,0.9), (b) (A, F) = (0.25,0.5), (a) (A, F) = (0.25,1.5), (a) (A, F) = (1.25,1.5) first case, the periodic signal is characterized by two modes with different magnitudes and period tb= 2 t0, where t0is the period of the imposed sinusoidal motion. This is 117
(a) (A, F ) = (0.25,0.9) (b) (A, F ) = (0.25,0.5) (c) (A, F ) = (0.25,1.5) (d) (A, F ) = (1.25,1.5) Figure 5.25 Power Spectral Density (PSD) for the flow past an oscillating cylinder at Re=100 for different amplitude Aand frequency Fratio confirmed by looking at the PSD in Figure 5.25(b), where two peaks associated with f=f0and f= 2f0can be observed in agreement with findings in [194]. The same beating phenomena [194] is also present for the F= 1.5 case, although in this case the lift coefficient has a period tb= 8t0. The last case simulates a more challenging configuration obtained by using (A, F) = (1.25,1.5). Figures 5.24(d) and 5.25(d) highlight a dominant frequency f=f0coupled with a second f= 2f0and third f= 3f0frequency, in agreement with the results reported by [194]. These small secondary frequencies are not associated with beating phenomena, rather they influence the vortex shedding patterns, resulting 118
Figure 6.1 Longitudinaland crosssections of the computational domain for the flow past a sphere. the sphere. At the same time, in the downstream direction, the computational domain is extended by 20Din order to resolve at least 3 vertical structures. A square cross-section with size 10Dis used, which results in a blockage ratio approximately equal to 0.7%. The definition of the boundary conditions follows the setup followed in subection 5.3.2: a no-slip boundary condition is applied to the sphere, while free-slip conditions are defined at the cross-section walls. Inlet-outlet boundary conditions and the initial conditions are imposed again in the same manner as in section 5.3.2. However, different from the two-dimensional study, the value of the smoothing length is h= 1.5∆xin order to reduce the computational cost. The DDT in Equation 3.66 and the PST are applied to the entire domain. For the setup of the refinement regions, a minimum resolution equal to D/∆x= 6.25 is chosen. As for the two-dimensional study on the flow past a cylinder, new sub-domains are nested inside the far field, each with a resolution that doubles as they get closer to the cylinder, with a maximum resolution equal to D/∆x= 100, which corresponds to a ratio between the boundary layer δbl and the particles size equal δbl/∆x≈6−5 for a range Re = 300 −500. It is noted that the total number of particles in these simulations is Np ≈11 ×106, while with a uniform resolution, the number of particles would have been equal to Np= 3 ×109. 125
Table 6.1 Numerical Results of the Drag Coefficient CD, Lift Coefficient CLand Strouhal Number St Values for the Flow Past a Sphere at Re=300 Case CD∆CDCL∆CLSt Present 0.667 0.0032 0.066 0.017 0.134 Costantinescu and Squires [46] 0.655 0.0032 0.065 – 0.136 Johnson and Patel [102] 0.656 0.0035 0.069 0.016 0.138 Tomboulides and Orszag [242] 0.671 0.0028 – – 0.136 Ploumhans et al. [195] 0.683 0.0025 0.061 0.014 0.13 6.1.2 Re=300 Herein, the first case considered is for a Reynold Number Re = 300. The flow is simulated for an overall duration of t∗= 250 time unit, where t∗=tU/D, and only the last 120 time units are considered for collecting the flow statistics. In Figure 6.2(a) the time history of the aerodynamical coefficients is shown: the average value of the drag CD=Fx/(0.5ρU2 ∞πD2/4) and lift CL=Fy/(0.5ρU2 ∞πD2/4) are respectively 0.667 and 0.066 with an oscillation amplitude ∆CDand ∆CL equal to 0.0032 and 0.17. These values agree with other numerical investigations, as seen in Table 6.1, where the present results are summarized and compared against the literature. From Figure 6.2(a) is also evident that the side coefficient CS=Fy/(0.5ρU2 ∞πD2/4) is equal to 0, indicating that the plane (x, y) is a plane of symmetry for the present flow, as expected from previous experimental and numerical investigations with Re = 300 [242, 102]. From the frequency analysis in Figure 6.2(b), where the Power Spectral Densities of CDand CLare reported, it can be found the value of the non-dimensional vortex shedding frequency St = 0.134, which is again in agreement with other numerical investigation, for example with the results in [242] where is predicted at a value St = 0.137. Interestingly, as in [102], the drag coefficient presents a peak in the spectrum at a value equal to 2St, which is not present in the lift coefficient spectrum. 126
(a) (b) Figure 6.2 (a) Time histories and (b) normalized power spectral densities of the aerodynamic coefficients for the flow past a sphere at Re = 300. To further investigate this aspect, Figure 6.3 displays the time history of the stream-wise velocity and its power spectral density at a point distant 5.75Dfrom the sphere in the downstream direction. The only peaks present are at a value equal to St = 0.134, along with its three super-harmonics at 2St, 3St, and 4St, as in [242]. Figure 6.4 presents the average streamwise velocity Uavg along the wake centerline, with a comparative analysis against the numerical outcomes from [242]. A good correlation is observed up to a point x/D = 7, beyond which the displayed results indicate lower values for Uavg. The root-mean-square value of the streamwise velocity on the same centerline is also highlighted in the figure, with the current simulation predicting a diminished peak value at x/D = 1.7. Further downstream there is a reasonable agreement between the two sets of results. The flow visualization displayed in Figures 6.5-6.6, where the Q criterion, a vortex identification method, is shown at every quarter of the vortex shedding period, helps understand the development of the vortical structures in the wake. In Figures 6.6(a)-6.5(a), a vortex is about to separate from the vortical structures surrounding the sphere. This coherent structure is convected downstream, as can be seen in Figures 6.6(b)-6.5(b)), where are now clearly visible the legs of the hairpin vortex. 127
(a) (b) Figure 6.3 (a) Time history and (b) normalized power spectral density of the value of the streamwise coefficient at a x/D = 5.75 on the wake centerline for the flow past a sphere at Re = 300. Figure 6.4 Comparison of the average and root-mean-square values of the streamwise velocity along the wake centerline with the numerical results in [242] 128
While the shed-wake structure moves further away from the sphere [Figures 6.6(c)- 6.5(c)], a new structure, resulting from the interaction between the wake and the outer flow, is formed around the legs of the hairpin vortex (6.6(d)-6.5(d)). The last figure also shows the development of a new coherent structure in the near-wake region. To summarize, the vortex-shedding mechanism is characterized by two hairpin structures oriented in opposite directions; one results from the wake-shed, the other from the interaction between the wake and the outer flow. 6.1.3 Re=500 To further validate the present multi-resolution algorithm, the Reynolds number is raised up to 500, where the flow past a sphere is characterized by the loss of the planar symmetry and a more complex vortex shedding. The numerical simulation is computed for a total of t∗= 300 time unit, and the last t∗= 180 are used to collect the statistics. In Figure 6.7(a), the time history of the force coefficient is reported. The average value of CDof the present simulation is equal to 0.566, which is in perfect agreement with the results in [160]. With respect to the case with Re = 300, the loss of planar symmetry can be discerned from the evolution of the side coefficient. The spectral analysis of the aerodynamical coefficients is reported in Figure 6.7(b), reveals the more complex nature with respect to the case with Re = 300, where only one dominant peak, corresponding with the vortex shedding frequency, was present in the spectrum of CD. Instead, in this case, the power spectral density of the drag coefficient is characterized by different peaks, with the dominant one at St2= 0.44, which is close to the value found in [126, 51]. As for the previous case, to further analyze the flow’s spectral characteristic, in Figure 6.8, the normalized power spectral densities for the stream-wise velocity are shown, and the pressure at three different locations in the near wake region, chosen to match the same measurement 129
(a) (b) (c) (d) Figure 6.5 Flow visualization by the Q criterion, colored by the velocity magnitude U, at every quarter of period from a view normal to the (x, z) plane for the flow past a sphere at Re = 300. 130
(a) (b) (c) (d) Figure 6.6 Flow visualization by the Q criterion, colored by the velocity magnitude U, at every quarter of period from a view normal to the (x, y) plane for the flow past a sphere at Re = 300. 131
(a) (b) Figure 6.7 (a) Time histories and (b) normalized power spectral densities of the aerodynamic coefficients for the flow past a sphere at Re = 500. Table 6.2 Numerical Results of the Drag Coefficient CD, Strouhal Number St1and St2Values for the Flow Past a Sphere at Re=500 Case CDSt1St2 Present 0.566 0.0167 0.044 Mittal and Najjar [160] 0.56 0.150 0.05 Lee [126] 0.54 0.164 0.045 Tomboulides and Orszag [242] - 0.167 0.045 Crivellini et al. [51] 0.558 0.157 0.044 points in [242, 51]. The first point on the wake centerline at (x, y, z) = (2.5D, 0,0) shows the presence of two different frequencies, one at St2= 0.44, the second at St1= 0.166, while at locations 2 and 3, the dominant frequency is the St1= 0.166. The values of the low and high frequencies found are in good agreement with the results in [242], as can be seen from Table 6.2, where the value of the drag coefficient and the Strouhal numbers found in this study are compared against other experimental and numerical investigations. 132
(a) Streamwise velocity (b) Pressure (c) Streamwise velocity (d) Pressure (e) Streamwise velocity (f) Pressure Figure 6.8 Normalized power spectral density for the streamwise velocity and pressure at (a)-(b) (x, y, z) = (2.5D, 0,0), (c)-(d) (x, y, z) = (2D, 0.3,0), (c)-(d) (x, y, z) = (2D, 0.0,0.3), for the flow past a sphere at Re = 500. 133
In Figure 6.9 the iso-surfaces for the Q criterion at every period T1= 1/St1are displayed: it is clear that the St1is the frequency of the vortex shedding mechanism, which resembles the same shedding process seen at Re = 300, although in this case, the vortex structure sequence is more irregular. Moreover, it can be seen that the vortex orientation changes from cycle to cycle. 6.2 3-D Dam-Breaking Flow Herein, the ability of the proposed multi-resolution algorithm to deal with threedimensional free-surface flows is validated by simulating a dam-breaking flow that impacts an obstacle. This test case is very popular in the SPH community and has been chosen as SPHERIC benchmark case #2. The configuration of the numerical experiment [113] is reported in Figure 6.10. A baseline spatial resolution ∆xmax = 0.005 is chosen, and a refinement region has been defined around the obstacle, wherein a finer resolution of ∆xmin = 0.0025 is applied, as shown in Figure 6.11. A uniform value across the different sub-domain for the smoothing length to particle distance is set equal to h/∆x= 1.5. The reference density has a value equal to ρ0= 1000 kg/m3, while the gravitational accelleration g=−9.81 m/s2. The sound speed is chosen equal to c0= 10√gH, where H= 0.55mis the height of the water column. To ensure the numerical stability, an artificial viscosity model with α= 0.01 is used, which corresponds to a physical kinematic viscosity ν= 8.6×1−−5, along with the DDT term in Equation 3.66. No-slip wall boundary conditions are enforced on the walls. Figure 6.12 shows the comparison between the numerical simulation and the experimental results of [113] at the pressure probes placed on the front side of the obstacles. A good agreement is found at gauges P1 and P2, where the numerical solution is able to predict quite well the evolution of the pressure at the first impact. Some discrepancies are instead found at P3 and P4, where the SPH simulation under134