Simulating unsteady conduit flows with smoothed particle hydrodynamics

(1)

Simulating unsteady conduit flows with smoothed particle

hydrodynamics

Citation for published version (APA):

Hou, Q. (2012). Simulating unsteady conduit flows with smoothed particle hydrodynamics. Technische Universiteit Eindhoven. https://doi.org/10.6100/IR733420

DOI:

10.6100/IR733420

Document status and date: Published: 01/01/2012

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

(2)

Simulating Unsteady Conduit Flows

with

(3)

All rights are reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without prior permission of the author.

A catalogue record is available from the Eindhoven University of Technology Library Proefschrift. - ISBN: 978-90-386-3167-7

NUR 919

Subject headings: initial boundary value problems; moving boundaries; channel flows; two-phase flows; waterhammer; pipe filling and emptying; isolated slug; smoothed particle hydrodynamics

2010 Mathematics Subject Classification: 65D07, 65D10, 65N15, 65N35, 76B10, 76D05, 76D50, 76M28, 76T10

The work described in this thesis has been financially supported by the China Scholarship Council (CSC).

(4)

Simulating Unsteady Conduit Flows

with

Smoothed Particle Hydrodynamics

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Eindhoven, op gezag van de rector magnificus, prof.dr.ir. C.J. van Duijn, voor een

commissie aangewezen door het College voor Promoties in het openbaar te verdedigen

op maandag 25 juni 2012 om 16.00 uur

door

Qingzhi Hou

(5)

prof.dr. R.M.M. Mattheij

Copromotor: dr.ir. A.S. Tijsseling

(6)

Summary

Pipelines are widely used for transport and cooling in industries such as oil and gas, chemical, water supply and sewerage, and hydro, fossil-fuel and nuclear power plants. Unsteady pipe flows with large pressure variations may cause a range of problems such as pipe rupture, support failure, pipe movement, vibra-tion and noise. The unsteady flow is generally caused by flow velocity changes due to valve or pump operation. Water hammer is the best known and exten-sively studied phenomenon in this respect. Fast transients may also occur in rapid pipe filling and emptying processes. Due to high driving heads, the ad-vancing liquid column may achieve a high velocity. When this high-velocity col-umn is blocked or restricted in its flow, high water-hammer pressures may result. Another scenario is that of slug flow, which arguably is the most dangerous type of two-phase pipe flow. Heavy isolated liquid slugs travelling at high speed be-have like cannonballs. Damage is likely to happen when these slugs impact on barriers such as pumps, bends and partially closed valves. Advancing liquid columns occurring in rapid pipe filling and emptying can be seen as a special case of isolated slugs.

In this thesis, we present a Lagrangian particle method for solving the Euler equa-tions with application to water hammer, rapid pipe filling and emptying, and iso-lated slugs travelling in an empty pipeline. As a meshfree method, the smoothed particle hydrodynamics (SPH) used herein is suitable for problems encompass-ing movencompass-ing boundaries and impact events, which are the common features of the concerned topics.

We first present the kernel and particle approximation concepts, which are two essential steps in SPH. Based on numerical approximation rules, the SPH dis-crete form of the Euler and Navier-Stokes equations are derived. To treat var-ious boundary conditions, we apply several types of image particles that are particularly designed to complete the kernels truncated by system boundaries. The global conservation of mass and linear momentum is then demonstrated. The SPH errors in the integral approximation and summation approximation are analysed based on given particle distribution patterns. Other problems such as particle clustering, tensile instability, particle boundary layer and lacking of poly-nomial reproducing abilities (incompleteness) are also discussed together with possible remedies.

(7)

Before applying the implemented particle solver to the thesis topics, we first thor-oughly test it against a selection of two-dimensionale benchmarks, which have a close relationship with the concerned problems. They include dam-break, jet im-pinging onto an inclined plane, emerging jet under gravity, free overfall and flow separation at bends. Good agreements with analytical and numerical solutions in literature are found. The convergence rate of SPH is shown to be of first order, which is consistent with the theoretical analysis.

For the rapid pipe filling problem, we apply the 1D SPH solver to the laboratory experiment of Liou & Hunt [116]. The velocity head at the inlet has to be taken into account to obtain a good agreement with the experiment. Water elasticity does not play a role and the friction formulation for steady state flows can be used. Head transition analysis provides deeper insight into the hydrodynamic behaviour of the filling process. As a special case of pipe filling, water hammer due to liquid impact at partially and fully closed valves is studied. The results agree well with standard MOC solutions. Similar observations are made for the rapid emptying process.

For the isolated slug travelling in a voided pipeline and impacting on a bend, we apply the 1D and 2D SPH solvers to the laboratory experiments of Bozkus [24]. To obtain the arrival velocity of the slug at the elbow, a 1D model including mass loss at the slug tail is used. In the slug impact, flow separation at the bend plays a vital role, which is typical 2D flow behaviour at a geometrical discontinuity. With a flow contraction coefficient obtained from 2D SPH solutions, the improved 1D model gives good results for the reaction force, not only in magnitude but also its duration and shape.

Finally, to study the evolution of air/water interfaces and its possible effect on fill-ing and emptyfill-ing processes, a new experimental study is performed in a large-scale pipeline. It is found that in filling the water front tends to split into two fronts propagating with different velocities. This results in air intrusion on top of a water platform. In emptying, flow stratification occurs at the water tail. Consequently, the validated 1D assumption of vertical air/water interfaces for small-scale systems with relatively high driving head may not be applicable to large-scale systems. The interface evolution does not play an important role in filling, the overall behaviour of which can be well predicted with 1D SPH solu-tions. However, flow stratification largely prolongs the overall draining process.

(8)

Samenvatting

Pijpleidingen worden gebruikt voor transport en koeling in de olie- en gasin-dustrie, in de chemische ingasin-dustrie, voor watervoorziening en riolering, en in en-ergiecentrales werkende op waterkracht, fossiele brandstof of kernreacties. In-stationaire buisstromingen met grote drukvariaties kunnen een scala aan prob-lemen veroorzaken, zoals lei-dingbreuk, schade aan verankering en ophanging, verplaatsing van buizen, trillingen en lawaai. De instationaire stroming wordt meestal veroorzaakt door veranderingen van de vloeistofsnelheid ten gevolge van het manipuleren van kleppen of pompen. Waterslag is het bekendste en meest bestudeerde verschijnsel op dit gebied. Waterslag kan ook optreden bij het snel vullen en legen van leidingen. Bij een hoge aandrij-vende druk kan de bewegende vloeistofmassa een grote snelheid bereiken. Wanneer zo’n hoge-snelheidsmassa botst op een blokkade of lokale weerstand, kan dit resulteren in waterslag en de bijbehorende hoge drukvariaties. Een ander scenario is dat van propstroming, mogelijk de meest gevaarlijke vorm van twee-fasenstroming. De zware vloeistofproppen die zich verplaatsen met hoge snelheid zijn te vergeli-jken met kanons-kogels. Schade is bijna onvermijdelijk wanneer een dergelijke prop botst op een pomp, een bochtstuk of een gedeeltelijk gesloten kleplichaam. De bewegende vloeistofmassa’s bij het vullen en ledigen van leidingen mogen beschouwd worden als een speciaal geval van de ge¨ısoleerde prop.

Dit proefschrift beschrijft een Lagrangiaanse deeltjesmethode voor het oplossen van de Euler vergelijkingen met toepassing op het gebied van waterslag, het snel vullen en legen van leidingen en de beweging van ge¨ısoleerde vloeistofproppen in een verder lege buis. De gebruikte roostervrije methode ”smoothed particle hydrodynamics” (SPH) is geschikt voor problemen met bewegende randen en botsingen met randen, zoals van belang voor de beschouwde onderwerpen. We introduceren eerst de begrippen kernen deeltjesbenadering, twee essenti¨ele stappen in de SPH methode. Gebaseerd op numerieke benaderingsregels wordt de discrete SPH vorm van de Euler en Navier-Stokes vergelijkingen afgeleid. Bij de behandeling van de gestelde randvoorwaarden passen we verschillende types gespiegelde deeltjes toe die speciaal zijn geconstrueerd om de door randen afge-broken kernen weer volledig te maken. Het globale behoud van massa en impuls wordt dan aangetoond. De fouten gemaakt in de SPH integraal- en sommatiebe-naderingen worden geanalyseerd voor gegeven verdelingen van de deeltjes. An-dere problemen, zoals het samenklonteren van deeltjes, treksterkte instabiliteit,

(9)

randlagen van deeltjes en de onmogelijkheid om polynomen te reproduceren (on-volledigheid) worden besproken samen met mogelijke oplossingen.

Voordat de ge¨ımplementeerde deeltjes-oplosser wordt toegepast op de onderwer-pen van het proefschrift, testen we deze eerst en grondig voor een aantal tweed-imensionale standaardproblemen die sterk gerelateerd zijn aan de beschouwde onderwerpen. Dit zijn: dambreuk, waterstraal botsend op een scheve plaat, fontein, waterkering en loslating van de stroming bij een bocht. De overeenkom-sten met analytische en numerieke resultaten uit de literatuur zijn goed. De eerste-orde convergentiegraad van SPH wordt aangetoond, wat in overeenstem-ming is met de theoretische analyse.

Voor het probleem van het snel vullen van een leiding, passen we de 1D SPH methode toe op het laboratoriumexperiment van Liou en Hunt [116]. Het snel-heidsafhankelijke drukverlies bij de buisinlaat wordt in rekening gebracht om goede overeenkomst met het experiment te krijgen. De elasticiteit van het water speelt zoals verwacht geen rol en een stationaire wrijvingswet kan worden ge-bruikt. Een analyse van energieomzettingen leidt tot dieper inzicht in het hydro-dynamisch gedrag van het vulproces. Als een speciaal aspect van vullen, is wa-terslag ten gevolge van de botsende vloeistofkolom op een gedeeltelijk of volledig gesloten klep bestudeerd. De resultaten berekend met SPH komen goed overeen met de conventionele oplossing verkregen met de karakteristiekenmethode. Voor het snelle ledigen van leidingen zijn vergelijkbare constateringen gedaan. Voor de ge¨ısoleerde prop zich voortbewegend in een lege leiding en botsend op een haaks bochtstuk, passen we de 1D en 2D SPH methodes toe op de laborato-riumexperimenten van Bozkus [24]. Om de snelheid waarmee de prop de bocht raakt te bepalen is een 1D model met massaverlies aan de achterkant van de prop gebruikt. Bij de botsing van de prop met het bochtstuk speelt loslating van de stroming een cruciale rol, wat kenmerkend is voor een 2D stroming rond een ge-ometrische discontinu¨ıteit. Met de invoering van een contractieco¨effici¨ent waar-van de waarde bepaald is met behulp waar-van 2D SPH simulaties, geeft het verbeterde 1D model goede resultaten voor de kracht op de bocht, niet alleen wat betreft de grootte, maar ook wat betreft duur en vorm.

Uiteindelijk, om de evolutie van een lucht-water front en de invloed daarvan op het vullen en legen te bestuderen, is er een nieuwe experimentele studie op grote schaal uitgevoerd in een lange pijpleiding. Bij het vullen is waargenomen dat het water-front zich splitst in twee fronten die zich met verschillende snelheden verplaatsen. Dit resulteerde in het indringen van lucht boven een min of meer vlakke laag water. Bij het legen treedt deze stratificatie op aan de achterkant van de waterkolom. Het gevolg is dat de 1D aanname van verticale lucht-water fron-ten, zoals gevalideerd in experimenten op kleine schaal met relatief hoge aan-drijvende druk, niet toepasbaar is voor grotere leidingen. De vervorming van het lucht-water front speelt geen belangrijke rol bij het vullen van leidingen: het globale gedrag kon goed worden voorspeld met de 1D SPH oplossing. Daarente-gen verlengt het ontstaan van een horizontale water-lucht laag wel de tijdsduur van het legen van leidingen.

(10)

Introduction

1.1 Motivation for the study

Piping systems are widely used for the transportation of fluids in different ap-plications such as power stations, water supply and sewerage works, oil and gas industry, and chemical plants. In these systems, many accidents and incidents occur [108]. For example, in 2009 a terrible accident happened in the Sayano-Shushenskaya hydropower plant, the largest power plant in Russia, and led to grave consequences. The turbine hall and engine room were flooded due to pipe rupture, the ceiling of the turbine hall collapsed, 9 of 10 turbines were damaged or destroyed, and 75 people were killed [189]. Frequent incidents include displaced and broken pipes, failure of supports and hydraulic machinery, noise, vibration and fatigue.

What are the reasons for these accidents and incidents happening around us? Except for external excitations, such as earthquakes, wind loads and vehicle im-pacts, the most common reason is unsteady fluid dynamics within the pipeline. There are three special cases that deserve attention due to their violent appear-ance:

Water hammer. When fluid flowing in a pipe is suddenly stopped due to valve operation or pump failure, a large pressure rise results and excites the system. The pressure waves in the system impact the valve or pump like a hammer [230]. For illustration, we refer to a personal experience. During the pipe filling and draining experiments, described in Chapter 7, the simultaneous closure of control valves caused severe water-hammer in the experimental test rig [22]. Several joints close to an elbow and to a U-bend were disconnected (see Figs. 1.1 and 1.2b) and pipe anchors screwed into the floor were lifted up (see Fig. 1.2a). Five steel bars supporting the U-bend were displaced and one of them broke off and

(15)

was knocked down to the ground floor 5 metres below (see Fig. 1.2b).

Figure 1.1: Joint failure in the PVC pipeline apparatus [22].

Figure 1.2: Support failures in the PVC pipeline apparatus [22].

Rapid filling and emptying of a pipeline. Take a sewage system as an exam-ple. During rain periods, water from roofs and roads flows through manholes to the underground sewer system. When the advancing water flow is blocked somewhere, a regime transition from open channel flow to pressurized flow takes place. Damage is easily induced [236, 237], because sewage systems are not de-signed to take any overpressure. Surface and basement flooding due to function loss may result. In addition, when a new pipe system is taken into operation, the

(16)

filling (system start-up) has to be done very carefully. The same holds for system cleaning. In rapid pipe filling and emptying, violent unsteady flows may occur. Isolated slugs in an empty line. The scenario of this problem is for example as fol-lows. When steam lines in power plants are stopped for maintenance, condensed steam accumulates in lower locations and forms water slugs. When the system is reactivated, the water slugs propelled by highly compressed air or steam may achieve very high velocities and propagate like ”bullets” in the pipeline. When such ”bullets” impact on an elbow, pump or partially closed valve, damage may easily be the result [170].

These three violent cases make up the topics in this thesis.

1.2 Problem statement and method justification

Water hammer, rapid filling and emptying, and travelling slugs in pipelines are sketched in Fig. 1.3. These are simplified physical models or laboratory test rigs. The three cases have several physical features in common. First, since the fluid velocity is generally high, the flows are turbulent with large Reynolds number. Second, as two components – air and water – are involved in pipe filling and emp-tying, slug flow and sometimes water hammer, they are two-phase flows with pos-sible flow-regime transitions. Third, they are unsteady flows with moving bound-aries confined within slender structures. Mathematically they can be modelled as initial boundary value problems governed by the Navier-Stokes (N-S) equations. What makes them special are the moving boundaries and contact discontinu-ities due to fluid-solid impact. Special attention is needed to understand local flow behaviour, such as moving and deforming water fronts and flow separation at bends.

For flows in slender pipes and channels, one-dimensional models, such as rigid-column models and elastic models, are preferred because of low computational cost and acceptable accuracy. Rigid-column models have been applied to rapid pipe filling [27, 116, 183] and emptying [105], and travelling slugs [24–26, 98, 170, 231]. When the flow impacts a solid or is suddenly accelerated, an elastic model has to be used [136, 137, 219].

One-dimensional models cannot adequately describe flow stratification, flow sep-aration, impact on a bend and water front evolution as encountered in the consid-ered problems. Two- and three-dimensional models are highly needed. However, flows in two and three dimensions with moving and deforming free boundaries are difficult to simulate using traditional mesh-based methods with either Eule-rian or Lagrangian grid. The fluid-solid impacts are hard to tackle. Therefore, our approach is based on moving ”particles”. In particular, we use the smoothed particle hydrodynamics (SPH) numerical method. SPH is a proper choice here

(17)

Figure 1.3: Sketches for the problems considered in this thesis: (a) water hammer [211], (b) rapid pipe filling, (c) rapid pipe draining and (d) isolated slug in a void pipe [24].

because of its Lagrangian, mesh-less and inherently elastic features. The exten-sion from one-dimenexten-sional SPH to two and three dimenexten-sions is straightforward, but expensive. SPH is a good method for impact simulations too [93, 117].

(18)

Using the SPH method, we model the water-hammer problem, rapid pipe filling and emptying, and isolated slug motion in a pipe and its impact on a bend. Even in 1D setting, the moving particle model is already remarkably versatile and gives satisfactory results when compared to physical experiments. To improve existing 1D models, 2D SPH simulations give crucial insight into local flow behaviour, such as free-surface evolution, flow separation, contraction and stratification, and velocity and pressure distribution. Flow stratification and deformation of water front and tail occur in rapid pipe filling and draining. This is confirmed by the novel large-scale laboratory measurements reported herein. Flow separation at a sharp bend turns out to be vital to accurately predict slug impact (magnitude and duration) in both the 2D and improved 1D models.

Important basic topics include mathematically deriving and evaluating the SPH method, applying it to the N-S equations to establish the SPH fluid dynamics equations, investigating and extending the SPH treatment of various boundary conditions, estimating the accuracy of the SPH approximations (both the function and its derivatives), analyzing the incompleteness of SPH and finding possible improvements.

1.3 Outline of the thesis

The outline of this thesis is as follows. In Chapter 2, we present the mathemat-ical modelling of the concerned problems, concentrating on the Navier-Stokes (N-S) equations in Lagrangian form and the formulations of various boundary conditions. The smoothed particle hydrodynamic (SPH) method is introduced in Chapter 3 for solving the N-S equations. Different approximation rules are evaluated and a formulation fulfilling the conservation properties is derived. To properly complete the kernels associated with particles close to the boundary, the implementation of the image particle approach is detailed. Based on given par-ticle distributions, the numerical error in SPH is estimated. New algorithms to improve the completeness condition are proposed. Also, in Chapter 3, essential features of the numerical implementation for the particle solver are presented. In Chapter 4, the coded SPH solver is verified against several selected 2D bench-marks that are closely related to the defined problems. The procedure to set up a numerical experiment is introduced, with focus on how to choose the key pa-rameters. The SPH numerical convergence rate is checked and compared with the theory. The application of the 1D SPH solver to rapid pipe filling and emp-tying is worked out in Chapter 5. The SPH results are validated against exper-iments of Liou & Hunt [116] and numerical solutions found in literature. The rigid-column model is established and validated. Fluid dynamics in the filling and emptying of conduits is examined through numerical experiments. The SPH results for the water-hammer problem are validated against the exact solution by the method of characteristics. Chapter 6 presents the application of the 1D and 2D SPH solver to the isolated slug experiment of Bozkus [24]. Pressure loss and

(19)

reaction force at the bend are estimated from the flow separation results obtained in Chapter 4. Chapter 7 deals with the large-scale laboratory experiments per-formed at Deltares, Delft. The water-hammer experiments have been published in [23, 99] and a general description of the filling and emptying experiments has been reported in [105]. The focus of Chapter 7 is on the water-air and air-water in-terface evolutions in the filling and emptying processes. The 1D SPH model and the rigid-column model developed in Chapter 5 are used to interpret the new experiments. The considered problems, models and methods, and the related topics are so many and diverse that an overall literature review for all of them is not possible. The review of a specific problem is given where it is studied.

(20)

Mathematical Modelling

2.1 Introduction

For the flow phenomena described in Chapter 1, the continuum mechanics prin-ciple can still be applied and the Navier-Stokes equations are valid. Before the derivation of the equations is carried out, several physical and mathematical con-cepts are introduced in Section 2.2. With the model of an infinitesimally small fluid element moving with the flow, the continuity equation and momentum equation are derived in Sections 2.3 and 2.4, respectively. The equations for vis-cous and inviscid flow are summarized in Section 2.5 in the notation used in SPH literature. The physical boundary conditions and the mathematical representa-tions are described in Section 2.6. Section 2.7 describes the treatment of a moving domain.

2.2 Basic concepts and assumptions

In obtaining the basic equations governing the fluid motion, the following phi-losophy is adopted:

• Choose the appropriate fundamental principles from the laws of physics. • Apply these physical principles to a suitable model of the fluid flow. • Extract the mathematical equations which state the physical principles. From the physical principles, the mass and momentum conservation laws will be used to derive the governing equations. Four alternative models of the flow

(21)

are constructed in [3]. Here the model of an infinitesimally small element moving with the fluid will be used, because this is in close relation with the basics of SPH. This will be further explored in Chapter 3.

2.2.1 Infinitesimal fluid element moving with the flow

In mathematical physics, a physical model is needed first, from which a corre-sponding mathematical model can be derived. In this section, we discuss a model of the flow which is in particular useful to understand the ”particle” concept that will be elaborated in Chapter 3.

Figure 2.1: A flow field represented by streamlines.

Consider a flow field represented by streamlines as shown in Fig. 2.1. An infinites-imal fluid element moving with the flow will be extracted and examined. The fluid element is very small but large enough to contain a huge number of molecules so that the physical principles for continuum mechanics are applicable. The fluid element has a fixed mass, but its volume will be time-dependent. Since the flow field is described by fluid elements, the physical principles will be applied to one infinitesimal fluid element instead of the whole flow field. The fundamental par-tial differenpar-tial equations (PDEs) will then be derived in what is generally called the non-conservation form.

2.2.2 The material derivative

The material derivative is an important notation and concept used in SPH. It has an important physical meaning that will be helpful to understand the basic

(22)

ideas of SPH. To understand its physical meaning, an infinitesimally small fluid element moving on a pathline in a three-dimensional domain is shown in Fig. 2.2. Then the material derivative is expressed in vector notation as:

D

Dt:=

∂

∂t+v· ∇, (2.1)

where t is time and∇ := i ∂

∂x+j

∂

∂y+k

∂

∂z in a Cartesian coordinate system.

Figure 2.2: Fluid element moving in a flowing fluid – sketch for the material derivative.

From the definition (2.1), the material derivative D/Dt is the time rate of change in a moving fluid element. The first term on the right hand side (RHS) of (2.1) is the local derivative, which is the time rate of change at a fixed point. The second term is the convective derivative, which is the spatial rate of change due to the movement of the fluid element in the flow field. The material derivative can be applied to any scalar field variable in the flow, such as density, velocity components, pressure, etc. The material derivative is essentially the same as the total derivative with respect to time [3] in calculus.

2.2.3 The divergence of the velocity

The physical meaning of the divergence of the velocity given in this section is helpful to understand the mass conservation property of the SPH equations. Assume that we are dealing with a velocity field v in a three-dimensional space

in a Cartesian coordinate system. The symbol∇ · v is defined by

(23)

where v :={u, v, w}T. Equivalently, it can be expressed as [3]:

∇ · v = 1

δV

D(δV)

Dt (2.2)

where δV is the volume of the infinitesimal fluid element as sketched in Fig. 2.2.

The relation (2.2) clearly demonstrates the physical meaning of∇ · v. The

diver-gence of the velocity is the time rate of change of the volume of a moving fluid element, per unit volume.

2.3 The continuity equation

The mass conservation law will be applied to the infinitesimally small fluid ele-ment moving with the flow, resulting in the continuity equation.

Consider the flow model sketched in Fig. 2.2. Assume that the fluid element has

a fixed mass δm and a time-dependent volume δV. Then

δm = ρδV, (2.3)

in which ρ is the mass density of the fluid. The mass conservation law is stated as the rate of change of δm being zero. According to the physical meaning of the material derivative introduced in Section 2.2.2, we have

D(δm)

Dt = 0. (2.4)

Replacing mass by density and volume gives Dρ Dt + ρ [ 1 δV D(δV) Dt ] = 0. (2.5)

The term in brackets in Eq. (2.5) is equal to ∇ · v given by Eq. (2.2). Hence,

combining Eqs. (2.2) and (2.5), we obtain the continuity equation Dρ

Dt + ρ∇ · v = 0. (2.6)

2.4 The momentum equation

Newton’s second law can be applied to the infinitesimal fluid element moving with the flow. The momentum equation is then derived. The moving fluid ele-ment with forces exerted on it is sketched in Fig. 2.3. Under action of these forces the fluid element will change its position, shape and volume. When applied to the moving fluid element in Fig. 2.3, Newton’s second law states that the net force

(24)

on the fluid element equals its mass times its acceleration. Although it is a vector relation, only the x component of Newton’s second law is considered, i.e.

Fx= δmax, (2.7)

where Fxis the x component of the net force, δm is the mass of the fluid element,

and axis the x component of the acceleration.

Figure 2.3: Forces on the moving fluid element used in the derivation of the x component of the momentum equation.

The forces acting on the fluid element include body forces and surface forces.

Assume that the body force per unit mass in x-direction is fx. The surface forces in

the x-direction exerted on the fluid element are depicted in Fig. 2.3. The symbol

τij denotes a stress in the j direction exerted on a plane perpendicular to the i

axis. The total net force in the x-direction is written as

Fx= [ −∂p ∂x + ∂τxx ∂x + ∂τyx ∂y + ∂τzx ∂z ] dx dy dz + fxρ dx dy dz. (2.8)

The fixed mass of the fluid element is δm = ρdxdydz. Since axis the material

rate of change of the velocity u, we have ax := Du/Dt due to the Lagrangian

approach. By substituting Fx, δmand axinto Eq. (2.7), we obtain

ρDu Dt = − ∂p ∂x + ∂τxx ∂x + ∂τyx ∂y + ∂τzx ∂z + ρfx. (2.9)

(25)

and z components of the momentum equation read ρDv Dt = − ∂p ∂y+ ∂τxy ∂x + ∂τyy ∂y + ∂τzy ∂z + ρfy (2.10) and ρDw Dt = − ∂p ∂z + ∂τxz ∂x + ∂τyz ∂y + ∂τzz ∂z + ρfz. (2.11)

Newtonian fluids are defined as fluids where the shear stress is proportional to the time rate of strain, i.e., the spatial derivative of velocity. The fluids considered in this thesis are assumed to be Newtonian fluids. For such fluids, we have

τxx= λ∇ · v + 2µ ∂u ∂x, (2.12) τyy= λ∇ · v + 2µ ∂v ∂y, (2.13) τzz= λ∇ · v + 2µ ∂w ∂z, (2.14) τxy= τyx= µ ( ∂v ∂x+ ∂u ∂y ) , (2.15) τxz= τzx= µ ( ∂w ∂x + ∂u ∂z ) , (2.16) τyz= τzy= µ ( ∂w ∂y + ∂v ∂z ) , (2.17)

where µ is the dynamic viscosity coefficient of the fluid and λ is the second vis-cosity coefficient. The dynamic visvis-cosity coefficient is a measure of the internal molecular friction of the fluid (momentum diffusion). The following hypothesis (Stokes) is generally used for λ:

λ = −2

3µ.

Substituting all the shear forces into (2.9) – (2.11), we obtain the three momentum equations ρDu Dt = − ∂p ∂x + ∂ ∂x ( λ∇ · v + 2µ∂u ∂x ) + ∂ ∂y [ µ ( ∂v ∂x+ ∂u ∂y )] + ∂ ∂z [ µ ( ∂w ∂x + ∂u ∂z )] + ρfx, (2.18) ρDv Dt = − ∂p ∂y + ∂ ∂x [ µ ( ∂v ∂x+ ∂u ∂y )] + ∂ ∂y [ λ∇ · v + 2µ∂v ∂y ] + ∂ ∂z [ µ ( ∂w ∂y + ∂v ∂z )] + ρfy (2.19)

(26)

and ρDw Dt = − ∂p ∂z + ∂ ∂x [ µ ( ∂w ∂x + ∂u ∂z )] + ∂ ∂y [ µ ( ∂v ∂z + ∂w ∂y )] + ∂ ∂z ( λ∇ · v + 2µ∂w ∂z ) + ρfz. (2.20)

2.5 Summary of the governing equations

By applying two physical conservation laws to a moving fluid element, the con-tinuity and momentum equations have been derived in Sections 2.3 and 2.4, re-spectively. The equations for viscous and inviscid flow are summarized in this section in the notation used in SPH literature.

2.5.1 The Navier-Stokes equations for viscous flow

For viscous flow the governing equations are the Navier-Stokes equations. The Greek superscripts α, β and ϑ are introduced to denote the coordinate directions. The summation in the equations is indicated by repeated indices. Then the governing equations as used for SPH in Chapter 3 are:

Continuity equation Dρ Dt = −ρ ∂Vα ∂xα, (2.21) Momentum equation DVα Dt = 1 ρ ∂σαβ ∂xβ , (2.22)

where Vα _{is the α-th component of the fluid velocity v, and σ is the total stress}

tensor consisting of the isotropic pressure p and the viscous stress τ,

σαβ= −pδαβ+ ταβ, (2.23)

where δαβ_{is the Kronecker delta and τ}αβ_{= µε}αβ_with

εαβ= ∂V α ∂xβ + ∂Vβ ∂xα − 2 3 ∂Vϑ ∂xϑδ αβ . (2.24)

2.5.2 The Euler equations for inviscid flow

The governing equations for inviscid flow are known as the Euler equations. By dropping all terms involving viscosity from the Navier-Stokes equations, the Eu-ler equations are obtained as

(27)

Continuity equation Dρ Dt = −ρ ∂Vα ∂xα, (2.25) Momentum equation DVα Dt = − 1 ρ ∂p ∂xα. (2.26)

2.5.3 Equation of state

To close the system of equations for either viscous or inviscid flow, additional information is provided by the equation of state. The Tait-Murnaghan equation of state [162] relating pressure p and density ρ is given by

p − p0= ρ0c20 γ (( ρ ρ0 )γ − 1 ) , (2.27)

where p0is the reference pressure, ρ0is the reference density at reference pressure,

c0 is the speed of sound, and γ is a positive constant. In this thesis atmospheric

pressure is taken as the reference pressure p0.

2.6 Initial and boundary conditions

The equations governing the flow have been presented and explained in the pre-vious sections. To obtain a particular solution to these equations, initial and boundary conditions are needed. In this section, the physical initial and bound-ary conditions used in later chapters are described.

Initial conditions

Initially (at t = 0), we impose a given velocity v0and pressure to the flow, i.e.

v = v0, p = p(t = 0). (2.28)

The initial density ρ(t = 0) is calculated from Eq. (2.27). Boundary conditions

Several types of boundaries are considered and four of them are described here for free-surface flows. They are the wall boundary, free surface boundary, in-flow boundary and outin-flow boundary. Other boundaries such as supply reser-voir, conduit bend, partially closed valve, etc will be described when they are encountered in specific problems. In order to efficiently formulate the boundary

(28)

is the exterior normal direction),Vtthe velocity component parallel to the

bound-ary (t is the tangential direction), and ∂Vn/∂nand ∂Vt/∂ntheir derivatives in

the normal direction. (i) Wall boundary

With respect to the fixed wall boundary, the condition is different for viscous and inviscid flows. For a viscous flow, the boundary condition is the no-slip condition, i.e. fluid does not penetrate through the boundary and the fluid at the wall is at rest, i.e.

Vn = 0 and Vt= 0. (2.29)

For an inviscid flow, there is no skin friction to induce its ’sticking’ to the surface. Hence the flow velocity at the wall is a finite, non-zero value. Moreover, there cannot be mass flow into or out of a non-permeable wall. The flow velocity vector immediately adjacent to the wall must therefore be tangent to the wall. Conse-quently, the wall boundary condition for an inviscid flow is given by

Vn = 0 and ∂Vt/∂n = 0. (2.30)

The first condition of (2.30) implies that there is no flow penetrating through the boundaries, whilst the second implies that there are no frictional losses at the boundary. Apart from being enforced for the frictionless wall boundary, the free-slip condition is often imposed along a line or plane of symmetry, thereby reducing the size of the domain where the flow needs to be calculated by a half.

(ii) Free surface boundary

On a free-surface a kinematic and a dynamic condition are required. The kinematic boundary condition is written as

∂Vt/∂n = 0, (2.31)

which is the same as the second condition of (2.30). The difference is that n is the exterior normal direction to the free-surface which is generally changes with time before a steady flow is reached. This condition mathematically states that in the

exterior normal direction the velocity component Vthas a maximum value at the

free surface. Physically, it means that the particles of air at the surface basically move at the same speed as the particles of water at the surface, because there is no shear on the surface. It also implies that a fluid particle initially on the surface will remains on it during the evolution of the surface. Since surface tension is not taken into account and the atmospheric pressure is taken as the reference pressure

p0, the dynamic boundary condition is

p = 0. (2.32)

(iii) Inflow boundary

For the inflow boundaries considered herein, the velocity is prescribed at the in-let. For an inflow boundary with a given velocity, the condition is

(29)

A constant and uniform pressure is also described at the inflow section. (iv) Outflow boundary

At the outlet, the velocity is assumed not to change in the direction normal to the boundary. Accordingly, the outflow conditions are given by

∂Vn/∂n = 0, and ∂Vt/∂n = 0. (2.34)

2.7 Moving fluid domain

For free-surface flows the domain occupied by the fluid may change with time. In addition to the pressure and the velocity, the domain Ω(t), which the fluid occupies at time t, is also to be determined. Besides the governing equations and the initial and boundary conditions, the problem description also requires

the initial configuration Ω0. Starting from it, the fluid domain evolves with time

(30)

Smoothed Particle

Hydrodynamics (SPH)

In this chapter, the SPH method and its variations are described. The errors in the SPH approximations of a function and its derivatives are estimated. A de-tailed outline is as follows: In Section 3.1, a brief introduction to meshfree meth-ods is given. As one of the earliest meshfree methmeth-ods, SPH and its variations are reviewed in Section 3.2. The basic concepts and formulations of SPH are scribed in Section 3.3. In Section 3.4 the SPH equations for fluid flows are de-rived and their conservation properties are shown. A numerical implementation is also derived with much attention to the treatment of boundary conditions. The numerical errors in the SPH integral approximation and in the summation ap-proximation are estimated in Section 3.5, where corrective methods aiming at the improvement of the consistency in SPH approximations are presented.

3.1 Introduction

Fluid flow can be descried by both Eulerian and Lagrangian approaches [117]. Conventional numerical methods such as the finite difference method (FDM), fi-nite volume method (FVM) and fifi-nite element method (FEM) partition the region

of interest into a finite number of (small) cells, associate degrees of freedom (DOFs)

with each cell and lead to a (non)linear system in theseDOFs. In an Eulerian

ap-proach the collection of cells and associatedDOFs (called mesh) is static in space.

In a Lagrangian approach the numerical mesh follows the material motion. For Eulerian grids, various difficulties appear in modelling complex geometries, free surface flows, moving interfaces, deformable boundaries, etc. For Lagrangian grids, difficulties occur in dealing with problems with large deformation,

(31)

pene-tration and crack propagation, as mesh connectivity and topology changes. For problems with complex geometry, the generation of a quality mesh and proper mesh refinements are computationally expensive. Moreover, mesh-based methods are not suitable for simulations of explosions, high velocity impacts and discrete particle systems.

Over the past decades, a new class of computational methods has been devel-oped, which approximate the governing equations only based on a set of uncon-nected nodes without the need for an additional mesh. They are usually referred to as meshfree or meshless methods and have been a major research focus starting from the middle of the 1990s. For the basics of the mathematical theory and nu-merical implementation of different meshfree methods, see the review papers by, among others, Monaghan (1992) [151], (2005) [156], Belytschko et al. (1996) [16], (1998) [17], Randles & Libersky (1996) [182], Li & Liu (2002) [111], Babuˇska et al. (2003) [7], (2004) [8], Fries & Matthies (2004) [63], Nguyen et al. (2008) [164], and Liu & Liu (2010) [121]. The work of Monaghan [151, 156] and Liu & Liu [121] is about smoothed particle hydrodynamics (SPH) and its applications, something that we focus on herein. Randles & Libersky [182] summarized the use of SPH in solid mechanics and recent developments in this field are reviewed by Liu & Liu [121]. Fries & Matthies [63] classified the meshfree methods and described them from their different origins and viewpoints. For meshfree methods based on weak forms, a numerical implementation similar to FEM is used. Babuˇska and his co-workers [7] discussed meshfree methods (weak form) from the mathemat-ical point of view. The core technique of meshfree methods (weak form) is the partition of unity of the shape functions.

Two regular international workshops on meshfree methods have been established. One is the annual workshop SPHERIC [197] which started in 2006. It is mainly on SPH and its variations (strong form). The other is the biennial workshop Meshfree Methods for PDEs [144] which started in 2003 and focuses on the rig-orous mathematical background of meshfree methods typically based on weak forms. The SPH method has been implemented in LS-DYNA [131] as a commer-cial toolbox for high-velocity impact (HVI) problems. There are at least three open source codes available around the world. The first one can be found via the link ’http://www.liugr.org’. This code is in Fortran 77 and developed by Liu & Liu accompanying their textbook [117]. The code is clear and useful for beginners, but less efficient because of an ineffective searching algorithm. The second code named HYDRA can be found from ’http://www.pearcef.org’. The code has been developed by Couchman et al. accompanying the paper [45] for applications in astrophysics. It is written in Fortran 95. For the application of SPH in ocean engi-neering with large-scale free surfaces, the code ”SPHERIC” has been developed by the SPHERIC committee and this can be found via ’http://www.spheric.org’. All three codes have their own features and specific application targets. Since none of them is proper for the problems considered in this thesis, we developed our own 1D codes and worked on the 2D codes developed by Dr. Kruisbrink and Dr. Pearce (one of the developers of HYDRA [45, 135, 172, 209]) of the University of Nottingham. It is an interaction between Matlab 2010b and Fortran 95.

(32)

3.2 Literature review on SPH

SPH methodology. SPH is a meshless, Lagrangian, particle method that uses

an approximation technique to calculate field variables like velocity, pressure, po-sition, etc. In SPH the governing PDEs for fluid dynamics are directly trans-formed into ordinary differential equations (ODEs) by constructing their integral forms with a kernel (or smoothing) function and its gradient. Unlike FDM, FVM and FEM, the SPH method uses a set of particles without predefined connec-tivity to represent a continuum system and thus it is easy to handle problems with complex geometries. It does not suffer from mesh distortion and refinement problems that limit the usage of FEM for large deformation problems, HVIs in solid mechanics (where solids are treated fluid-like but with material strength), and hydrodynamic problems with free surfaces and moving boundaries. As a Lagrangian method, SPH naturally tracks material history information due to movement of the particles. It is an ideal alternative for attacking fluid dynam-ics problems. It has the strong ability to incorporate complex physdynam-ics into the SPH formulations. It is easy to work with and acceptable accuracy is usually obtained. The main disadvantage is that it is more time consuming than conven-tional methods and parallel computation techniques are needed if many particles are used. SPH differs from other meshfree methods in two aspects. First, ”par-ticles” (nodes) are moving with the fluid, thus representing a fluid flow. They offer nice visualization (animation) possibilities. The particles satisfy the Navier-Stokes equations, so that SPH is a proper method for fluids. In other meshfree methods [63, 164], the weak form of a governing equation and a Galerkin proce-dure of FEM are generally applied. In these other methods, fixed points with-out priori connectivity are treated as nodes as in FEM. They are mainly used in solid mechanics with small deformations [164]. Second, SPH approximation tech-niques are directly used to discretize the governing equations (Euler or Navier-Stokes equations) and all field information such as position, velocity and density are carried by particles. However, the calculated nodal values in other meshfree methods for solids are not the required displacements and a back substitution is usually needed, which is a typical FEM procedure.

SPH history. SPH was invented independently by Lucy [133] and Gingold & Monaghan [69] for three-dimensional astrophysical problems, such as the forma-tion and evoluforma-tion of proto-stars and galaxies. Since the moforma-tion of astrophysical particles satisfies the same conservation laws as gas flows, it was modelled by gas dynamics. SPH is a truly meshfree Lagrangian method based on strong for-mulations. During its early development (1980s), it was mainly applied to com-pressible astrophysical and cosmological flows. Monaghan [151] presented an el-egant review for applications in this field. Today, SPH is being used to model the collapse and formation of galaxies [20, 45, 173], coalescence of black holes with neutron stars, detonations in white dwarfs and even the evolution of the uni-verse. Recently Rosswog [185] presented a detailed review of the SPH method with particular focus on its astrophysical applications. At the beginning of the 1990s, Libersky & Petschek [114] extended SPH to solid dynamics, especially to

(33)

HVI problems. Monaghan (1994) [152] developed a concept of artificial com-pressibility and used an equation of state for liquids with artificially low speed of sound to model incompressible free-surface flows. After that, SPH was used to model a vast range of fluid dynamics problems. The applications of SPH in the fluid dynamics field include low Reynolds number flows [159], viscous and heat conducting flows [41, 158], multiphase flows [87, 88, 153], coastal hydrody-namics including water wave impact and sloshing [71] and many more [121]. Throughout the development of the SPH method, several benchmarks have been systematically investigated such as Sod’s shock tube, blast waves and dam break. SPH drawbacks. Regardless of a great amount of success achieved with SPH, several technical drawbacks had to be overcome before it was fully developed. Among the shortcomings, four aspects, namely, boundary deficiency, nodal in-completeness, zero–energy mode and tensile instability are well known. These drawbacks often result in a loss of accuracy and stability, and a failure of con-vergence in the SPH calculations. To address these problems, a number of useful techniques have been developed, such as symmetrising [150], normalization [93] and regularization [182] of the first-derivative approximations without changing the kernels, employment of ghost-particles [205], use of extra stress points [56], and introducing artificial pressure terms [72, 155]. Interested readers are referred to Belyschko et al. [17] and Liu et al. [119, 120] for boundary deficiency and nodal completeness, Liu & Liu [117] and Vignjevic [223] for zero–energy mode, and Monaghan [155] and Liu & Liu [121] for tensile instability.

3.3 SPH approximations

3.3.1 Function approximation

As a meshfree method, SPH is built on a set of distributed particles in a domain without a grid or mesh. It consists of two steps of approximations. One is the kernel approximation and the other is the particle approximation.

A continuous scalar field f : Ω⊂ R3_{→ R can be written in integral form as}

f(x) =

∫

Ω

f(x′)δ(x − x′)dx′, ∀ x ∈ Ω, (3.1)

where x → δ(x − x′)is the Dirac delta distribution. The domain Ω has boundary

∂Ω. The Dirac delta distribution is difficult to use as a function in a collocation process [112]. As a main technical ingredient of SPH, a continuous function W(x−

x′, h)known as kernel is chosen to mimic the Dirac delta, such that for h → 0,

W(x − x′, h)→ δ(x − x′). The kernel W is assumed to be radial symmetric in its

support domain Γ , which is totally within Ω (Γ ⊂ Ω) and has boundary ∂Γ. The

(34)

size of the kernel W. The approximated function is then written as ⟨f⟩ (x) :=

∫

Ω

f(x′)W(x − x′, h)dx′, (3.2)

where⟨·⟩ denotes an approximation operator. For a vector function f = [f1, . . . , fn],

the approximation is componentwise and we have⟨f⟩ = [⟨f1⟩ , . . . , ⟨fn⟩]. This step

is known as kernel approximation or integral approximation.

Suppose that space Ω is partitioned into N parts and each part has a volume Ωb.

We have Ω =∪_bΩband Ωa ∩ Ωb=∅, a ̸= b. From (3.2) we get ⟨f⟩ (x)=. ∑ b f(xb)W(x − xb, h)Ωb=: bf(x). _(3.3)

Note that ”volume” Ωbis an area segment in two-dimensions and a line segment

in one-dimension. In much of the SPH literature, the ’hat’ in the approximated

function bf(x)is often deleted without further explanation, but herein we

distin-guish them for the convenience of error estimation in Section 3.5.1.

Assume that the continuum has a mass density distribution ρ(x). A point ’b’ is chosen to represent one part, which is referred to as a particle hereafter. Then the mass of particle b is

mb= ρbΩb ⇔ Ωb= mb/ρb. (3.4)

where ρb= ρ(xb). Using (3.4) and collocation, (3.3) is rewritten as

bfa= ∑ b fbWabΩb= ∑ b mb ρb fbWab, (3.5)

where bfa := bf(xa), fb := f(xb)and Wab := W(xa−xb, h). This step is called

particle approximation or summation approximation.

The idea of the SPH particle approximation for a 2D function is sketched in Fig. 3.1. Particles b within domain Γ are referred to as neighbours of particle a. One particle and its neighbour are interacting particles and form an interaction pair.

Remark 3.1. If we have _∫

Ω

W(x − x′, h)dx′= 1, (3.6)

then (3.2) always holds for constant functions. Condition (3.6) is referred to as normal-ization condition of the kernel.

Remark 3.2. The factor W(x−xb, h)Ωbis a weight function in the approximation (3.3)

of the function values at scattered points, i.e. f(xb).

With respect to the kernel approximation, it is easy to prove that the operator⟨·⟩

is linear. If f1and f2are smooth functions and C is a constant, then

(35)

Figure 3.1: Sketch for SPH particle approximation: (a) side view and (b) top view.

⟨Cf2⟩ = C ⟨f2⟩ . (3.8)

The operator b· is also linear.

Theorem 3.1. The SPH kernel approximation of the product of two smooth functions f1and f2is generally different from the product of the kernel approximation of these two

functions1_{, i.e.}

⟨f1f2⟩ ̸= ⟨f1⟩ ⟨f2⟩ . (3.9)

(36)

Proof. Apparently when h→ 0 or any one of two functions f1and f2is constant,

two sides of (3.9) will be exactly equal. However, when h > 0,⟨f1f2⟩ = ⟨f1⟩ ⟨f2⟩

does not hold even for a linear function. Take scalar functions f1(x) = f2(x) =

x as an example. As W(x − x′, h)is radial symmetric in Ω, together with the

normalization condition (3.6), we have ⟨f1f2⟩ (0) = ∫ Ω (f1f2)(x′)W(0 − x′, h)dx′ = ∫ Ω (x′)2W(0 − x′, h)dx′ = ∫ Ω (0 − x′)2W(0 − x′, h)dx′ > 0, and ⟨f1⟩ (0) = ⟨f2⟩ (0) = ∫ Ω x′W(0 − x′, h)dx′= 0. Consequently, ⟨f1f2⟩ ̸= ⟨f1⟩ ⟨f2⟩ .

3.3.2 Gradient and divergence

Aside from satisfying the conditions such as even function and normalization, it

is assumed that kernel W(x − x′, h)has some other properties:

• It has a compact support Γ such that

W(x − x′, h) = 0 ∀ x /∈ Γ. (3.10) • Its gradient satisfies the antisymmetry condition [112,117]

∫

Ω

∇W(x − x′_{, h)d}_x′₌_0. _(3.11)

For a 3D domain ∂Γ is the surface of a sphere. Define η = x − x′, then in

compo-nentwise, we have ∇xW(η, h) =∇ηW(η, h) dη dx =∇ηW(η, h), and ∇x′W(η, h) =∇ηW(η, h) dη dx′ = −∇ηW(η, h). Thus ∇xW(x − x′, h) = −∇x′W(x − x′, h). (3.12)

(37)

If (∇f)(x) is taken as a function and substituted into (3.2), utilising integration by parts, one can show that

⟨∇f⟩ (x) = ⟨∇xf⟩ (x) = ∫ Ω ( ∇x′f)(x′)W(x − x′, h)dx′ = f(x′)W(x − x′, h)_∂Ω− ∫ Ω f(x′)∇x′W(x − x′, h)dx′ = (3.12) ∫ Ω f(x′)∇xW(x − x′, h)dx′ . = (3.4) ∑ b f(xb)∇W(x − xb, h)Ωb = ∑ b f(xb)∇W(x − xb, h) mb ρb =: c∇f(x), (3.13)

If we differentiate both sides of (3.2), the gradient of the approximated function is obtained as ∇ ⟨f⟩ (x) := ∫ Ω f(x′)∇W(x − x′, h)dx′, (3.14) or in summation form ∇bf(x) :=∑ b mb ρb f(xb)∇W(x − xb, h). _(3.15)

Comparing (3.13) and (3.14), one sees that

⟨∇f⟩ = ∇ ⟨f⟩ ⇒ c∇f = ∇bf, (3.16)

which means that the SPH approximation of the gradient is commutative. This is

an important property and it is an exact relation for Γ ⊂ Ω.

Remark 3.3. The factor∇W(x − xb, h)Ωbis a weight function in the approximation

(3.13). It implies that the derivative of a function can be approximated from neighbouring function values.

SPH Rule I. A concise notation of (3.15) is c ∇fa= ∑ b mb ρb fb∇aWab, (3.17)

where c∇fa := c∇f(xa)and ∇aWab := ∇W(xa−xb, h). Although the gradient

approximation has been given by (3.17), special formulations closely related to hydrodynamics are more often used in SPH. They are given below.

SPH Rule II. With (3.11) and (3.13), we have ⟨∇f⟩ (x) = ∫ Ω f(x′)∇W(x − x′, h)dx′− f(x) ∫ Ω ∇W(x − x′_{, h)dx}′_, _(3.18)

(38)

which in particle approximation form is c ∇fa= ∑ b mb ρb ( fb− fa ) ∇aWab. (3.19)

Remark 3.4. The SPH rule II is a partition of nullity (PN) (this means that the sum of weights of particle a and its neighbours is zero), while the original SPH gradient (3.17) is not. The SPH gradient (3.19) is exact for constant functions, but the original (3.17) is not.

Similarly, another SPH gradient is written as [112, 150]

∇fa= ∑ b mb ρb ( fb+ fa ) ∇aWab. (3.20)

Although this gradient approximation is symmetric, it is not a PN. When it was applied to the 2D heat equation, the numerical error was much larger than that given by its counterpart (3.19) [30].

SPH Rule III. As a primary variable, mass density is an important quantity in standard SPH. To obtain higher accuracy, the triviality [117]

ρ∇f = ∇(ρf)− f∇ρ (3.21)

is often used. Then, applying (3.17) to the gradients in (3.21) yields

( dρ∇f)_a=∑ b mb ( fb− fa ) ∇aWab. (3.22)

SPH Rule IV. When the gradient of a scalar field is written as ∇f ρ = f ρ2∇ρ + ∇ (_f ρ ) , (3.23)

one can construct another SPH gradient:

(d_∇f ρ ) a= ∑ b mb (_f a ρ2 a + fb ρ2 b ) ∇aWab. (3.24)

The SPH gradient (3.24) can also be derived naturally from a variational principle as shown in [112, 177].

SPH Rule V. Many more alternative formulations of the SPH gradient are

possi-ble. A general identity for ρ∇f can be written as

ρ∇f = ρ2−σ ( f∇ρσ−1+ ρ2σ−2∇(fρ1−σ)), (3.25) with SPH equivalent ( dρ∇f)_a= −ρ2−σa ∑ b mb ρ2−σ_b ( fa− fb ) ∇aWab, (3.26)

(39)

where σ is a constant. When both sides of (3.25) is divided by ρ2_{, we obtain} ∇f ρ = f ρσ∇ ( ₁ ρ1−σ ) + 1 ρ2−σ∇ ( _f ρσ−1 ) , (3.27)

and its SPH equivalent is

(d_∇f ρ ) a =∑ b mb ( _f a ρσ aρ2−σb + fb ρσ bρ2−σa ) ∇aWab. (3.28)

With respect to σ = 2, SPH gradient (3.26) reduces to (3.22), and (3.28) reduces to (3.24). Several gradient formulations with different values of σ such as σ = 1/2 and 3/2 are highlighted for problems with large density gradients [177]. The case with σ = 1 in (3.28) is (d_∇f ρ ) a =∑ b mb (_f a+ fb ρaρb ) ∇aWab. (3.29)

More generalized SPH approximations of ρ∇f and ∇f/ρ can be derived (see

e.g. [178]), which is the reason why so many different variants of SPH formula-tions have been proposed [112, 117, 177]. These variants generally have different conservation properties and calculation efficiencies.

The above approximation rules hold for gradient of a scalar field is valid for the divergence of a vector field. From SPH rule II, we have

d ∇ · fa= − ∑ b mb ρb fab· ∇aWab. (3.30)

The notation fab := fa−fbwill be used for vector quantities throughout this

chapter. From SPH rule III, we get (\ρ∇ · f)_a= −∑

b

mbfab· ∇aWab. (3.31)

Divergence approximations from other SPH rules can be similarly derived. The SPH approximation of the curl of a vector field can be found in [178].

3.3.3 Second derivatives

A straightforward SPH approximation uses the second derivative of the kernel, but this proves to be very noisy and sensitive to particle disorder (more details can be found in e.g. [178]). It is better to approximate the second derivatives utilising only the first derivative of the kernel [151, 178]. Various forms are proposed and used for discretisation of the Laplacian operator in an SPH context [47, 107, 151,

(40)

158,160,178,191,224]. The Laplacian operator may be estimated using the integral approximation (see e.g. [151, 177] for the derivation)

⟨⟨∆f⟩⟩ (x) = 2 ∫ Ω ( f(x) − f(x′))r · ∇W ||r||2 dx′, (3.32)

giving the SPH Laplacian

(c∆f)_a= 2∑ b mb ρb ( fa− fb )rab· ∇aWab r2 ab , (3.33)

where r := x−x′, rab:=xa−xband rab :=||rab|| (Euclidean distance between two

particles a and b). This formulism has been used for heat conduction problems where a second-order parabolic PDE is directly solved (see e.g. [41, 95]), but it is not needed when a first-order system is solved as in [92].

3.3.4 Kernels and their properties

The kernel function W(x − x′, h)determines the accuracy of function

approxima-tion (3.5), while the kernel gradient determines the approximaapproxima-tion accuracy of the first and second derivatives given in Sections 3.4.2 and 3.4.3. To achieve the required accuracy, the kernel and its gradient have to be chosen properly. Be-fore going further, several kernels commonly used in SPH are summarized here, because they are a key element in SPH. The properties will be demonstrated by taking the widely used cubic spline kernel as an example. It is believed that this will be helpful in understanding the SPH approximations and related topics dis-cussed later in this chapter.

The kernel function Wab := W(xa−xb, h)can be written in a general manner as

Wab:= W(rab, h) =

1

hdw (q) , (3.34)

where d is the dimension of the system and q := rab/h.

Standard kernels. Since the kernel is supposed to mimic the Dirac delta distribu-tion, the first option that comes to one’s mind could be the Gaussian

W(q, h) := λ

hde

−q2_, _(3.35)

where the normalization factor λ = π−d/2_{. The Gaussian has the advantage that the}

spatial derivative is infinitely smooth and therefore exhibits good stability prop-erties. The inventors of SPH, Gingold & Monaghan (1977) [69], used the Gaussian kernel in their original work on a three-dimensional astrophysical problem. For practical applications involving boundaries, however, a Gaussian kernel is rarely used these days because it is not compactly supported (the approximation spans

Simulating unsteady conduit flows with smoothed particle hydrodynamics

Simulating unsteady conduit flows with smoothed particle

hydrodynamics

Simulating Unsteady Conduit Flows

with

Simulating Unsteady Conduit Flows

with

Smoothed Particle Hydrodynamics

PROEFSCHRIFT

Summary

Samenvatting

Contents

Introduction

1.1

Motivation for the study

1.2

Problem statement and method justification

1.3

Outline of the thesis

Mathematical Modelling

2.1

Introduction

2.2

Basic concepts and assumptions

2.2.1

Infinitesimal fluid element moving with the flow

2.2.2

The material derivative

2.2.3

The divergence of the velocity

2.3

The continuity equation

2.4

The momentum equation

2.5

Summary of the governing equations

2.5.1

The Navier-Stokes equations for viscous flow

2.5.2

The Euler equations for inviscid flow

2.5.3

Equation of state

2.6

Initial and boundary conditions

2.7

Moving fluid domain

Smoothed Particle

Hydrodynamics (SPH)

3.1

Introduction

3.2

Literature review on SPH

3.3

SPH approximations

3.3.1

Function approximation

3.3.2

Gradient and divergence

3.3.3

Second derivatives

3.3.4

Kernels and their properties