125,99 €
Mehr erfahren.
- Herausgeber: John Wiley & Sons
- Kategorie: Wissenschaft und neue Technologien
- Sprache: Englisch
Unveiling the Secrets of Plasma Physics: A Practical Guide to Computational Simulations Plasma physics focuses on the most abundant state of matter in the universe, corresponding to ionized gas comprising ions and electrons. It can be created artificially and has a huge range of technological applications, from television displays to fusion energy research. Every application of plasma technology requires its own numerical solution to the complex physical and mathematical equations which govern the research field of plasma physics. Modelling and Simulation in Plasma Physics for Physicists and Mathematics offers an introduction to the principles of simulating plasma physics applications. It provides knowledge not only of the fundamental algorithms in computational fluid mechanics, but also their specific role in a plasma physics context. In addition, the book dissects the challenges and advancements, unveiling the delicate balance between accuracy and computational cost. Modelling and Simulation in Plasma Physics for Physicists and Mathematics readers will also find: * Cutting-edge computational insights where powerful simulations meet theoretical complexities, providing physicists and mathematicians a gateway to cutting-edge research. * An overview of programming language-agnostic code generation and the construction of adaptable models that resonate with the intricate dynamics of plasma physics, ensuring precision in every simulation. * Advanced simplification strategies, including time splitting, analytic models, averaged rates, and tabular material, offering scientists and engineers a roadmap to balance computational demands with scientific rigor. Modelling and Simulation in Plasma Physics for Physicists and Mathematics is ideal for plasma physicists, students, and engineers looking to work with plasma technologies.
Sie lesen das E-Book in den Legimi-Apps auf:
Seitenzahl: 437
Veröffentlichungsjahr: 2024
Ähnliche
Table of Contents
Cover
Table of Contents
Title Page
Copyright
Preface
Note
Acknowledgements
Preamble
References
Note
Part I: Continuum Methods
1 Foundations of Computational Fluid Mechanics
1.1 Basic Concepts of Fluid Mechanics
1.2 The Basic Equations of Fluid Mechanics
1.3 Ideal (Dissipationless) Flow – Hyperbolic Equations
1.4 Formal Solution
1.5 Discontinuities
1.6 Plasma Fluid Dynamics
1.7 Basic Principles of Finite Differencing
1.8 Conservative Finite‐Difference Approximations
1.9 Numerical Fluid Approximations
1.10 Grid Geometry
1.11 Control Volume Differencing
1.12 Mesh Types
References
Notes
2 Analytic and Quasi‐Analytic Approximations
2.1 Analytic and Quasi‐Analytic Methods
Appendix 2.A The Nemchinov Conjecture
Appendix 2.B The Energy Integral
References
Notes
3 Numerical Fluid Dynamics in the Eulerian Scheme
3.1 An Introduction to Steady‐State Engineering Design Models
3.2 Spatial Differencing
3.3 Generalised Euler Schemes
Appendix 3.A Energy Conservation in Fluid Codes
References
4 Lagrangian Systems
4.1 Lagrangian Fluid Dynamics
4.2 One‐Dimensional von Neumann–Richtmyer Algorithm
4.3 Multidimensional Lagrangian Schemes
4.4 Choice of Method
References
5 Arbitrary Lagrangian–Eulerian Schemes
5.1 Introduction
5.2 Step 1: The Lagrangian Stage
5.3 Step 2: The Iteration Stage
5.4 Step 3: Mesh Generation
5.5 Step 4: Re‐zoning
References
6 Hybrid or d Schemes
6.1 Introduction
References
7 Magneto‐Hydrodynamics
7.1 Introduction
7.2 The MHD Equations
7.3 MHD Equations
7.4 Conjugate Conservative Operators
7.5 Finite‐Difference Approximations
7.6 Cylindrical Geometry – Self‐Generated Magnetic Fields
References
Notes
Part II: Particle Methods
8 Particle in Cell Simulations
8.1 Introduction
Appendix 8.A Buneman‐Boris Algorithm
Appendix 8.B Stability of the 1D Electrostatic Model
References
Notes
Part III: Stochastic Methods
9 Monte Carlo Schemes
9.1 Monte Carlo Methods
9.2 Monte Carlo Integration
9.3 Random Walks
9.4 Nuclear Reactor Criticality
9.5 Thermodynamic Properties and Equation of State
Supplementary Material
Appendix 9.A Deviates to Specific Distributions
Appendix 9.B Kinematics of Elastic Scattering
References
Notes
10 Numerical Diffusion Schemes
10.1 Introduction
10.2 Split‐Time‐Step and ADI Methods for Solving Diffusion Problems in Orthogonal Cartesian Grid Systems
10.3 The Diffusion Matrix
Appendix 10.A Thomas Algorithm ‐ Tri-diagonal Matrix EquationSolver
References
Notes
11 Ion–Electron Equilibration
Reference
12 Ionisation–Recombination Models
12.1 Introduction
12.2 Collisional‐Radiative Model
12.3 Two‐Stage Model
Supplementary Material
Appendix 12.A Solution of the Collisional-Radiative Equations
Appendix 12.B A Theorem on Determinants
Appendix 12.C An Algorithm Using the Exact Solution
M.1 Quasi‐Linear Equations
M.2 Stiff Equations
M.3 Weak Solutions
M.4 Operator Splitting
M.5 Statistics Primer
M.6 Numerical Solution of Poisson's Equation
M.7 Compressible Gas Potential Flow
M.8 Viscous Incompressible Flow
M.9 Equation of State
References
Notes
Index
End User License Agreement
List of Tables
Chapter 2
Table 2.1 Matching parameters.
Chapter 12
Table M.6.1 Number of unknowns required per cycle.
List of Illustrations
Preamble
Figure I Sketch to illustrate the key interactions in the simulation of a la...
Chapter 1
Figure 1.1 Simple mesh element in cylindrical geometry. The cell is centred ...
Chapter 3
Figure 3.1 Arrangement of a four‐step multi‐grid stage with a sequence of Eu...
Figure 3.2 Representative forms of the standard grid configurations used for...
Figure 3.3 Schematic representation of a block structured grid around an air...
Figure 3.4 Advection of a propagating square wave pulse of height 1 unit wit...
Figure 3.5 Advection of a propagating square wave pulse of height 1 unit wit...
Figure 3.6 Advection of a square pulse of height 1 unit with a Courant numbe...
Figure 3.7 Advection of a square wave pulse of height 1 unit with a Courant ...
Figure 3.8 Advection of a square pulse of height 1 unit with a Courant numbe...
Figure 3.9 Advection of a square pulse of height 1 unit with a Courant numbe...
Figure 3.10 Advection of a square pulse of height 1 unit with a Courant numb...
Figure 3.11 Energy spectrum due to non‐linear hydrodynamic coupling.
Figure 3.12 The development of a Mach 4 shock wave in a polytropic gas of in...
Chapter 4
Figure 4.1 Artificial viscosity stabilised shock profile for gas of polytrop...
Figure 4.2 The shock wave test calculated with the von Neumann–Richtmyer Lag...
Figure 4.3 Basic hourglass mode.
Figure 4.4 Creation of a ‘bowtie’ cell formed by mesh point overtaking.
Figure 4.5 Relationship between the laboratory ( and the Lagrangian ) co‐o...
Figure 4.6 General integration surface within the neighbouring cells surro...
Figure 4.7 Sub‐cell mesh with rectangular sub‐zones in () space.
Figure 4.8 Sub‐cell mesh with diamond sub‐zones in () space.
Figure 4.9 Components of the artificial viscosity produced into response to ...
Figure 4.10 Response of artificial viscosity applied to long, thin zones.
Figure 4.11 Tensor viscosity components for each side of the cell prevents t...
Chapter 8
Figure 8.1 Staggered grid for calculation.
Figure 8.2 Sketch to show the two time level mesh employed for integrating t...
Figure 8.3 Spatial arrangement of variables on a two‐dimensional grid.
Chapter 9
Figure 9.1 The particle following elastic scattering at leaving at an angl...
Figure 9.2 Sketch of the Lennard‐Jones interparticle potential.
Figure 9.3 Periodic grid with randomly positioned molecules.
Figure 9.A.1 Plot of the Watt spectrum and its composite comparison function...
Figure 9.A.2 Comparisons of rejection method Maxwellian distribution functio...
Chapter 10
Figure 10.1 Flux across cell boundaries in one‐dimensional diffusion on a on...
Chapter 12
Figure 12.1 Diagram to show how the total set is the union over all the io...
Figure M.6.1 The progressive reduction of the number of equations and unknow...
Figure M.7.1 Diagram illustrating the relationship of the points and on ...
Guide
Cover
Table of Contents
Title Page
Copyright
Preface
Acknowledgements
Begin Reading
Index
End User License Agreement
Pages
iv
xiii
xiv
xv
xvii
1
2
3
5
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
105
106
107
108
109
110
111
112
113
114
115
116
117
119
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
251
252
253
254
255
256
257
258
259
260
261
262
263
264
Modelling and Simulation in Plasma Physics for Physicists and Mathematicians
Geoffrey J. Pert
University of York, Heslington, YorkUnited Kingdom
This edition first published 2024
© 2024 John Wiley & Sons Ltd
All rights 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, except as permitted by law. Advice on how to obtain permission to reuse material from this title is available at http://www.wiley.com/go/permissions.
The right of Geoffrey J. Pert to be identified as the author of this work has been asserted in accordance with law.
Registered Offices
John Wiley & Sons, Inc., 111 River Street, Hoboken, NJ 07030, USA
John Wiley & Sons Ltd, The Atrium, Southern Gate, Chichester, West Sussex, PO19 8SQ, UK
For details of our global editorial offices, customer services, and more information about Wiley products visit us at www.wiley.com.
Wiley also publishes its books in a variety of electronic formats and by print‐on‐demand. Some content that appears in standard print versions of this book may not be available in other formats.
Trademarks: Wiley and the Wiley logo are trademarks or registered trademarks of John Wiley & Sons, Inc. and/or its affiliates in the United States and other countries and may not be used without written permission. All other trademarks are the property of their respective owners. John Wiley & Sons, Inc. is not associated with any product or vendor mentioned in this book.
Limit of Liability/Disclaimer of Warranty
While the publisher and authors have used their best efforts in preparing this work, they make no representations or warranties with respect to the accuracy or completeness of the contents of this work and specifically disclaim all warranties, including without limitation any implied warranties of merchantability or fitness for a particular purpose. No warranty may be created or extended by sales representatives, written sales materials or promotional statements for this work. This work is sold with the understanding that the publisher is not engaged in rendering professional services. The advice and strategies contained herein may not be suitable for your situation. You should consult with a specialist where appropriate. The fact that an organization, website, or product is referred to in this work as a citation and/or potential source of further information does not mean that the publisher and authors endorse the information or services the organization, website, or product may provide or recommendations it may make. Further, readers should be aware that websites listed in this work may have changed or disappeared between when this work was written and when it is read. Neither the publisher nor authors shall be liable for any loss of profit or any other commercial damages, including but not limited to special, incidental, consequential, or other damages.
Library of Congress Cataloging‐in‐Publication Data
Names: Pert, G. J. (Geoffrey J.), author.
Title: Modelling and simulation in plasma physics for physicists and mathematicians / Geoffrey J. Pert, University of York, Heslington, York, United Kingdom.
Description: Hoboken, NJ : Wiley, 2024. | Includes bibliographical references and index.
Identifiers: LCCN 2024008761 (print) | LCCN 2024008762 (ebook) | ISBN 9781394239207 (hardback) | ISBN 9781394239221 (adobe pdf) | ISBN 9781394239214 (epub)
Subjects: LCSH: Plasma (Ionized gases) | Plasma (Ionized gases)–Mathematical models. | Plasma (Ionized gases)–Computer simulation.
Classification: LCC QC718 .P3755 2021 (print) | LCC QC718 (ebook) | DDC 530.4/4011–dc23/eng/20240324
LC record available at https://lccn.loc.gov/2024008761
LC ebook record available at https://lccn.loc.gov/2024008762
Cover Design: Wiley
Cover Image: © BokehStore/Shutterstock
Preface
This volume is the fruit of a career embracing the development of computational modelling from the 1960s to the present, mainly involving plasma physics. The author has been active as both user and developer of sophisticated computer models for plasma over the past 50 years, starting in the late 1960s. This book is intended as a continuation of our earlier study of the fundamental principles underlying the present understanding of the physics of plasmas and ionised gases in laboratory environments. The author started to study plasma using computational models over 50 years ago and has maintained an activity at the forefront of the field.
The theoretical development of the behaviour of such systems involves complex sets of equations whose calculation is beyond current analytic methods. As a consequence, the theory has increasingly relied on computational modelling methods. These generally involve complex programmes with many interacting elements. The development of these codes has relied heavily on some standard procedures, whose foundations were developed in the early days of computational techniques by pioneers such as von Neumann, Richtmyer, Metropolis and Harlow principally at Los Alamos. Building on these methods has resulted in effective methods for solving the partial differential equations involved in studying the dynamics of plasmas, such as computational fluid dynamics, Monte‐Carlo, finite‐difference methods. In this text, we review these methods applied to the multiple problems involved in solving the behaviour of plasma systems.
At low densities, typical of most terrestrial plasma, the space charge field between the positively charged ions and negatively charged electrons is sufficiently strong at quite small separations that the two ionic species are ‘tied’ together, and the plasma behaves as a single fluid. In this case, the fluid has the particularly simple equation of state of a monatomic polytropic gas, with polytropic index . In the absence of external magnetic fields, and when the ion and electron densities and temperatures are equal, the plasma dynamics may be treated as a standard gas by the well‐established routines of compressible gas dynamics. The early sections see Chapter 1 of this book are therefore devoted to a review of the development of computational gas dynamics. The distinctive application to magneto‐hydrodynamics is considered later. We note however that energy transfer and consequent equilibration is much slower between electrons and heavy particles and that consequently the energy balance of each particle must be treated separately.
The modelling described above is known as the continuum approximation, which yields a method of sampling the pseudo‐distribution function of the computed plasma. The output is the local thermodynamic state and fluid mean velocity.
Alternatively, we may sample the spatial and temporal velocity of a representative group of particles by calculating the dynamics of a group and the electromagnetic fields associated with them. Such an approach is known as a particle method, of which particle‐in‐cell is the most important.
The development of small, inexpensive yet powerful personal computers and laptops has allowed modelling simulations of complex interactive systems at a scale that a decade ago required the use of large multi‐processor vector and parallel processing computer machines. Correspondingly, the expected accuracy, number and complexity of different elements treated in a single computer run have greatly increased. The laboratory facilities and manpower required to run such codes have increased and the costs correspondingly have grown markedly. Thus, we have the situation where we can perform the most accurate calculations which can be compared directly with experiment but at a high price and as a result, only at a limited number of facilities. Alternatively, we can relax the accuracy and perform representative calculations using desktop computers. The enormous increase in CPU power and memory has allowed the most complex calculations that were possible in the 1960s and 1970s to be performed on top‐of‐the‐range desktop and laptop machines. As a consequence, many of the routines developed in the earliest days of computational physics can now be used as the foundation for practical desktop calculations. A very useful reference for these programmes is the set of annual texts ‘Methods in Computational Physics’.1 It is this latter requirement that this book will address. We will not consider the most sophisticated models or the construction of an accurate picture, but only the elements which tied together enable us to easily run and manipulate programmes, thereby allowing us to quickly identify the basic scaling of potential plasma configurations, i.e. as a design tool.
Like many aspects of kinetic physics, there are simplifications that make theoretical studies tractable, for example the markedly different time scales in the Bogobuliev treatment of gas kinetic theory and the Chapman–Enskog approximation in gas transport theory (§4.5 4.11). In this case, it is the different scales relative to the Debye length, where scale lengths are short compared to the Debye length or plasma frequency.
It is a requirement for the existence of the electric field cancellation associated with the Debye length that the plasma be dilute, i.e. the number of electrons contained within the Debye sphere is large. This requirement ensures that the electron distribution in the neighbourhood of any charge imbalance is Maxwellian. The Debye length is then the distance over which the field exists given by
We define a number of parameters which characterise the condition of the plasma from the Debye length.
The Debye parameter is defined by
The plasma or Spitzer parameter by The Spitzer parameter plays an important role as the cut‐off for long‐range weak collisions. It differentiates between the effects of single strong collisions and weak multiple ones – weak multiple collisions being times more effective than strong.
Characteristic times are defined by the plasma frequency
The plasma frequency and the Debye length are related by an average electron velocity.
Laboratory plasma physics in almost all experimental configurations involves dilute plasma with a relatively large plasma parameter () whose values greatly limit the range of practical computational complexity. Values outside this range tend to be limited to astronomical plasma.
The short Debye length in experimental situations has another important effect, namely since electric fields are screened out by the electron mobility it follows from Poisson's equation that the space charge , i.e. that the plasma is quasi‐neutral. This has a very important consequence on the plasma dynamics, namely the plasma behaves as a single fluid, with the electrons and ion tied together by their mutual electric field. We may therefore treat the simulation of the plasma using a single fluid model with two separate energy relations if the electron and ion temperatures are unequal. We will therefore examine in detail the computer models for fluid simulation.
The quasi‐neutral behaviour of the plasma results in the density of ions and electrons being equal. However, it does not prevent the mobile electrons from carrying current and generating magnetic fields, behaviours taken into account in the magneto‐hydrodynamics relationships.
It follows from the small value of the Debye length that the plasma (or Spitzer) parameter is small. Consequently, the fraction of weak multi‐particle collisions compared to strong single interactions is small. This result ensures the effectiveness of multiple weak particle interactions in cell models, where the interaction occurs via the self‐generated fields rather than direct collisions.
Note
1
Methods in Computational Physics, edited by B. Alder, S. Fernbach and M. Rotenberg: published by Academic Press, New York, 1963–1977.
Acknowledgements
I am grateful to several people who have contributed to the work reported in this book:
Firstly, my many graduate students and research fellows who have acted as ‘guinea pigs’ for many of the methods and routines described;
Secondly, my editorial team at Wiley who saved the day when my health failed;
Finally, my colleague and friend Prof G. J. Tallents for enabling me to maintain a continued interaction with the real world of experimental physics.
Preamble
Plasma is defined as a conducting fluid. As such, it is an extensive state of matter within the universe, generally consisting of a fully ionised gas. It is found in many situations both cosmically and terrestrially. In this work, we will be principally concerned with only terrestrial applications, related to controlled fusion research. This encompasses magnetically and inertially confined plasma, together with magnetospheric and applied (engineering) applications.
Experimental studies of plasma in these regimes involves large and complex apparatus, such as tokamaks and high‐power lasers. As a result, theoretical studies play an important role. Again, since the plasma is a material state with a large number of degrees of freedom, theoretical study generally falls into one of two methodologies:
Analytic studies
: These are relatively simple models which can be described by a limited number of parameters and involve a restricted number of coupled equations. They generally describe fundamental activity within the plasma, for example Landau damping.
Computer simulation
: These are more complex models attempting to give an accurate picture of the behaviour of plasma in an experimental environment. The large number of coupled equations demands the application of substantial computer power and storage. Originally this required top‐of‐the‐range machines, but increasingly simple desktop personal computers have had sufficient power to perform substantial realistic simulations. A major advance has been the development of parallel processing systems, essentially a series of coupled PCs operating simultaneously in parallel.
The nature of the theoretical model is determined by the environment in which it is applied. This is characterised by the scale of the system with respect to the Debye length . Again the models fall into one of two classes:
Debye length small (). In this case, the inter‐particle electrostatic forces dominate the motion of the charged particles, whose overall motion is, in turn, dominated by their mutual coupling. This is the region of collective motions characterised by plasma waves. In the context of plasma theory, it is a regime of individual single‐particle motion generally known as
particle‐in‐cell
. At this level of approximation, the self‐consistent inter‐particle fields are treated directly through Maxwell's equations.
This level of approximation corresponds to dilute plasma where the plasma parameter (namely the ratio of the Debye length to electron mean separation) and is the normal plasma state [see (Pert 2021; §1.3.2)].
Debye length large (). The space charge fields established by even a relatively small, spatial charge separation are so large that no space charge density can exist (). The negative electrons and positive ions therefore tend to move as a single unit with mass approximately equal to the ion mass (), i.e. the plasma dynamics are that of a single fluid. In this regime, the plasma is described by standard fluid mechanics, and the simulation follows the methods of computational physics fluid dynamics (CFD).
Under this level of approximation, the basic equations are supplemented by two further conditions.
If the characteristic time scale of the plasma development is long compared to the mean free time (ion/electron collision time), the heavy particles (ions) and electrons are in thermal equilibrium. The particle energy is equilibrated and the energy is calculated by a single first law of thermodynamics relation –
single fluid model
.
On the other hand, if equilibration is slow, the energy relation for the heavy particles and ions must be calculated separately with an equalisation term between the two –
two fluid model
.
Finally, if there are magnetic fields produced by plasma currents or by external fields or pressure gradients, transport coefficients and energy relations are modified. The magnetic fields can be calculated from the Faraday and Biot–Savart laws. When magnetic fields are important, the magneto‐hydrodynamic approximation (MHD) may be valid.
By the early 1960s, the basic equations describing the motion of plasma had been well established and it was clear that in the fluid regime, plasma acted as a continuum and could be described by modified fluid dynamics equations. Under these conditions, computational methods based on CFD were appropriate. By the end of the 1960s to the early 1970s these were well established.
On the other hand, in the particle regime, the dynamic equations of motion are simply established based on Newton's laws of motion and Maxwell's laws of electromagnetism. The scale of modelling was however limited by the speed and storage of the basic computer.
The development of small, inexpensive yet powerful personal computers and laptops has allowed modelling simulations of complex interactive systems at a scale that a decade ago required the use of large multi‐processor vector and parallel processing computer machines. Correspondingly, the accuracy, number and complexity of different elements treated in a single computer run have greatly increased. The laboratory facilities and manpower required to run such machines have increased and the costs correspondingly have grown markedly. Thus, we have the situation where we can perform the most accurate calculations which can be compared directly with experiment but at a high price and as a result only at a limited number of facilities. Alternatively, we can relax the accuracy and perform representative calculations using desktop computers. The enormous increase in CPU power and memory has allowed the most complex calculations that were possible in the 1960s and 1970s to be performed on top‐of‐the‐range desktop and laptop machines. As a consequence, many of the routines developed in the earliest days of computational physics can now be used as the foundation for practical desktop calculations. A very useful reference for these programmes is the set of annual texts ‘Methods in Computational Physics’.1 It is this latter requirement that this book will address. We will not consider the most sophisticated models, which generally require a high‐end parallel machine for the construction of an accurate picture, but only the elements which tied together enable us to easily run and manipulate programmes, thereby allowing us to quickly identify the basic scaling of potential plasma configurations i.e. as a design tool. As a particular case, it is convenient to imagine the plasma formed by the interaction of a powerful laser beam with a solid target (see Pert, 2013). Many of the elements and the interactions involved in such a model are identified schematically in Figure I.
Figure I Sketch to illustrate the key interactions in the simulation of a laser–plasma interaction indicating how different physical processes interact.
References
Alder, B, S Fernbach and M Rotenberg. 1963–1977.
Methods in Computational Physics
. New York: Academic Press.
Courant, R, K Friedrichs and H Lewy. 1967. “On the partial difference equations of mathematical physics.”
IBM J Res Dev
11(2):215–234.
Pert, G J. 2013.
Introduction to Fluid Mechanics
. London: Wiley.
Pert, G J. 2021.
Introduction to Plasma Physics
. London: Wiley.
Richtmyer, R D and K W Morton. 1967.
Difference Methods for Initial Value Problems
. New York: Interscience.
Note
1
Alder, Fernbach and Rotenberg (
1963
–
1977
).
Part IContinuum Methods
1Foundations of Computational Fluid Mechanics
1.1 Basic Concepts of Fluid Mechanics
Computational fluid dynamics attempts to solve the dynamical structures of systems within the fluid approximation, namely when the behaviour of the fluid is dominated by collisions between the constituent molecules. This implies that the mean free path is everywhere smaller than the characteristic dimensions of the problem, and that the mean collision interval short compared to the characteristic time so that the system may be treated as a continuum and individual molecular motions are treated by averaging. These averaging techniques may be applied to plasma when the Debye length is small and the averaging leads to the magneto‐hydrodynamic (MHD) approximation (Pert, 2021). The derivation and validity of the dynamical MHD equations are developed in Pert (2021). Only plasma within the MHD approximation will be treated here. The collision‐free model of plasma are fully treated in Section 8.1 and in more detail by Birdsall and Langdon (1991) and Hockney and Eastwood (1988).
The key underlying equation of motion in the MHD equations are those of compressible fluid mechanics: MHD requiring the introduction of additional effects. Making use of the concept of time splitting discussed earlier these additional behaviours may be treated by separate additional routines. The key fluid routines thus form a basic core to which additional physics is tied. It is therefore appropriate to examine with some care the treatment of the compressible fluid initially and subsequently the additional terms.
Compressible fluid is described by its thermodynamic state and by the velocity of flow. These values must be specified at each point within the flow field. The thermodynamic state is specified by two state variables (e.g. pressure and density ), the remainder being then known from the equation of state. The velocity is defined by the three‐vector components . There are thus five independent quantities, which vary in space and time , and must be determined subject to the constraints on the system. The equations of fluid mechanics contain five equations, which are simply the basic conservation laws of physics applied to moving fluid, whose derivation is well known (Pert, 2013).
Typical problems investigated using these basic fluid techniques include aerodynamics (e.g. aeronautical design or car design), meteorology, oceanography and plasma dynamics.
Let us now examine the characteristic structure of compressible fluid mechanics which, as we shall see, are hyperbolic in form, which leads to problems in formulating numerical approximations.
1.2 The Basic Equations of Fluid Mechanics
1.2.1 Eulerian Frame
The simplest co‐ordinate system, and in many cases the simplest, in which to define the flow of fluid is the stationary one which is relative to the laboratory, known as the Eulerian frame.
In the laboratory (or Eulerian) frame, the general equations of fluid mechanics can be written in several forms, of which the general conservation form is most appropriate for a general discussion (Pert, 2013). In Cartesian tensor form, these are:
where is the specific internal energy of the fluid, the specific enthalpy and the heat flux vector. is the momentum flux tensor
is the total stress tensor
where is the viscous stress tensor and the Kronecker delta (1 if and 0 otherwise).
The heat flow vector may be written in terms of the temperature gradient and the thermal conductivity ()
The temperature is determined from the equation of state in the form .
In the laboratory (or Eulerian) frame, the general equations of fluid mechanics can be written in several forms, of which the general conservation form is most appropriate for a general discussion (Pert, 2013).
where is the specific internal energy of the fluid, the specific enthalpy, and the heat flux vector. is the momentum flux tensor
is the total stress tensor
where is the viscous stress tensor and the Kronecker delta
These five equations contain six unknown quantities, namely , to which we must include an additional term, e.g. entropy to obtain closure. This is normally the equation of state : for example for a polytropic gas , where is the adiabatic index of the gas. The equation of state normally includes coefficients like viscosity and thermal conductivity, which allow the viscous stress and thermal conduction flux to be calculated knowing the spatial gradients of velocity and internal energy.
Neglecting the viscous stress , the momentum balance equation may be cast in a more familiar form, Euler's equation, by subtracting the first equation from the second
where is the appropriate sound speed, (namely Conservation Law form.
The set of equations are of the general conservative form
where is the time derivative moving with fluid, is a general conserved quantity and the flux associated with it.
Clearly, for any volume V enclosed by a bounding surface S,
which expresses the same result in an integral form expressing global balance within the system.
1.2.2 Lagrangian Frame
Alternatively, it is frequently convenient to consider the fluid particles stationary within their own frame, known as the Lagrangian frame. The transformation from one frame to the other is accomplished by the advection operator. Thus, if the particle velocity is , it can be easily seen that if the time derivative in the Lagrangian frame, namely with respect to the stationary particle, is , then the transformation from Lagrangian to Eulerian systems for the general quantity is
which we call the advection equation.
1.3 Ideal (Dissipationless) Flow – Hyperbolic Equations
In many fluids, viscosity and thermal conduction are weak (in particular in gases), i.e. the fluid behaves as an ideal dissipationless continuum. To a first approximation in plasma, viscosity can also be neglected and the medium treated as a two‐component fluid (electrons and ions), but thermal conduction is an important element and dissipational, which may be treated separately within our numerical scheme. We will therefore seek to develop the numerical approximations in an ideal inviscid fluid and subsequently introduce corrections to the flow model necessary to take account of dissipative effects. In these circumstances, the equations of fluid dynamics (1.1) form a set of quasi‐linear hyperbolic equations.
1.3.1 One‐Dimensional Isentropic Flow
One‐dimensional isentropic flow
Let us consider the simple one‐dimensional case of the time‐dependent Eulerian non‐dissipative system, which can be expressed in the general form
where and are general ‘vector sets’ (column matrix) of the fluid variables
and is the matrix . This set of three equations represents the general conservation equation in one dimension, where is explicitly a function of the variables alone. It therefore has the form of three first‐order quasi‐linear differential equations in two variables and .
Since the mathematics is messy if the full set of Eqs. (1.54) are used, it is convenient to reduce them to a simpler form by explicitly forming Euler's equation (1.14) from the mass (1.1a) and momentum conservation (1.1b) equations using the pressure/density relation for the adiabatic sound speed , and taking account of the entropy constancy. We correspondingly replace the energy conservation Eq. (1.1c) by the isentropic condition, such that the entropy of a fluid particle remains constant, Eq. (1.14b).
where is the entropy and is the sound speed squared. For this set
The array is given by substitution
The fluid equations are hyperbolic if the eigenvalues of are real. In this case of an ideal fluid, they can be shown to be , where is the isentropic sound speed: . To find the eigenvalues , solve the determinantal equation
to give . The eigenvectors are
respectively. is the identity matrix: (.
Since the eigenvalues are distinct we may diagonalise by a similarity transformation so that
where is diagonal and is a function of the variable only. Hence, substituting for
for the eigenvalues . The functions , obtained from , are constant on the lines , and are known as Riemann invariants. For this case, they take the form
for the ideal fluid equations.1
In principle, a knowledge of the initial values of these quantities will allow us to integrate along the characteristics in space‐time and determine the flow. Since given the initial values at time zero , we may form an expansion in terms of the Riemann invariants and follow them in time along the characteristics. These results form the basis for much of the study of compressible fluids and by extension of simple plasma neglecting thermal conduction.
1.4 Formal Solution
At any point , we may determine the three values of the variables if the Riemann invariants are known, provided we know the equation of state, i.e. or more succinctly by tracing the development of the fluid variables in space and time along the characteristics. Thus, in principle, we integrate along the field of characteristics to determine the flow. Clearly, the value of at the point is determined only by points which lie within the outermost characteristics through – its domain of dependence. Similarly, it can only influence future events within its outermost outgoing characteristics – range of influence. The existence of such solutions has been formally demonstrated (Courant and Hilbert, 1962).
For this system, the following uniqueness theorem is valid (Courant and Friedrichs, 1948; Courant and Hilbert, 1962).
Consider a solution with continuous second‐order derivatives in the region ABP bounded by two characteristics through P, and the domain of dependence AB cut‐off by them on the initial line I. Suppose another solution with continuous second derivatives is given in ABP, which assumes the same values on AB as the first. Then, the second solution is identical to the first within the domain ABP.
1.5 Discontinuities
The uniqueness theorem quoted earlier requires that the solution be continuous up to the second order everywhere within the domain of dependence. However, discontinuities arise naturally within fluid mechanics by virtue of the non‐linearity of the governing equations (Courant and Friedrichs, 1948; Landau and Lifshitz, 1959). These occur in two forms: strong and weak discontinuities. Strong discontinuities involve a jump in the value of the flow variables and include shocks and contact discontinuities. Shocks involve flow through the discontinuity, whereas tangential discontinuities do not. Weak discontinuities are continuous, but with a change of gradient. Weak discontinuities, which include the head of rarefaction waves, propagate along the characteristics (Landau and Lifshitz, 1959).
If two characteristics were to intersect, the Riemann invariant would become indeterminate. Such intersections can occur in compression, where they lead to shocks. In this neighbourhood, the ideal flow condition breaks down. When the characteristics approach over distances of the order of the mean free path, large gradients in the flow are established, and dissipation becomes important. This leads to an entropy increase in accordance with the second law of thermodynamics across the shock. The thickness of the ‘discontinuity’ is of the order of the mean free path, which is vanishingly small in the continuum theory. In general fluid theory, the collisional behaviour at a molecular level gives rise to dissipational processes, namely viscosity and thermal conduction. Within the fluid description, the shock ‘discontinuity’ is controlled by dissipation, principally viscosity. The shocks therefore appear as narrow zones within which viscous effects are dominant, which in molecular terms are typically a few mean paths thick (Zel'dovich and Raizer, 1966).
Considering the case of a general fluid with dissipation, the general conservation law in integral form (1.10) requires that across a normal shock in the laboratory frame
where and , and , and and and are, respectively, the density, flow velocity, pressure and enthalpy upstream () and downstream () of the shock in the laboratory frame. Alternatively and more succinctly in the shock frame
where is the shock speed, the shock moving from left to right, are the fluxes normal to the shock and and the values of the components of on the upstream (left) and downstream (right) sides of the shock, respectively. Note that this relation must be obeyed for all the elements ofwith the same value of. In one dimension these are the familiar Rankine–Hugoniot relations for a normal shock (Eq. (1.27)).
Discontinuities arise naturally in the general solution of the hyperbolic equations. At such points, the differential equations break down, but these are expressions of more general integral relations, namely the integral forms of the conservation laws. We may construct generalised solutions which satisfy these relations called weak solutions by Lax (1954) (see Supplement Section M.3). Clearly, they satisfy the differential equations where the derivatives exist, and the Rankine–Hugoniot relations across discontinuities. However, in contrast to the actual physical situation, rarefaction as well as compression discontinuities can occur.
Two questions arise if we wish to use the dissipationless equations for calculation.
Is the solution unique? Clearly, our previous uniqueness theorem fails when discontinuities occur.
Is the solution without dissipation the limit of vanishing viscosity to the complete fluid dynamic problem?
In fact both questions are related. The discontinuities in ideal flow can be shown not to be uniquely specified by the Rankine–Hugoniot condition. In the real fluid, such uniqueness is imposed by the additional requirement for an entropy increase. This is established by the following results proved for one‐dimensional systems (Oleinik, 1963a).
Solutions of dissipative equations converge to weak solutions of the hyperbolic system in the limit of vanishing viscosity (Oleinik,
1963a
).
Such solutions of the viscous system satisfy
known as the entropy condition, where is the shock speed (
1.27
) and
where lies between and (Hopf,
1969
; Quinn,
1971
).
Weak solutions of the hyperbolic system satisfying the entropy condition are unique (Oleinik,
1963b
; Quinn,
1971
).
The entropy condition has a close relationship with many of the familiar properties of shocks. Thus, since the limit of a weak shock is a sound wave, which propagates at the local sound speed , the shock speed is supersonic with respect to the upstream flow and subsonic with respect to the downstream flow . Indeed the entropy condition reflects the facts that shock adiabat is generally convex, and that physical shocks are compressive as required by the second law of thermodynamics (Courant and Friedrichs, 1948). The entropy condition thus ensures that the shock satisfies the second law of thermodynamics, and that there is an entropy increase across the shock.
1.6 Plasma Fluid Dynamics
Since the viscosity of plasmas is very small, we may neglect it and the dynamics can be treated as those of an ideal fluid, with the proviso that the ions and electrons may have different temperatures and , respectively. As a consequence of ion/electron collisions, there is an energy exchange between the two sets of particles. The rate at which this takes place is governed by the equilibration time ,
where and are the specific internal energies and specific heats of ions and electrons, respectively, is the ion/electron relaxation time given by Spitzer (1956).
As a result of the different internal energies, the pressures of the ions and electrons, and , respectively, are different; the total pressure follows from Dalton's law . The individual species act as almost ideal gases, so that
Fortunately, provided the Debye length is much smaller than the characteristic lengths in the fluid, the ion and electron charge densities are equal. As a result,
where is the degree of ionisation. If is the fractional ionisation of the th stage of ionisation, the fractional ionisation where is the fractional charge on the ion . The fractional ionisation in stage is calculated from the balance between its loss by ionisation and recombination and gain from the ionisation of stage and recombination of stage , namely
where and are the rates of ionisation and recombination, respectively.
1.7 Basic Principles of Finite Differencing
1.7.1 Introduction
We shall develop numerical solutions for the case of ideal dissipationless flow, which, as we shall see, presents serious difficulties in view of its non‐linear hyperbolic nature. Such solutions are of great importance in both physics and engineering encompassing many of the applicable classical disciplines, namely aerodynamics and plasma physics. Since the central problem is most easily attacked in one dimension, we will firstly present the analysis considering only one spatial variable and the temporal . In general, the exact solution is defined at each point . The numerical solution will take the form of an approximation function , which is only defined at a set of discrete points in . These points conveniently form a mesh in , which in a simple form have a constant separation, in space and in time. We may identify these points by integer co‐ordinates .
Values of at points intermediate to the mesh are constructed by interpolation.
We require a finite‐difference scheme to be constructed to allow the approximating function to be determined from the boundary conditions. These are usually initial value conditions, i.e. we know the initial fluid state and need to calculate its subsequent history subject to boundary conditions in space. Thus, it is natural in many cases to seek to develop the solution in time. Hence, we require a relation of the type
where is the value of at and . The mesh runs over cells from the ‘left’ boundary at to the right at . The scheme is called implicit if the right‐hand side of the equation contains terms evaluated at time , and is explicit otherwise. Implicit schemes generally involve complex matrix inversion and may have other disadvantages for fluid dynamics.
In general, storage limitations usually allow only one time step to be stored, which restricts to a one‐step operator
For an explicit two‐step operator, there is a fundamental time‐step constraint determined by the characteristics of the flow. Since information is carried along the characteristics with a characteristic speed , and since within an explicit scheme information must be carried to the neighbours only, it is seen that the time step must be limited to a maximum value . This is known as the Courant–Friedrichs–Lewy limit following the pioneering paper by Courant, Friedrichs and Lewy (1928) in which it was first derived. The parameter is known as the Courant number.
Since the exact solution is conserved, it is natural to expect that obeys a similar relation, which on this simple mesh requires
where and are approximations to the boundary fluxes of and , respectively.
This requires
where denotes a time difference, and a space difference. is a function of the variables and of .
This integral (conservation) representation of the finite‐difference equation allows another interpretation of , which may be useful, in terms of the volume average of the quantity taken over the cell, , at time .
In order that the approximation be valid as an expression of the exact solution of the quantity , we clearly require that as . Indeed is a weak solution in this limit provided (Lax and Wendroff, 1960).
These solutions do not, however, necessarily satisfy the entropy condition. This is known to be rigorously satisfied only for one class of finite‐difference schemes (namely monotone) Jennings (1974). Commonly used approximations can be made to converge onto the correct weak solution, with the discontinuities matching the shocks found in practice – a procedure known as shock capture.
1.7.2 Necessary Conditions in the Finite‐Difference Operator
It is important to specify some necessary conditions for the general finite‐difference operator to give satisfactory approximations to the solutions of the corresponding partial differential equations. If is such an operator, we may identify the following necessary conditions on (Richtmyer and Morton, 1967).
Consistency:
As the spatial and temporal steps are reduced together to give a finer mesh, the solution of the difference equations must reduce to the exact solution of the corresponding differential equation. Thus, if as and is any genuine solution of the differential operator to which approximates, then is a consistent approximation to if
For a conservative system, this is equivalent to the previous result (Eqs. (1.37) and (1.38)).
Convergence:
Operating times with the finite‐difference operator on should yield an approximation to . If , then the approximate solution converges to for a sequence of time steps : tending to zero, and such that as for all possible initial states . Thus,
for all possible initial states .
Stability:
It is essential that during the calculation no component be amplified without limit. Thus, an approximation is said to be stable, if for some time , and the infinite set of operators
is uniformly bounded for specified and .
Stability and convergence are related by Lax's equivalence theorem (Richtmyer and Morton, 1967).
Given a properly posed initial‐value problem and a finite‐difference approximation to it that satisfies the consistency condition, then stability is the necessary and sufficient condition for convergence.
Accuracy:
For smooth solutions , the accuracy is limited by the truncation of the Taylor's series at any point in generating the finite‐difference form. In the case of constant coefficients, we define the order of the difference scheme by the largest real number for which
for all sufficiently smooth solutions of the differential equation. Usually, is an integer, in which case consistency requires that the equations be at least ‘first‐order accurate’. Clearly, gives a measure of the missing terms from the Taylor expansion of the function in the difference scheme.
Positivity:
An additional physical constraint is positivity, i.e. that some quantities, for example density and energy, are always non‐negative. We may require that the finite‐difference scheme retain this property.
Monotonicity:
Some physical behaviour maintains a monotonic profile. Finite‐differences schemes with this property are known as
monotonicity preserving
.
Closely related is a monotone finite‐difference operator defined by
Monotone operators are always monotonicity preserving (although the converse is not true), and are necessarily first order (Jennings, 1974). They have the important and useful property that they always satisfy the entropy condition, and thus guarantee uniqueness, i.e. they guarantee correct shock capture.
1.7.3 Constant Coefficients
If the matrix is constant, then the governing differential equations are linearised, and the characteristic equation (Eq. (1.12)) reduces to a particularly simple form, which can be solved exactly by Fourier transformation.
As an example, consider the one‐dimensional single‐variable case, the generalisation to many variables and/or dimensions is straightforward. The governing equation is
subject to an initial condition . The formal solution is
corresponding to a transport (advection) of the quantity with velocity .
Alternatively, let be the spatial Fourier transforms of . Then
Since is a linear operator, this solution also carries over into the finite‐difference form
Since is defined only by values at the spatial points , its Fourier transform is discrete with values of in the range .
In this form, the function , with corresponding arguments, is the amplification factor of the Fourier component . Its magnitude, , determines the amplitude growth of the Fourier component with wave number , and the argument, arg(), the phase shift from step to step at a fixed point. As a result, determines the stability of this system. Hence, we deduce that for stability
for all .
This result can be extended to give a very useful necessary condition for stability due to von Neumann (Richtmyer and Morton, 1967), namely .
When this result is generalised to several variables, becomes a matrix, amplification matrix, with eigenvalues . In this case, von Neumann's condition is generalised to
Although the Von Neumann condition is not sufficient, it has in practice proved extremely useful in judging the stability of finite‐difference schemes. When a difference scheme is found to be conditionally stable (i.e. for a restricted range of , for example) the von Neumann condition nearly always gives the correct range. Sufficient conditions for constant coefficients have been derived but are less convenient to use.
In practical problems, the differential equations do not have constant coefficients. However, practical experience shows that instabilities usually start as local phenomena. In a small region, the differential equation can be linearised to a form with constant coefficients, which can be tested for local stability. This is a useful procedure and does give necessary conditions for stability. However, it must be remembered that non‐linear instabilities, which are more difficult to assess, can and do occur.
1.8 Conservative Finite‐Difference Approximations
The basic fluid hydrodynamic (and more generally magneto‐hydrodynamic) are in the form of a conservation law (Eq. (1.57)). For application in numerical simulation, it is natural that these should be reproduced in the finite‐difference form. Thus, consider a closed system (with no boundary fluxes) in which the quantities generated by the numerical procedure at a mesh point at the conclusion of the th time step are to be conserved, i.e. is constant in time. In general, for a closed system
In many cases, the physical quantity is inherently non‐negative, for example density, energy, etc. in which case the finite scheme should be non‐negativity maintaining
As a result of truncation errors introduced thereby, the resulting equations are no longer exact conservation laws. We identify different types:
Strong conservation:
If the finite‐difference equations obey condition (
1.47
) within the round‐off error, then the system is
strongly conservative
.
Theorem 1.1 A strongly conservative, non‐negativity maintaining system is always stable if the input data are non‐negative.
Proof 1.1 The proof follows directly from the definition of stability since any element is bounded from above by
