UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Investigation into the feasibility and operation of a magnetized target fusion reactor : insights from… Lindstrom, Michael 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2015_september_lindstrom_michael.pdf [ 2.23MB ]
JSON: 24-1.0167722.json
JSON-LD: 24-1.0167722-ld.json
RDF/XML (Pretty): 24-1.0167722-rdf.xml
RDF/JSON: 24-1.0167722-rdf.json
Turtle: 24-1.0167722-turtle.txt
N-Triples: 24-1.0167722-rdf-ntriples.txt
Original Record: 24-1.0167722-source.json
Full Text

Full Text

Investigation into the Feasibility andOperation of a Magnetized TargetFusion ReactorInsights from Mathematical ModellingbyMichael LindstromB.Sc., The University of British Columbia, 2008M.Sc., The University of British Columbia, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)May 2015c© Michael Lindstrom 2015AbstractMagnetized target fusion reactors are a modern idea to generate hydrogenfusion energy on earth. The design entails confining a plasma with a mag-netic field and crushing it in an imploding shell of molten metal. Such adesign has many unresolved questions in terms of its feasibility as a powersource and ways to make it efficient. In this thesis, we will look into two ofthe approaches undertaken to explore these questions. Firstly, we use a coor-dinate transformation and implement a novel flux-limited, split-step, finitevolume scheme for nonlinear coupled hyperbolic partial differential equa-tions. With this numerical scheme, we do a parameter sensitivity analysisfor the design performance. Secondly, by a careful series of asymptotic ar-guments, we establish a leading order asymptotic expression for the plasmacompression. This expression is qualitatively consistent with the numericalwork, but it also gives new insights into how the device operates. Togetherthese approaches allow us to infer key design parameters for the success ofmagnetized target fusion. We will conclude with a look into the viability ofmagnetized target fusion and some problems for future work.iiPrefaceRelative ContributionsVarious chapters are derived from peer-reviewed journal articles (onepublished, and one under revision). The numerical finite volume work istaken from “Investigation into Fusion Feasibility of a Magnetized Target Fu-sion Reactor: A Preliminary Numerical Framework” (published with Jour-nal of Fusion Energy, 2014) [26] with co-authors Brian Wetton and SandraBarsky. The asymptotic analysis comes from “Asymptotic Analysis of Mag-netized Target Fusion Reactors” (under revision) [25] where I was the soleauthor. At least 95% of the writing in the finite volume paper was doneby myself; I did all of the numerical simulations, in consultation with theco-authors.Publications• “Investigation into Fusion Feasibility of a Magnetized Target FusionReactor: A Preliminary Numerical Framework” (Journal of FusionEnergy, 2014, published) [26]• “Asymptotic Analysis of Magnetized Target Fusion Reactors” [25]Most of section 1.2, part of chapter 2, and much of chapter 4 are takenfrom [26]. Chapter 5 is taken from a revision of [25].iiiPrefaceA previous version of this manuscript included work on modelling super-conductivity, but this work will appear as its own publication and not partof this thesis.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxI Introduction and Motivation . . . . . . . . . . . . . . . . . 11 Context for This Work . . . . . . . . . . . . . . . . . . . . . . 21.1 Scope of This Thesis . . . . . . . . . . . . . . . . . . . . . . . 21.2 Magnetized Target Fusion . . . . . . . . . . . . . . . . . . . . 31.2.1 Fusion Motivation . . . . . . . . . . . . . . . . . . . . 31.2.2 Magnetized Target Fusion . . . . . . . . . . . . . . . 4II Background to the Model for Magnetized Target Fu-vTable of Contentssion and Finite Volume Methods for Hyperbolic Conserva-tion Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.1 Overall Components . . . . . . . . . . . . . . . . . . . . . . . 92.2 Deriving the Euler Equations . . . . . . . . . . . . . . . . . . 92.2.1 Mass Conservation . . . . . . . . . . . . . . . . . . . 102.2.2 Momentum Conservation . . . . . . . . . . . . . . . . 112.3 Equation of State for Lead-Lithium . . . . . . . . . . . . . . 132.4 Equation of State for Plasma . . . . . . . . . . . . . . . . . . 142.4.1 Plasma Gas Pressure . . . . . . . . . . . . . . . . . . 152.4.2 Plasma Magnetic Pressure . . . . . . . . . . . . . . . 192.4.3 Total Plasma Pressure . . . . . . . . . . . . . . . . . 243 Brief Background into Finite Volume Methods for Conser-vation Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.1 Finite Volume Schemes for Conservation Laws . . . . . . . . 253.1.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . 263.1.2 Finite Volume Setup . . . . . . . . . . . . . . . . . . 263.2 Linear Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 293.2.1 Low Resolution . . . . . . . . . . . . . . . . . . . . . 293.2.2 CFL Condition . . . . . . . . . . . . . . . . . . . . . 333.2.3 High Resolution . . . . . . . . . . . . . . . . . . . . . 333.2.4 Flux Limiters . . . . . . . . . . . . . . . . . . . . . . 343.2.5 Extrapolation . . . . . . . . . . . . . . . . . . . . . . 373.3 Extensions to Nonlinear Systems . . . . . . . . . . . . . . . . 373.3.1 General Idea . . . . . . . . . . . . . . . . . . . . . . . 37viTable of Contents3.3.2 Example of Different Methods as Applied to Burger’sEquation . . . . . . . . . . . . . . . . . . . . . . . . . 393.3.3 Example of the Importance of Conservative Form andProblems with Geometric Source Terms . . . . . . . . 45III Journal Articles . . . . . . . . . . . . . . . . . . . . . . . . 494 Numerical Finite Volume Investigation . . . . . . . . . . . . 504.1 Physical Assumptions . . . . . . . . . . . . . . . . . . . . . . 504.2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . 514.2.1 Parameters, Initializations, and Setup . . . . . . . . . 534.3 Numerical Computations . . . . . . . . . . . . . . . . . . . . 544.3.1 Change of Variables to Map Problem to Fixed Domain 544.3.2 Finite Volume Discretization . . . . . . . . . . . . . . 554.3.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . 574.4 Model Predictions . . . . . . . . . . . . . . . . . . . . . . . . 594.4.1 Plasma Compression and Fusion Possibility . . . . . . 594.4.2 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . 614.5 Summary and Future Work . . . . . . . . . . . . . . . . . . . 635 Formal Asymptotic Analysis . . . . . . . . . . . . . . . . . . . 655.1 Context for the Analysis . . . . . . . . . . . . . . . . . . . . 655.2 Physical Assumptions . . . . . . . . . . . . . . . . . . . . . . 665.3 Equations and Nondimensionalization . . . . . . . . . . . . . 675.3.1 Euler Equations . . . . . . . . . . . . . . . . . . . . . 685.3.2 Boundary Conditions . . . . . . . . . . . . . . . . . . 695.3.3 Initial Conditions . . . . . . . . . . . . . . . . . . . . 70viiTable of Contents5.3.4 Nondimensionalization . . . . . . . . . . . . . . . . . 705.3.5 Selecting an Asymptotic Parameter . . . . . . . . . . 715.4 Formal Asymptotic Analysis . . . . . . . . . . . . . . . . . . 735.4.1 Phase I . . . . . . . . . . . . . . . . . . . . . . . . . . 745.4.2 Phases II and III . . . . . . . . . . . . . . . . . . . . 825.4.3 Phases IV and V . . . . . . . . . . . . . . . . . . . . . 955.5 Comparison with Numerics . . . . . . . . . . . . . . . . . . . 1035.5.1 Pulse Formation . . . . . . . . . . . . . . . . . . . . . 1045.5.2 Focusing . . . . . . . . . . . . . . . . . . . . . . . . . 1055.5.3 Residual Velocity Field . . . . . . . . . . . . . . . . . 1065.5.4 Outer Region Describing Inner Wall Velocity . . . . . 1065.5.5 Qualitative Agreement of Inner Wall Velocity . . . . 1075.5.6 Predictions for Plasma Compression . . . . . . . . . . 1085.6 Conclusions and Future Work . . . . . . . . . . . . . . . . . 111IV Further Details and Conclusions . . . . . . . . . . . . . 1136 Additional Insights Concerning the Numerical and Asymp-totic Modelling of the Fusion Reactor . . . . . . . . . . . . . 1146.1 Fusion Reactor Numerical Details . . . . . . . . . . . . . . . 1146.1.1 Transformed Coordinate System . . . . . . . . . . . . 1146.1.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . 1166.2 Asymptotic Commentary . . . . . . . . . . . . . . . . . . . . 1176.2.1 Different Model than Numerical Simulations . . . . . 1176.2.2 Suitability of Linear Lead-Lithium Equation of State 1176.2.3 Inclusion of Gas Pressure . . . . . . . . . . . . . . . . 1186.2.4 Justification of Equation (6.3) . . . . . . . . . . . . . 120viiiTable of Contents6.2.5 Numerical Validation . . . . . . . . . . . . . . . . . . 1237 Summary and Future Work . . . . . . . . . . . . . . . . . . . 1257.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1257.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133A Detailed Pseudo-Code of Fusion Finite Volume Scheme . 133ixList of Tables2.1 Data for quadratic equation of state. Using these data points,we are able to construct an empirical fit to the equation ofstate for lead. . . . . . . . . . . . . . . . . . . . . . . . . . . . 13xList of Tables3.1 An illustration of the grid variables. The solution and truespatial mesh are defined at indices i = 0, 1, ..., N . By constantextrapolation, we extend the solution and the spatial meshgiving the E-variables with indices i = −2,−1, 0, ..., N,N +1, N+2. By subtracting these E−values and averaging them,we obtain difference values D and interpolated values I bothwith indices −3/2,−1/2, ..., N + 1/2, N + 3/2. With the I−values, we can also interpolate a linearized system and deter-mine the components of the D-values in the coordinates of theright eigenvectors, giving us the w variables. By upwinding,we determine the fluctuation ratios and the limiters θ and φwith indices i = −1/2, 1/2, ..., N − 1/2, N + 1/2, and we canalso define limited w-variables, the w˜-variables. In makinguse of the matrices formed at the I−cells, we can also obtainthe upwind E-values, uˆ, relevant to define the fluxes fˆ . Asecond-order flux limited method would then determine thelow resolution numerical flux, fˆ along with a second-orderflux correction that is suitably limited. . . . . . . . . . . . . . 393.2 The L1-numerical errors are computed using the trapezoidalrule where LF denotes the Lax-Friedrichs scheme, G denotesthe Godunov scheme (first-order), LTD denotes the limitedsecond-order method, and LW denotes the unlimited second-order method (Lax-Wendroff). Based on the slope of the bestfitting line of the log-log plot of the error vs N , we also notethe convergence rate. . . . . . . . . . . . . . . . . . . . . . . . 43xiList of Tables3.3 The L1-numerical errors are computed using the trapezoidalrule where G1 denotes a first-order scheme with a geometricsource term; G2 denotes a non-limited second-order schemewith a geometric source term; GL2 denotes a flux-limitedsecond-order scheme with a geometric source term; C1 de-notes a conservative first-order scheme; C2 denotes a non-limited, conservative second-order method; and CL2 denotesa conservative, flux-limited second-order method. Based onthe slope of the best fitting line of the log-log plot of the errorvs N , we also note the convergence rate. . . . . . . . . . . . . 484.1 A list of the parameters used in the numerical simulations. . 544.2 Convergence of Numerical Scheme at t = 1.243× 10−4 s . . . 594.3 A sensitivity analysis of the reactor performance upon vari-ations in the model parameters. We denote the minimumradius by rm, the maximum total pressure by PM, the com-pression time by tc, the maximum temperature by θM, andthe Lawson triple product by ΠL. The final three rows of thetable illustrate the most promising conditions for fusion. . . . 635.1 Table of physical constants/parameters. . . . . . . . . . . . . 685.2 Characteristic dimensions for the system. The first three en-tries are chosen by the operating conditions. The last twoentires can be derived from the first three. . . . . . . . . . . . 715.3 The dimensionless parameters of the asymptotic model. . . . 73xiiList of Tables5.4 Verifying the numerical validity of the approximation givenin (5.44) by doubling/halving the governing parameters andcomputing I∗ with the approximation of v20 given in (5.44)and the numerically integrated value. . . . . . . . . . . . . . . 945.5 Asymptotic and numerical predictions of minimum radius ofplasma for different values of  with b = 1.05, χ = 0.937 andγ = pi. Based on a least squares linear regression in a log-logplot of the difference versus , the convergence rate appearsto be O(2.3) = o(). The convergence rate is excellent likelyeither due to small coefficients in the o() asymptotic seriesterms or fortuitous cancellations of higher order asymptoticerror terms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1096.1 Asymptotic and numerical predictions of minimum radius ofplasma for different values of  with b = 1.05, χ = 0.937,γ = pi, and µ = 4. The asymptotic minimum radius shouldbe 1.4415. The apparent convergence rate is O(1.7). . . . . . 124xiiiList of Figures1.1 Magnetized Target Fusion Reactor Design . . . . . . . . . . . 72.1 Quadratic fit to experimental equation of state data. . . . . . 142.2 The geometry of our simple sheet of plasma. The currentdensity is λ. The closed loops we consider are labelled, as isthe differential surface element considered for computing thepressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3 The r−direction runs horizontally and the currents run intothe page as indicate by the ⊗. The magnetic fields generatedby the currents are drawn in blue and within the red boxes,the r−components of the field cancel. . . . . . . . . . . . . . 213.1 Diagram of basic finite volume setup. . . . . . . . . . . . . . . 273.2 A sketch of the solution to the Burger’s equation illustratingthe shocks and expansion fans. . . . . . . . . . . . . . . . . . 42xivList of Figures3.3 An illustration of different numerical solvers applied to theinviscid Burger’s equation. The exact solution is plotted ingreen; the first-order Godunov method, based on upwinding,is plotted in red; the first-order Lax-Friedrich’s scheme is plot-ted in dark blue; the second-order flux-limited method withapproximate Riemann solver is plotted in light blue; and thesecond-order Lax-Wendroff method without limiters is plot-ted in black. N = 50 mesh points were chosen. . . . . . . . . 443.4 We observe that the numerical scheme in conservative formis slightly better than the numerical scheme with geometricsource terms. The plots depict the solution at t = 0.5 withN = 25 mesh points. . . . . . . . . . . . . . . . . . . . . . . . 474.1 Spherically Symmetric Model . . . . . . . . . . . . . . . . . . 534.2 Profiles of velocity (left) and density (right) at various timevalues compared to the 1/rs growth rate predicted in theacoustic limit. . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.3 Profiles of velocity (left) and density (right) at t = 7.03 ×10−4s for N = 500 and N = 4000. The plots are very similar. 584.4 The radius of the plasma vs time in the simulations. . . . . . 604.5 The radial velocity of the plasma-lead-lithium interface vstime in the simulations. . . . . . . . . . . . . . . . . . . . . . 604.6 The total plasma pressure vs time in the simulations. . . . . . 615.1 The geometry we consider is spherically symmetric, with rLand rR free boundaries. . . . . . . . . . . . . . . . . . . . . . 67xvList of Figures5.2 Relevant space-time scales and boundary motion. In phase I,pulses are formed, moving towards the plasma. In phase II,the pulses move radially inward, growing in amplitude due tofocusing. In phase III, the pulses interact with the plasmawith much of the energy reflecting but a small portion ofuseful energy remaining. In phase IV, the plasma slowly be-gins to collapse and in phase V, the maximum compression isreached. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.3 Depiction of phase I. This is the formation phase. In thisphase, ρ ∼ 1 + ρ1 + 3/2ρ2 and v ∼ 1/2v0 + v1. . . . . . . . 775.4 Illustration of the right- and left-going characteristic pathsand parameterizations. . . . . . . . . . . . . . . . . . . . . . . 785.5 Plots depicting the profiles of the velocity and density at avery small value of t to validate the asymptotics of phaseI. In the case of the second row, the two term asymptoticexpansions agree so well with the numerics the plots cannot bedistinguished. Parameters: b = 0.557, χ = 1.19, γ = 3.52,  =0.0126, t = /2 (top row); b = χ = γ = 1,  = 0.001, t = /2(bottom row). . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.6 Depiction of phase II. This is the focusing phase. In thisphase, we initially have ρ ∼ 1 + ρ1 + 3/2ρ2 and v ∼ 1/2v0 +v1, and finally obtain ρ ∼ 1 + 1/2ρ1 + ρ2 and v ∼ v0 + 1/2v1. 835.7 The pulses are moving to the left. The pressure peak growslike 1/r as the pulses move inwards. The velocity has a peakwith negative value and one with positive value. The negativevalue grows roughly like 1/r. Parameters: b = χ = γ = 1, = 0.0025,t = 0.0262. . . . . . . . . . . . . . . . . . . . . . . 91xviList of Figures5.8 Depiction of phase III. This is the reflection phase. In thisphase, we start with ρ ∼ 1 + 1/2ρ1 + ρ2 and v ∼ v0 + 1/2v1and finally obtain ρ = 1 + o() and v = O(1/2). . . . . . . . . 925.9 Depiction of phase IV. This is the slow compression phase. Inthis phase, ρ ∼ 1 +O(2) and v = O(1/2). The peak velocityand pressure are growing in this phase. . . . . . . . . . . . . . 975.10 Depiction of phase V. This is the rapid compression phasewhere the plasma reaches its minimum radius. In this phase,ρ ∼ 1 +O(1/2) and v = O(−1/4). . . . . . . . . . . . . . . . 1005.11 Numerical and asymptotic descriptions of the growth in peakvalues for the minimum value of velocity and peak value ofpressure. Based on a best fit, the numerically observed scal-ings are 1/r0.99 and 1/r0.90 for the left and right plot respec-tively. Parameters: b = χ = γ = 1,  = 0.001. . . . . . . . . . 1055.12 A plot of the numerical and asymptotic velocity profile aftermuch of the pulses have reflected off of the inner wall. Thetwo are nearly indistinguishable. Parameters: b = χ = γ = 1, = 0.0025, t = 0.0784. . . . . . . . . . . . . . . . . . . . . . . 1065.13 A plot of the numerical and asymptotic velocity profile forhow the inner wall velocity implicitly depends on the innerwall position. Parameters: b = χ = γ = 1,  = 0.0025, tranging from 0.0784 to 0.11. . . . . . . . . . . . . . . . . . . . 107xviiList of Figures5.14 A numerical validation that the asymptotics have qualita-tively captured the behaviour of the inner wall. Parameters:b = χ = γ = 1,  = 0.0025. The different qualitative regimesare labelled in the diagram. The r∗-value that we seek oc-curs at the bottom spike of the curve where the maximumcompression label appears. . . . . . . . . . . . . . . . . . . . . 1086.1 Plot of quadratic and linear equation of state models withpeak pressure as computed numerically. . . . . . . . . . . . . 118xviiiAcknowledgementsA gigantic thanks to my thesis supervisor Brian Wetton for all his guidanceand mentorship. It has been heaps of fun working with Brian on very coolproblems. I’ll miss our casual research meetings and amusing non-academicconversations.Also, many thanks to my advisory committee of Sandra Barsky, LeahKeshet, Rob Kiefl, Michael Ward, and Dominik Schoetzau (my best wishesto you in your health). The helpful feedback, mathematical insights, andprofessional support you have offered throughout my graduate studies ismost appreciated!Thank you to the university examiners, Philip Loewen and Carl Ollivier-Gooch, and the external examiner, Thomas Witelski, for your feedback onthis manuscript.I’d also like to thank George Bluman, Randy LeVeque, and Aaron Froesefor helpful discussions.Thank you to Daniel Coombs and David Wilkinson for collaborations oninteresting and rewarding side projects and publications.xixDedicationThere are many people I’d like to thank on a personal level for their supportand helpful input throughout the toils of graduate school.Thank you to my family, and especially my Mom (rest in peace) andDad for your enthusiasm and encouragement.Many thanks to my fellow Institute for Applied Math graduates: IainMoyles, Kai Rothauge, Bernhard Konrad, Alastair Jamieson-Lane, ColeZmurchok, Jia Gou, Jielin Zhu, Simon Tse, William Thompson, and oth-ers; and thanks to my fellow graduate students both inside and outside ofmathematics who I’ve gotten to know and socialize with over the past years:Carmen Bruni, Kyle Hambrook, Sean Mehta, Cindy Blois, Thomas Wong,Matt Coles, and dozens more. It has been loads of fun to share research andteaching woes, to (jokingly?) discuss how bad grad school is, and to go forhikes, play games, and eat sushi together.Thanks to Alex Fang for your enthusiasm in assisting in my research.Our work on superconductivity will make a nice publication. It was a realpleasure to work with you as both a collaborator and a friend.To my students over the years: thanks for being so much fun. It hasbeen great to see you learn and I’ve quite enjoyed getting to know many ofyou. Being able to teach has really made grad school a more enjoyable time- seriously!xxDedicationThanks to my friends at Sprouts: Adrien, Mustafa, Harlan, Mikaela,Josh, Sean, Steve, Chelsea, Tim, Sarah, Lianne, Luc, and many others; andmy group of meditation friends: JP, Katherine, Al, Isabelle, and Hrvoje.A healthy mix of non-academic activities has benefited my sanity over theyears...And thanks to many others out there!xxiPart IIntroduction and Motivation1Chapter 1Context for This Work1.1 Scope of This ThesisThis thesis focuses upon one of the bigger projects I have been workingon during my PhD studies. It relates to modelling the design work of amagnetized target fusion reactor - a new and promising means for creatingconditions for and harvesting energy from nuclear fusion on earth.The work presented here comes from a variety of sources. The vast ma-jority or the entirety of papers I have authored appear in this thesis, withadditional information provided for a more detailed exposition of the math-ematical thought processes. Due to the fact that some journals have verytargeted audiences and/or page restrictions, there are important backgroundelements that did not/will not appear in the published form, and the otherchapters of this thesis aim at providing this extra background and otherinsights.Section 1.2 provides an introduction and motivation; chapters 2 and 3provide a mathematical and physical background to the analysis undertaken;chapters 4 and 5 are taken from publications with further details and exam-ples included in chapter 6. Finally, chapter 7 provides a summary and anoutlook at possible future work.21.2. Magnetized Target Fusion1.2 Magnetized Target FusionMuch of this section is taken from [26]. The purpose is to introduce the con-cept of magnetized target fusion and some important design considerations.1.2.1 Fusion MotivationFinding sustainable, clean energy sources is of the utmost importance forsociety. The use of fossil fuels has harmful effects upon the environmentin releasing greenhouse gases [29] and the damages done to the earth andecosystems in their harvesting; they are becoming more and more costly; andtheir supply is limited [31]. The investment in alternative energy sources isone of the most exciting areas of research in recent decades [11] includingthe rapid development of fuel cells [4], wind-driven power stations [2], solarpower [24] and geothermal power [16]. These technologies are improvingeach day, but there is another source of power that as of yet has not beenrealized on earth: nuclear fusion [18].Nuclear fusion is a complicated process in which lighter particles arefused together, often at very high temperatures and pressures, to produceheavier particles. Through a change in rest mass between the reacting andresultant particles, energy is released through radiation. The major chal-lenge in producing fusion is forcing the lighter particles to interact enoughwith each other so that they can fuse. Modern ideas of fusion are often basedon using deuterium (2 H) and tritium (3 H) plasmas as fuel and fusing themto produce helium and free neutrons [3] according to21H +31H→42He +10n+ 17.6MeV.31.2. Magnetized Target FusionOther fuel sources have also been considered, such as using only deuterium,or a special process known as Muon Catalyzed Fusion µCF. With µCF, anegatively charged muon facilitates the bonding of a proton and deuteriumto form 32 He and the release of the muon [15]. Unfortunately, the muons canalso be trapped in the bonds they form and thereby cease in their catalysisof fusion.Nuclear Fusion is how the sun generates its energy, with the conversionof hydrogen into helium [32]. Much of the work being done upon nuclearfusion for energy on earth is at the theoretical level, or attempting a proof-of-concept. The engineering of devices theoretically capable of fusing atomicnuclei is complex, time-intensive, and costly. As a result, it is important toproceed with as thorough an understanding of potential devices as possible,often through modelling work. The motivation for this work is in the math-ematical modelling and analysis of the design of a magnetized target fusionreactor being designed at General Fusion in Burnaby, BC, Canada.1.2.2 Magnetized Target FusionThere are a number of ongoing large-scale efforts to achieve nuclear fusion,including work at the National Ignition Facility (NIF)[7] and the Interna-tional Thermonuclear Experimental Reactor (ITER) project [9]. This paperreports on an effort by a Canadian nuclear fusion research company GeneralFusion [19]. The main idea is to compress a deuterium-tritium (DT) plasmalocated in a cylindrical cavity by an imploding shell of lead-lithium. Withsufficient compression, the plasma will heat to ignition, and a recovery sys-tem will convert the heat and radiation generated to electricity. There area number of factors to consider in the design of such a system. The work41.2. Magnetized Target Fusiondone on this topic in this thesis focuses on developing a basic model of thecompression so that we can apply finite volume numerical schemes for non-linear hyperbolic conservation laws with moving boundaries to make basicpredictions of the reactor being designed and thereby do a parameter searchto optimize the conditions under which fusion can occur. We will also applyasymptotic analysis to gain analytic insight into the factors influencing thedesign performance.One engineering obstacle is that heating the fuel to a high temperaturegreatly increases convective cooling, and any contact a plasma makes withthe walls of its confinement yields significant energy losses [3]. Fortunately,it is possible to suspend a plasma with magnetic fields and minimize itscontact with its enclosure. This insulation helps to reduce energy losses.Beyond the mere concern of energy losses is the necessary particle den-sity and pressure necessary for fusion to occur. One useful criterion is theLawson condition for fusion (whereby more energy is released than goes intoheating), which states that during a compression the product of number den-sity and time should exceed 1014 s/cm3 for deuterium-tritium fusion and thetemperature should exceed 6 − 10 keV [21]. This is sometimes generalizedinto a triple product where temperature times density times confinementtime must exceed ∼ 4× 1015 s keV / cm3 [28]. The product is even higherfor deuterium-deuterium fusion. In the numerical work that we present, wecompare our predictions of plasma densities and pressures with this tripleproduct by integrating the product of number density and temperature overthe time the plasma is being compressed.The design we consider consists of a 3-m diameter spherical region ofmolten lead-lithium, rotating about a 0.4-m diameter central vertical axisat high velocity such that a nearly empty, vertical cylindrical region is formed51.2. Magnetized Target Fusiondown the central axis. Within this axis, at a very precise time, a toroidalplasma of deuterium and tritium is released and reaches the very centreof the spherical region (we assume the plasma can survive for hundreds ofmicroseconds). Also at a precise time, pistons, delivering impact pressuresof 2 GPa hit the outside of the sphere, and a pressure wave travels throughthe lead-lithium until reaching the central axis at the precise moment theplasma arrives. The intense pressure and radially inward particle velocitycompress the plasma to cause fusion and release energy. The use of lead-lithium is appropriate as lead has a high cross section for neutron captureand absorbs radiation.The initial total pressure (gas plus magnetic pressure) of the plasma is 5MPa with particle density 3.2×1016 particles per cm3 and temperature of 100eV. The intense pressure pulse is modelled by a Gaussian with characteristictime 45 µs. A diagram illustrating the reactor is given in figure 1.1.Given the design as described, we wish to gain insights into how it willoperate in practice and what design specifications will make it more effective.This is taken up in the subsequent chapters.61.2. Magnetized Target FusionFigure 1.1: Magnetized Target Fusion Reactor Design7Part IIBackground to the Model forMagnetized Target Fusionand Finite Volume Methodsfor Hyperbolic ConservationLaws8Chapter 2Modelling2.1 Overall ComponentsSections 2.3 and 2.4.3 are derived from [26].The system introduced in section 1.2 has a variety of components thatneed to be treated together. The pistons impact the outer wall of the lead-lithium sphere and a localized disturbance in the fluid propagates; for thisunderstanding we need basic fluid dynamics equations. In our case, we willdeal with the compressible Euler equations. With the Euler equations, it isnecessary to know the relationship between the fluid pressure and its density,and this requires a model for the so called equation of state for lead-lithium.The plasma is another component to model, where we need to calculateits pressure and basic properties. This chapter will derive the equationsrelevant for our study.2.2 Deriving the Euler EquationsFor these derivations, we will consider a coordinate system (x, t) describingposition x ∈ R3 at time t ∈ [0,∞). We let the scalar ρ(x, t) ∈ [0,∞) bethe mass density (mass per unit volume), the vector field v(x, t) ∈ R3 bethe local fluid velocity, and P (x, t) ∈ R be the pressure. We denote Ω asa control volume in space with ∂Ω denoting is boundary with unit normal92.2. Deriving the Euler Equationsvector n pointing outwards. Note that n depends on the position upon thesurface.2.2.1 Mass ConservationIn this setting, the total amount of mass within Ω isM =∫Ωρ(x, t)dx.As mass cannot be spontaneously created or destroyed in our setting, theonly possible cause for M to vary is for mass to enter or leave Ω.At a point on ∂Ω, the amount of mass flowing into Ω per unit area perunit time, the mass flux, is −ρv · n. To understand this, we first considerthat only velocities that are parallel to n will carry fluid into or out of Ωsince any fluid flowing parallel to the surface of ∂Ω will not cross it. Thevelocity in the unit normal direction is v ·n. A positive value of v ·n indicatesthat fluid is leaving because it is flowing in the outward direction. In effect,v · nda, for a differential area element da, describes the volume of fluid thatis expelled, normal to the surface, per unit time. We can multiply this byρ to obtain ρ(v · n)da to find mass that leaves ∂Ω through the area elementda per unit time. Similarly, the mass flux that enters Ω is −ρv · n.By integrating this mass flux over the boundary of Ω, we can find thenet influx of mass. This can be equated to the time derivative of the totalmass. We have:ddt∫Ωρdx =∫∂Ω−ρ(v · n)daand by moving the time derivative inside the first integral by making it a102.2. Deriving the Euler Equationspartial derivative and using the divergence theorem on the second integral,we have∫Ωρtdx = −∫Ω∇ · (ρv)dx=⇒∫Ω(ρt +∇ · (ρv))dx = 0.For this to hold everywhere in space, for all control volumes Ω, it must holdpointwise so thatρt +∇ · (ρv) = Momentum ConservationThis derivation is similar, but there are subtleties, most notably in that weneed to consider a Lagrangian reference frame. At one point in this deriva-tion we will take a time derivative of an integral whose domain and integrandare time-dependent and we will use the Reynolds Transport Theorem for this[22].We consider a control volume Ω(t) that is time-dependent. We defineΩ(t) as the set of all positions of fluid particles that originated in the setΩ(0) at t = 0. We will consider a single component of the momentum inthe e-direction where e is a unit vector, and denote φ = ρv · e to be thee-momentum density per unit volume. We denoteΦ(t) =∫Ω(t)φdxto be the total e−momentum. As pressure is a force per unit area normal toan interface, we can compute the total force in the e−direction acting upon112.2. Deriving the Euler Equationsthe ensemble of particles within Ω(t) by a surface integralF (t) =∫∂Ω(t)−P (e · n)da(where the negative sign accounts for the fact that the outward normalpoints in the opposite direction to the force imparted by the boundary ofthe control volume upon the contained fluid). Newton’s second law says thatthe net force acting on a body is the rate of change of its total momentumand thusddt∫Ω(t)φdx︸ ︷︷ ︸Use Reynolds Transport=∫∂Ω(t)−P (e · n)da∫Ω(t)∂tφdx+∫∂Ω(t)φv · nda = −∫Ω(t)∇ · (Pe)dx∫Ω(t)∂tφdx+∫Ω(t)∇ · (φv)dx = −∫Ω(t)∇P · edxwhich finally implies, since the relation holds for all control volumes Ω(t)thatφt +∇ · (φv) +∇P · e = 0and upon using the fact that φ = ρv · e we havee · (ρv)t +∇ · (ρ(v · e)v) + e · ∇P = 0.Since this holds for all directions e, in particular e = xˆ, yˆ, zˆ, then(ρv1)t(ρv2)t(ρv3)t+∇ · (ρv1v)∇ · (ρv2v)∇ · (ρv3v)+PxPyPz=000.122.3. Equation of State for Lead-LithiumTable 2.1: Data for quadratic equation of state. Using these data points, weare able to construct an empirical fit to the equation of state for lead.Density Value11340 kg m−3 Pressure = 101325 N m−211340 kg m−3 Sound Speed = 2090 m s−116017 kg m−3 Pressure = 4.06185× 1010 N m−2Here we are using the notation that vj is the jth component of v.Often this is written in a more compact form as(ρv)t +∇ · (ρv ⊗ v) +∇P = 0where we interpret the tensor product ⊗ by a⊗ b = abT, a matrix, if a andb are column vectors.2.3 Equation of State for Lead-LithiumThe Euler equations for mass and momentum conservation are only a closedset of equations if there is a known function that relates the pressure P tothe density ρ.Based on the experimental data for lead [34], we fit a quadratic equationof state for lead-lithium using the density ρ0 and sound speed√∂P∂ρ [39] at1 atm, and the density at a pressure of 4.06185× 1010 Pa. This gives us twoprescribed points in the density-pressure plane and a slope at a point, whichdefines a quadratic equation of state P = A(ρ− ρ0)2 +B(ρ− ρ0) +C. Thedata we use in this quadratic fit are given in table 2.1. The fit is plotted infigure 2.1.It is important to note that such an equation of state does yield nega-tive absolute pressures for density values much less than those at standard132.4. Equation of State for PlasmaFigure 2.1: Quadratic fit to experimental equation of state data.conditions. Interpreting a negative pressure is rather delicate, but someexperiments have suggested that for some liquids, a negative pressure canbe sustained for fast time scales [13]. Under these negative pressures, thefluid will cavitate (i.e., it will turn into a gas or there will be empty space).Due to this physical observation, we allow for the negative pressures. Weremark further that as can be seen in table 4.3, the fit parameters (suchas the sound speed used) do influence the model predictions. Further ex-perimental investigation into the equation of state of lead-lithium would beuseful in improving the modelling.2.4 Equation of State for PlasmaBecause the plasma consists of charged particles moving in a magnetic field,the pressure of the plasma has two components: a gas pressure and a mag-142.4. Equation of State for Plasmanetic pressure due to the conservation of magnetic flux. We consider thesecomponents separately now.2.4.1 Plasma Gas PressureThis derivation makes use of basic thermodynamics, a good reference forwhich is [35].Because the process is very rapid, and with a first-order approximationthat energy losses can be neglected, we will assume the plasma compressionis adiabatic i.e. no heat is transferred into or out of the plasma during itscompression. Under this assumption, the gas pressure Pg and volume V ofthe plasma will obeyPgVγ = constantwhere γ = 5/3.This is a standard result assuming an ideal gas law and adiabatic con-ditions, but for a thorough coverage, we derive it below. Specific Heat RelationsThe first law of thermodynamics states thatdQ = dU + PdV. (2.1)i.e., that the differential heat exchange between a system and its environmentdQ is the sum of the small change in its internal energy dU and the pressure-volume work done in expanding PdV where P is the pressure of the systemand dV is the differential change in volume [35].Thermodynamically, we can choose two variables as independent, allow-152.4. Equation of State for Plasmaing the rest of the system to be implicitly determined by these two. Bytaking P and the temperature T as independent, then we can writedQ = (∂U∂P)TdP + (∂U∂T)PdT︸ ︷︷ ︸dU+P ((∂V∂P)TdP + (∂V∂T)PdT︸ ︷︷ ︸dV) (2.2)where the subscripts indicate which variable was held constant. Likewise,by taking V and T as independent, we can writedQ = (∂U∂V)TdV + (∂U∂T)V dT︸ ︷︷ ︸dU+PdV. (2.3)The specific heats cP and cV are defined by (∂Q∂T )P and (∂Q∂T )V , i.e., theyare respectively the instantaneous rate of heat exchange between a systemand its surroundings with respect to temperature under constant pressureand constant volume conditions. These can be readily expressed by usingdP = 0 in equation (2.2) and dV = 0 in equation (2.3) and dividing theequations by the differential dT giving uscP = (∂U∂T)P + P (∂V∂T)P (2.4)andcV = (∂U∂T)V . (2.5) Ideal Gas Specific HeatsAn ideal gas has an internal energyU =f2nRT (2.6)162.4. Equation of State for Plasmawhere f is the number of degrees of freedom, n is the number of moles ofgas, and R is the ideal gas constant. This follow from the equipartition ofenergy whereby each degree of freedom for a particle in a system has anaverage energy of 12kT where k is the Boltzmann constant [17]. In this case,U = U(T ). This immediately implies that dU = f2nRdT and that regardlessof what is held constant ∂U∂T =f2R and so that with (2.5) we obtaincV =f2nR. (2.7)Having found a simple expression for cV , we also seek one for cP . By theideal gas law,PV = nRT (2.8)where n is the number of moles of gas. This means that differentially PdV +V dP = nRdT so that in holding P constant and dividing by dT ,P (∂V∂T)P = nRor, using this in conjunction with (2.4) and (2.5),cP = nR+ cV . (2.9)We have established that cP = nR+ cV for an ideal gas. Adiabatic ConditionsAdiabatic conditions occur when there is no heat exchange between a systemand its environment i.e. dQ = 0. Under these conditions, we will establisha relationship that is held between the pressure and the volume.172.4. Equation of State for PlasmaWe first note that from equations (2.6) to (2.8), U = cV T =cVnRPV sothatdU =cVnRPdV +cVnRV dP.If this is used in (2.1), for adiabatic conditions thencVnRPdV +cVnRV dP︸ ︷︷ ︸dU+PdV = 0and therefore(cV + nR)PdV = −cV V dPor that upon using cP = cV + nR and rearranging,cVdPP= −cPdVV.Integration yields P cV = constant× V −cP or thatP ∝ V −cP /cV . (2.10)In the case of a plasma, comprised of monatomic species, there are 3degrees of freedom, one for each independent spatial direction, and hencewe model the gas pressure byPg ∝ V−5/3since cP = 52nRT and cV =32nRT.182.4. Equation of State for Plasma2.4.2 Plasma Magnetic PressureOur model of a plasma is quite rudimentary and as a first approximationwe consider a solenoidal plasma. The essence of this model is that a plasmaconsists of moving charged particles, moving in circular orbits. The designof General Fusion is that of a spheromak - a toroidal plasma. What weaim to establish here is a qualitative model for what is termed the magneticpressure of the plasma. This is the pressure associated with the magneticfield induced by the motion of the charged particles.An adiabatic plasma preserves its magnetic flux [1]. This key fact is es-sential in deriving the pressure. We begin by considering a sheet of plasmawithin a cylindrical geometry where the particles travel along a path ofradius R. In cylindrical coordinates (r, θ, z), the current density is 0 every-where in space except at r = R where it takes the value λ Coulombs/unitlength of z in the unit vector direction θˆ. Throughout this derivation, thefigure 2.2 will be relevant. Magnetic Field CalculationFrom Maxwell’s equations [12], we have that ∇× B = µ0J where B is themagnetic field, J is the current density, and µ0 is the permeability of freespace. Then, by Stokes’ Theorem∮CB · dr =∫∫S(∇×B) · nda =∫∫Sµ0J · nda = µ0I (2.11)where C is a closed loop with differential path displacements dr, S is aspanning surface to C with normal vector n given by the right-hand rule, dais the differential area, and I is the total amount of current passing through192.4. Equation of State for PlasmaFigure 2.2: The geometry of our simple sheet of plasma. The current densityis λ. The closed loops we consider are labelled, as is the differential surfaceelement considered for computing the pressure.S in the normal direction.We consider 3 rectangular loops Ci, i = 1, 2, 3, with vertices V i1 (r =r−i , θ = 0, z = −`/2), Vi2 (r = r+i , θ = 0, z = −`/2), Vi3 (r = r+i , θ = 0, z =`/2), V i4 (r = r−i , θ = 0, z = `/2) where for loop C1, r−1 = r−, r+1 = r(0);for loop C2, r−2 = r−, r+2 = r+; and for loop 3, r−3 = r(0), r+3 = r+ withr− < −R < r(0) < R < r+. Note that by the axial symmetry, there is noθ−dependence in any of these computations so we chose θ = 0, and dueto the infinite length there is no z−dependence, but for greater physicalinterpretation we have chosen the z−range to go from −`/2 to `/2 with` > 0 arbitrary. As each individual current component is in the θˆ-direction,this gives rise to magnetic fields in the z− and r−directions but not theθ−direction. Furthermore, by considering any point in space and looking atthe contribution of the magnetic field r−component due to currents above202.4. Equation of State for PlasmaFigure 2.3: The r−direction runs horizontally and the currents run into thepage as indicate by the ⊗. The magnetic fields generated by the currentsare drawn in blue and within the red boxes, the r−components of the fieldcancel.and below it, there is perfect cancellation, and there should be no magneticfield in the r−direction. See figure 2.3.We will write B = Brrˆ +Bθθˆ +Bz zˆ. Then, in exploiting the symmetryof the system, we have that B = Bz zˆ and∮C1B · dr = −Bz(r = r−)`+Bz(r = r(0))`+∫ r(0)r−0dr +∫ r−r(0)0dr= `(Bz(r = r(0))−Bz(r = r−)) = µ0`λ. (2.12)212.4. Equation of State for PlasmaAlso,∮C2B · dr = −Bz(r = r−)`+Bz(r = r+)`+∫ r+r−0dr +∫ r−r+0dr= `(Bz(r = r+)−Bz(r = r−)) = 0 (2.13)since the net current is zero for this case (same amount going into and outof the surface), and∮C3B · dr = −Bz(r = r(0))`+Bz(r = r+)`+∫ r+r(0)0dr +∫ r(0)r+0dr= `(Bz(r = r+)−Bz(r = r(0))) = −µ0`λ (2.14)where the negative sign in the current arises from the current pointing inthe opposite direction to the outward normal vector for the loop integral.Solving equations (2.12) through (2.14) tells us that Bz(r = r−) = Bz(r =r+) = 0 and Bz(r = r(0)) = µ0λ. Besides their range, r−, r(0), and r+ werearbitrary so the z−component of the magnetic field must be µ0I everywhereinside the solenoid and zero everywhere outside. Preservation of Magnetic FluxThe magnetic flux is defined byΦB =∫∫SB · nda222.4. Equation of State for Plasmaso that if the magnetic field is uniform everywhere inside the solenoid andwe choose a surface with normal vector zˆ, thenΦB = piR2µ0λ.If the magnetic flux is to be conserved, we must have thatλ =ΦBpiR2µ0. (2.15) Calculating the Magnetic PressureWe consider a small segment of the solenoidal pathway of z−length ` sub-tending an angle dθ so that the arc length subtended is Rdθ. See figure 2.2.If the magnetic field acting upon this segment is Bz zˆ then for a currentIθˆ, the force acting upon the infinitesimal piece of the surface is, via theBiot-Savart law,dF = IRdθθˆ ×Bz zˆ = IBzRdθrˆ = λ`BzRdθrˆ. (2.16)If the magnetic field just inside the solenoid at r = R− due to theinfinitesimal piece is B(i) (so that at r = R+ it would be −B(i)) and themagnetic field due to the rest of the solenoid is B(s) then just outside thesolenoid at r = R+, we have −B(i) + B(s) = 0, and just inside r = R−, wehave that B(i) +B(s) = µoλzˆ. We find that the magnetic field acting on theinfinitesimal slice is B(s) = 12µ0λzˆ. Using this result in (2.16), we find thatdF =λ2µ0`Rdθ2rˆ.As pressure is normal force dF · rˆ divided by area `Rdθ, we find there is a232.4. Equation of State for PlasmapressurePm =λ2µ02=Φ2B2pi2R4. (2.17)Equation (2.17) tells us that in our simple model, the magnetic pres-sure of a plasma should scale like 1/R4. This holds in more generality foradiabatic compressions of plasmas with torus-like structures [40].2.4.3 Total Plasma PressureWe take the total pressure PL as the sum of the two component pressuresPL = Pg+Pm. A more detailed description of the plasma would be an impor-tant improvement for using this modelling in making accurate predictionson the device’s operation.24Chapter 3Brief Background into FiniteVolume Methods forConservation LawsThe purpose of this chapter is to provided a more detailed mathematical de-scription and background to the numerical method implemented in chapter4. There are many books and papers written explaining the implementationof numerical methods for hyperbolic conservation laws, but one particularlyuseful book is one written by Randy LeVeque[23]; this book gives a verythorough treatment of the topic and any interested reader should refer tosaid book.3.1 Finite Volume Schemes for ConservationLawsA finite volume method is one that describes cell averages in a meshing; afinite difference scheme describes the pointwise values, and a finite elementmethod transforms a differential equation into its weak form and uses a finiteset of simple functions. The theory behind finite volumes in application toconservation laws is vast, with deep results on the stability of methods, the so253.1. Finite Volume Schemes for Conservation Lawscalled total variation diminishing property (preventing spurious oscillations),the fact that linear methods that preserve monotonicity are limited to first-order accuracy, the importance of writing schemes in conservative form, andmuch more. While all of these items are important, this chapter will coverthe basic theory necessary for developing and understanding our numericalmethod presented in chapter ConvergenceBecause hyperbolic problems have solutions that naturally exist in a weaksense, due to the possibility of shock solutions, convergence only occurs inthe L1-norm. This means that if unum(x, t) denotes the numerical solutionand uex(x, t) denotes the exact solution, both at a fixed end-time t, then asthe mesh is refined, the errorE =∫|unum(x, t)− uex(x, t)|dx→ 0where the integral ranges in x over the set upon which we are solving. Apth-order method is one where E = O(hp) where the spatial mesh size isO(h).3.1.2 Finite Volume SetupWe will consider a simple conservation law, in one spatial dimension x withtime t, for a quantity u ∈ RD:ut + (f(u))x = 0. (3.1)263.1. Finite Volume Schemes for Conservation LawsFigure 3.1: Diagram of basic finite volume setup.This describes the time-evolution of a quantity with flux f(u(x, t)) ∈ RD.Let us suppose further that we wish to solve this hyperbolic partial differ-ential equation on the domain (x, t) ∈ [0, 1]× [0,∞) so that D = 1.We mesh the interval [0, 1] into points xi = ih, i = 0, 1, ..., N , whereh = 1/N is the mesh spacing and we allow half-integer indices such that,for example, x3/2 = 3h/2, etc. Let us define u¯(xi, tj) =1h∫ xi+1/2xi−1/2u(x, tj)dxto be the average value of the quantity u on an interval of width h centredat xi at time tj . See figure Finite Volume Schemes for Conservation LawsWe can describe the time evolution of u¯(x, t) bydu¯(xi, t)dt=1h(ddt∫ xi+1/2xi−1/2u(x, t)dx)=1h(∫ xi+1/2xi−1/2ut(x, t)dx)= −1h(∫ xi+1/2xi−1/2(f(u))xdx)where the last equality came utilizing the original equation (3.1). We canthus write thatdu¯(xi, t)dt= −f(u(xi+1/2, t))− f(u(xi−1/2, t))h.This simple evolution is the basic idea behind many finite volume meth-ods for conservation laws. Ultimately, we seek to find ways of stepping for-ward in discrete time steps from a state u¯ji , our approximation to u¯(xi, tj),to u¯j+1i in ways that are consistent with the physical system with high ac-curacy. This amounts todu¯(xi, tj)dt= −F ji+1/2 −Fji−1/2h(3.2)such that F ji−1/2 = F(u¯ji−1, u¯ji ) is a discretized flux that depends upon theapproximated average cell values to its left and right. For simpler notation,we will drop the overbars so that uji represents u¯ji .In the following sections, we consider linear and nonlinear systems, andhow low and high resolution schemes can be combined to yield a flux-limitedhigh-resolution finite volume method.283.2. Linear Systems3.2 Linear SystemsWe consider the evolution of the vector u(x, t) ∈ RD according tout + (Au)x = 0 (3.3)for a diagonalizable, real, constant matrix A ∈ RD×D with distinct realeigenvalues. Let us suppose that A can be diagonalized according toA = RΛL (3.4)where R is a matrix whose columns are the right eigenvectors, L = R−1 is thematrix whose rows are the left eigenvectors ofA, and Λ = diag(λ1, λ2, ..., λD).For the sake of providing a simple example, we will assume that D = 2and A has two eigenvalues λ± with λ+ = λ1 > 0 and λ− = λ2 < 0 withright eigenvectors r± and left eigenvectors `±.3.2.1 Low ResolutionHere we consider a scheme that provides first-order accuracy for linear sys-tems. Using equation (3.4) and multiplying equation (3.3) by L, we have:vt + Λvx = 0 (3.5)where v = Lu gives a representation of u with respect to the basis col(R).This gives us two rather trivial equationsv±t + (λ±v±)x = 0293.2. Linear Systemswhere v± = `±u are scalars (recall `± are row vectors and u is a columnvector).From the method of characteristics, observe that the system informationpropagates along the lines dxdt = λ±. This is critical to our next step. Theequation (3.2) tells us that if we have discretized u and v according tothe finite volume approach with uji , vji = Luji , then a suitable first-orderapproximation to estimating vj+1i with a forward time step size k is:v±j+1i = v±ji −kh(F±ji+1/2 −F±ji−1/2) (3.6)At each point, the true flux is λ±v±, however from the characteristicinformation it makes sense to defineF+ji−1/2 = λ+v+ji−1F−ji−1/2 = λ−v−ji (3.7)so that we use the flux as defined on the left side of the interface wheninformation is right-going and the flux as defined on the right side of theinterface when the information is left-going. The justification of this ideacomes when considering each cell as a piecewise constant and our solving theRiemann problem between two adjacent cells. This is solving the Riemannproblem at each interface is referred to as Godunov’s method.Finally we can combine equations (3.6) and (3.7) along with uj+1i =303.2. Linear SystemsRvj+1i to writeuj+1i = uji −khR(λ+(v+ji − v+ji−1), λ−(v−ji+1 − v−ji ))T= uji −kh(λ+(v+ji − v+ji−1)r+ + λ−(v−ji+1 − v−ji )r−)= uji −kh(λ+w+ji−1/2 + λ−w−ji+1/2) (3.8)where w±ji−1/2 is the projection of uji−uji−1 onto r±. Often the w’s are calledfluctuations.Keeping with this idea of projection, there is a nicer way to write equa-tion (3.8) in definingA+ = Rdiag(λ+, 0)L (3.9)andA− = Rdiag(0, λ−)L. (3.10)We writeuj+1i = uji −kh(FLji+1/2 −FLji−1/2) (3.11)whereFLji−1/2 = A+uji−1 +A−uji (3.12)is the low resolution flux. In general we can writeFLji+1/2 −FLji−1/2 = A+(ui − ui+∗) +A−(ui−∗ − ui) (3.13)where i+∗ = i−1 is an upwind based on a positive eigenvalue and i−∗ = i+1is the upwind index based on a negative eigenvalue.313.2. Linear SystemsNote thatA = A+ +A−. (3.14)Later on, we will make use of the notation|A| = A+ −A−. (3.15)Because A+A− = Rdiag(λ+, 0)diag(0, λ−)L = 0 and likewise for A−A+ wehaveA2 = |A|2. (3.16)With D > 2, we can define A+ = Rdiag(λ1sgn(λ1), ..., λDsgn(λD))L andA− = Rdiag(λ1sgn(−λ1), ..., λDsgn(−λD))L.The following remark is more relevant to nonlinear systems but we men-tion it here: the low resolution flux (3.12) can be expressed byFLji−1/2 = A+uji−1 +A−uji= (A−A−)uji−1 +A−uji= Auji−1 +A−(uji − uji−1)= Auˆji−1/2where uˆji−1/2 = ui−1 + RI−L(uji − uji−1) is an upwinded u-value, and I− =diag(0, 1). Clearly,uˆji−1/2 = uji−1 + w−ji−1/2. (3.17)What this rewriting tells us is that we can obtain an upwind (low) resolutionflux by evaluating our flux function Au at an upwinded value of u.323.2. Linear Systems3.2.2 CFL ConditionThe low resolution method above is stable provided the CFL condition isobeyed, namely that the distance information could propagate at maximumover the time step of size k, max{|λ+|, |λ−|}k, is smaller than the spatialmesh size. For nonlinear systems for any D, the same condition must beupheld but we seek the maximum |λi| over all eigenvalues λi and all meshlocations.3.2.3 High ResolutionHaving established a first-order method given in equation (3.8), we willobtain a method that has a higher order accuracy by making use of Taylorseries. In what follows, we assume that in our region of interest, the solutionis smooth. Given the nature of hyperbolic problems whereby shocks canform, this is not always true.Beginning with (3.3), under the assumption of smoothness of u and thefact that A is a constant matrix, we can writeut = −Aux (3.18)andutt = −Auxt= −A(ut)x= −A(−Aux)x= A2uxx. (3.19)333.2. Linear SystemsBy Taylor series, we also have u(x, t+ k) = u(x, t) + kut(x, t) + k22 utt(x, t) +O(k3) so thatu¯(xi, tj+1) = u¯(xi, tj) +kh∫ xi+1/2xi−1/2(−Aux)dx+k221h∫ xi+1/2xi−1/2A2uxxdx= u¯ji −kh(Au(xi+1/2, tj)−Au(xi−1/2, tj))+k221h(Aux(xi+1/2, tj)−Aux(xi−1/2, tj)) (3.20)where in computing the cell averages we replaced the partial time derivativeswith partial space derivatives. To second-order, we have that u(xi−1/2, tj) =12(uji−1 +uji ) and ux(xi−1/2, tj) =1h(uji −uji−1) and thus with (3.20) we haveuj+1i = uji −kh(FHji+1/2 −FHji−1/2) (3.21)where we have that in general the high resolution flux can be expressed asFHji−1/2 =A2(uji + uji−1)−k2hA2(uji − uji−1). (3.22)This form is also known as the Lax-Wendroff method.This is second-order accurate, provided the solution is sufficiently smooth.In the next section, we will consider how to combine the low and high reso-lution fluxes, taking advantage of the good properties of both.3.2.4 Flux LimitersThe idea behind flux limiters is to use the high resolution flux when thenumerical solution is smooth and to use the low resolution flux, which isnot based in Taylor series, when the solution is not smooth or when it isdiscontinuous. First, we note that there is a simple way of expressing the343.2. Linear Systemsdifference between the high and low resolution fluxes (see equations (3.12),(3.9), (3.10), (3.15), (3.16), (3.14), and (3.22)):FHji−1/2 −FLji−1/2 =A2(uji + uji−1)−k2hA2(uji − uji−1)−A+uji−1 −A−uji=A+ +A−2(uji + uji−1)−k2hA2(uji − uji−1)−A+uji−1 +A−uji=A+2uji −A−2uji −A+2uji−1 +A−2uji−1−k2h|A||A|(uji − uji−1)=12|A|(I −kh|A|)(uji − uji−1).Let us define∆ji−1/2 =12|A|(I −kh|A|)(uji − uji−1). (3.23)Then, depending on how smooth the solution is, the ideal numerical fluxcan be written asF ji−1/2 = FLji−1/2 +12|A|(I −kh|A|)∆˜ji−1/2 (3.24)with∆˜ji−1/2 = RΘji−1/2L(uji − uji−1) (3.25)and the matrixΘji−1/2 = diag(θ+ji−1/2, θ−ji−1/2) (3.26)limits the fluxes based on how rapidly the solution is changing in its differenteigencomponents. If the solution is perfectly smooth in the r+ space then353.2. Linear Systemswe should take θ+ = 1; and if, for example, it is highly discontinuous thenθ+ = 0. We can think of RΘji−1/2L(uji − uji−1) as a limited net fluctuation.In our example, we consider the ratioρ±ji−1/2 =v±ji−1−v±ji−2v±ji−v±ji−1, if λ± > 0v±ji+1−v±jiv±ji−v±ji−1, if λ± ≤ 0=r±·w±ji−3/2r±·w±ji−1/2, if λ± > 0r±·w±ji+1/2r±·w±ji−1/2, if λ± ≤ 0.(3.27)If this ratio is negative or near zero then it means the derivative at our celli is zero or has a discontinuity; if this ratio is near 1 then the derivativeis not changing very quickly. If the ratio is very large and positive then itmeans that the solution is rapidly changing. One particular choice of limiter,known as the minmod limiter, is given byθ±ji−1/2 = max{0,min{1, ρ±ji−1/2}}. (3.28)This limiter works in our setting, but many other limiters are possible suchas, for example, the Van Leer limiter withθ±ji−1/2 =ρ±ji−1/2 + |ρ±ji−1/2|1 + |ρ±ji−1/2|.Clearly the minmod limiter will give full weight to the low resolution fluxif there is a spike in the solution and near a local max/min; however, ifthe solution changes rapidly but maintains the sign of its derivative, it ispossible for the limiter to give a strong weighting to the high resolution flux.363.3. Extensions to Nonlinear Systems3.2.5 ExtrapolationBoundary conditions are an important part of solving partial differentialequations. If we assume that we have unknown cell values uji where i =0, 1, ..., N then the physical boundary occurs at x = 0 where uj0 is located.In order to update u0, we need the flux function Fj−1/2 and this requiresus to know ρ±j−1/2, which in turn requires a value for vj−2. An identicalissue occurs at the right boundary. In practice, a constant extrapolation issufficient to maintain at least first-order accuracy at the boundary. Thus,at each time step when required, if we don’t otherwise have a component ofuji specified, we can define uj−2 = uj−1 = uj0. We can use the same approachat the right boundary.3.3 Extensions to Nonlinear SystemsHaving established the basic numerical approach to linear systems, this sec-tion is devoted to extending the methods presented to nonlinear systemsand giving examples.3.3.1 General IdeaNonlinear systems can be dealt with in much the same way as we haveconsidered in the case of a linear system with a few small differences. Thesystem takes the formut + (f(u))x = 0.When the solution is smooth, and even in the case of shocks the solutiontends to be smooth everywhere except for isolated points, we can define alinearized matrix A = ∂f∂u so that ut +A(u)ux = 0. Because the eigenvectors373.3. Extensions to Nonlinear Systemsand eigenvalues are changing at each grid point in general (the matrix A isdependent upon u), we need to consider the cell updates locally. Again weconsider the case where A is a 2×2 matrix with eigenvalues λ± with λ+ > 0and λ− < 0.In averaging the eigenvalues of the linearized system at cells i and i− 1,we determine the direction of the system in the basis of eigenvectors ofthe matrix A(ui−1/2) where ui−1/2 =12(ui + ui−1). Note that we use theaverage of the eigenvalues and not the eigenvalues of the average of A(uji−1)and A(uji ). We can define an upwind value of uji−1/2 by (see equation (3.17))uˆji−1/2 = uji−1+A−ji−1/2(uji−uji−1) and a corresponding upwind flux functionfˆ ji−1/2 in evaluating f at uˆji−1/2. This establishes a first-order method. Notethat in evaluating this flux function at the upwinded value of u, we areimplementing an approximate Riemann solver: instead of solving the fullRiemann problem at each interface, we only solve an approximate problem.In practice this works quite well.In the local basis, we can also compute w±ji−1/2 and, based on the di-rections, w±ji∗ (the upwind fluctuations), allowing us to define a local fluxlimiter. Note that because the w±’s are vectors and they may change fromone cell to the next, to compare them, we use a projection and defineρ±ji−1/2 =w±ji∗ · w±ji−1/2|w±ji−1/2|2.Then, we can define local flux limiters and include the high resolution fluxdefined by (3.24), (3.25), and (3.26), replacing all matrices that were pre-viously held as constant by their corresponding values at i − 1/2. For fullexposure, we are computing w±i−1/2 and w±i±∗ where i∗ is either i − 3/2 ori + 1/2, depending on the sign of the eigenvalue λ±i−1/2. The i±∗-value is383.3. Extensions to Nonlinear Systemsx0 x1 ... xi−1 xi xi+1 ... xN−1 xNu0 u1 ... ui−1 ui ui+1 ... uN−1 uNE−2 E−1 E0 E1 ... Ei−1 Ei Ei+1 ... EN−1 EN EN+1 EN+2D− 32 D− 12 D 12 ... Di− 12 Di+ 12 ... DN− 12 DN+ 12 DN+ 32I− 32 I− 12 I 12 ... Ii− 12 Ii+ 12 ... IN− 12 IN+ 12 IN+ 32w− 32 w− 12 w 12 ... wi− 12 wi+ 12 ... wN− 12 wN+ 12 wN+ 32θ− 12 θ 12 ... θi− 12 θi+ 12 ... θN− 12 θN+ 12φ− 12 φ 12 ... φi− 12 φi+ 12 ... φN− 12 φN+ 12w˜− 12 w˜ 12 ... w˜i− 12 w˜i+ 12 ... w˜N− 12 w˜N+ 12uˆ− 12 uˆ 12 ... uˆi− 12 uˆi+ 12 ... uˆN− 12 uˆN+ 12fˆ− 12 fˆ 12 ... fˆi− 12 fˆi+ 12 ... fˆN− 12 fˆN+ 12Table 3.1: An illustration of the grid variables. The solution and true spatialmesh are defined at indices i = 0, 1, ..., N . By constant extrapolation, weextend the solution and the spatial mesh giving the E-variables with indicesi = −2,−1, 0, ..., N,N + 1, N + 2. By subtracting these E−values and av-eraging them, we obtain difference values D and interpolated values I bothwith indices −3/2,−1/2, ..., N + 1/2, N + 3/2. With the I− values, we canalso interpolate a linearized system and determine the components of the D-values in the coordinates of the right eigenvectors, giving us the w variables.By upwinding, we determine the fluctuation ratios and the limiters θ andφ with indices i = −1/2, 1/2, ..., N − 1/2, N + 1/2, and we can also definelimited w-variables, the w˜-variables. In making use of the matrices formedat the I−cells, we can also obtain the upwind E-values, uˆ, relevant to definethe fluxes fˆ . A second-order flux limited method would then determine thelow resolution numerical flux, fˆ along with a second-order flux correctionthat is suitably limited.always chosen to take data from the side, from which, the information iscoming. The second-order flux correction makes use of the limited w-values.To illustrate the locations of the various grid variables, we present table Example of Different Methods as Applied to Burger’sEquationTo illustrate the implemention of 3.3.1, we consider the inviscid Burger’sequationut + (12u2)x = 0393.3. Extensions to Nonlinear Systemson (x, t) ∈ [−3, 3]× [0, 3/2] with initial conditionsu(x, 0) = H(x+ 3/2)−H(x− 3/2)where H denotes the Heaviside step function. From the analytic solution,we can consider the x−range to be a subset of R where the true solution lies.We focus on the qualitative behaviour of the different numerical methods.Numerically our boundary conditions are implemented by constant extrap-olation as described above in 3.2.5. We will begin by solving the equationanalytically; then we will plot the solutions obtained in applying a first-orderupwinding method, a well-known method that is extremely low resolution(Lax-Friedrichs), the second-order flux limited scheme explained in 3.3.1,and finally, and second-order method that is not limited.From characteristics, in smooth regions, we have thatut + uux = 0 (3.29)so that along curves such thatdxdt= u,thendudt= 0.Thus, alongx = ξ + u(ξ, 0)t,u = u(ξ, 0)403.3. Extensions to Nonlinear Systemsfor ξ ∈ R. This breaks down in two areas. At t = 0, the characteristic curvescoming from (3/2, 0) overlap: at x = 3/2+, the characteristics are verticaland at x = 3/2−, the characteristics move to the right with unit speed.This makes u multivalued and there is a shock. Also, the curves comingfrom (−3/2, 0) fail to define values for u in the sector (x + 3/2) = ct for0 < c < 1. To deal with this, we insert an expansion fan.With a flux of f(u) = 12u2, the Rankine-Huginoit conditions [36] tell usthat the shock emitted from (3/2, 0) travels at speed v = [f(u)][u] =1202− 12120−1 =12 where [u] = u+−u− measures the discontinuity in the quantity u across theshock where u+ = limx→x+s u(x, t) and u− = limx→x−s u(x, t) are the limitingvalues of u to the right and left of the shock at x = xs(t) respectively.For the expansion fan, we fit a similarity solution u = u(η) where η =(x + 3/2)/t. Then ∂t = −(x + 3/2)/t2Dη and ∂x = (1/t)Dη. Using this in(3.29) gives−(x+ 3/2)/t2u′(η) + uu′(η)/t = 0oru′(η)(−η + u) = 0so that either u′(η) = 0 or u = η. A constant value for u does not work,but u = η = (x + 3/2)/t does. This fits the expansion fan. A diagramillustrating the qualitative behaviour of the solution is found in 3.2.413.3. Extensions to Nonlinear SystemsFigure 3.2: A sketch of the solution to the Burger’s equation illustrating theshocks and expansion fans.At t = 3/2, the analytic solution is thatu(x) =0 if x ≤ −3/223x+ 1 if − 3/2 < x < 01 if 0 ≤ x < 9/40 if x ≥ 9/4.Comparison of several computational approximations to this exact solu-tion are shown in figure 3.3. For our implementation, we used an approxi-mate Riemann solver when required to solve the Riemann problem and notthe exact Riemann solution. The other methods have been previously ex-plained, but we wish to state the Lax-Friedrichs algorithm [23]. This methodis popular due to its extreme simplicity in implementation, but it fails to423.3. Extensions to Nonlinear SystemsN LF G LTD LW25 0.96 0.38 0.22 0.4250 0.62 0.20 0.12 0.46100 0.38 0.12 0.06 0.37200 0.22 0.06 0.02 0.33Rate 0.7 0.9 1.1 0.1Table 3.2: The L1-numerical errors are computed using the trapezoidalrule where LF denotes the Lax-Friedrichs scheme, G denotes the Godunovscheme (first-order), LTD denotes the limited second-order method, and LWdenotes the unlimited second-order method (Lax-Wendroff). Based on theslope of the best fitting line of the log-log plot of the error vs N , we alsonote the convergence rate.provide an accurate solution without refining the spatial domain immensely.This method tends to smear out shocks and is, at best, first-order accurate.Given a mesh of cell averages uji , the Lax-Friedrichs update is:uj+1i =12(uji−1 + uji+1)−k2h(f(uji+1)− f(uji−1)).Although simple, it performs quite poorly as seen in figure 3.3. Also worthnoting in the figure is that the second-order flux-limited scheme matchesthe solution the closest, and that without applying a flux-limiter, there areoscillations in the computed solution.In table 3.2, we plot the convergence rates of the different methods. Assecond-order methods can only be second-order away from shocks and withsmooth solutions, we do not obtain second-order convergence even with thesecond-order method. However, it does outperform the other methods. Wesee that the Lax-Friedrichs method is the least accurate, and that Lax-Wendroff may even fail to converge.433.3. Extensions to Nonlinear SystemsFigure 3.3: An illustration of different numerical solvers applied to the in-viscid Burger’s equation. The exact solution is plotted in green; the first-order Godunov method, based on upwinding, is plotted in red; the first-order Lax-Friedrich’s scheme is plotted in dark blue; the second-order flux-limited method with approximate Riemann solver is plotted in light blue;and the second-order Lax-Wendroff method without limiters is plotted inblack. N = 50 mesh points were chosen.443.3. Extensions to Nonlinear Systems3.3.3 Example of the Importance of Conservative Form andProblems with Geometric Source TermsWe consider here a simple example to illustrate how a scheme being inconservative form is more true to the analytic solution than one that isnot in conservative form. This also touches on the topic of geometric sourceterms and how they tend to reduce accuracy. We consider a simple signallingproblem in a sphere where r ∈ (0, 1] and t ≥ 0:(r2u)t − (r2u)r = 0withu(r, 0) = 0, u(1, t) = h(t) = t2(1− t).To solve this, we consider w = r2u so thatwt − wr = 0, with w(r, 0) = 0, w(1, t) = h(t).If w = w(r(ξ, η), t(ξ, η) then along ∂r∂ξ = −1 and∂t∂ξ = 1, so thatr = −ξ + r0(η)andt = ξ + t0(η),we find∂w∂ξ= wt − wr = 0 =⇒ w = w(ξ = 0, η).As η is a free-parameter and we know w(1, t), we choose r0 and t0 to be 1and η, respectively so that at ξ = 0 we can parameterize all of the line r = 1453.3. Extensions to Nonlinear Systemswith t ≥ 0 by (ξ = 0, η). This means that r = 1− ξ and t = ξ + η and thusη = r+ t− 1 and w(0, η) = w(0, r+ t− 1) = h(r+ t− 1). We thus find thatu(r, t) =(r − 1 + t)2(1− (r − 1 + t))r2H(r − 1 + t)where the Heaviside function H arises because there is no forcing for t < 0so that h(t) = H(t)h(t).Numerically, we consider a variety of approaches, some conservativewhere we solvewt − wr = 0with w = r2u, and some where we solveut − ur −2ru = 0by using a split-step strategy [23] where at each time step we first stepforward for the equation ut − ur = 0 and then with the updated u-value,we step forward with ut = 2ru where u is the updated value. The term2ruis known as a geometric source term; when such terms cannot be avoidedthey often make a numerical method less accurate, formally reducing it fromsecond-order to first-order. The plots of the analytic solution and second-order flux-limited numerical solutions with and without a geometric sourceterm where N = 25 and t = 0.5 are plotted in figure 3.4. In our example,we just see a slightly less strong agreement with the numerics when thegeometric source term is not removed.As a point of interest, we also include table ?? showing the L1-convergence(as computed with the trapezoidal rule) of various numerical methods. Weindeed notice that with the geometric source terms, the second-order meth-463.3. Extensions to Nonlinear SystemsFigure 3.4: We observe that the numerical scheme in conservative form isslightly better than the numerical scheme with geometric source terms. Theplots depict the solution at t = 0.5 with N = 25 mesh points.ods are mostly first-order in nature, whereas in the conservative form, themethods are both second-order.two approaches. For the conservative approach, we solvewt − wr = 0where w = r2u using a flux-limited second order finite volume method, andfor the non-conservative approach we solveut − ur −2ru = 0by using a split-step strategy [23] where at each time step we first stepforward for the equation ut − ur = 0 and then with the updated u-value,we step forward with ut = 2ru where u is the updated value. The term2ruis known as a geometric source term; when such terms cannot be avoided473.3. Extensions to Nonlinear SystemsN G1 G2 GL2 C1 C2 CL225 1.1× 10−3 5.6× 10−4 4.4× 10−4 1.5× 10−3 3.6× 10−4 1.6× 10−450 6.0× 10−4 2.4× 10−4 1.9× 10−4 7.6× 10−4 1.1× 10−4 4.7× 10−5100 2.8× 10−4 9.0× 10−5 7.3× 10−5 3.5× 10−4 2.8× 10−5 1.1× 10−5200 1.4× 10−4 3.9× 10−5 3.4× 10−5 1.8× 10−4 7.6× 10−6 2.8× 10−6Rate 1.0 1.3 1.2 1.0 1.9 2.0Table 3.3: The L1-numerical errors are computed using the trapezoidal rulewhere G1 denotes a first-order scheme with a geometric source term; G2 de-notes a non-limited second-order scheme with a geometric source term; GL2denotes a flux-limited second-order scheme with a geometric source term;C1 denotes a conservative first-order scheme; C2 denotes a non-limited,conservative second-order method; and CL2 denotes a conservative, flux-limited second-order method. Based on the slope of the best fitting line ofthe log-log plot of the error vs N , we also note the convergence rate.they often make a numerical method less accurate, formally reducing it fromsecond-order to first-order. In our example, we just see a slightly less strongagreement with the numerics. See figure 3.4.48Part IIIJournal Articles49Chapter 4Numerical Finite VolumeInvestigationThis chapter is derived from [26]. A context to this work is found in chapter1.2. Details of the model equations used here are found in chapter 2. A morethorough explanation of numerical methods for hyperbolic conservation lawsis found in 3, and further details of the numerics implemented here are alsofound in section 6.1 and in the appendix A.4.1 Physical AssumptionsGeometrically, we adopt a spherically symmetric model in which the plasmais treated as a spherical gas, immersed in the lead-lithium whose outerboundary is a sphere. The pistons hit the lead-lithium from all sides, gen-erating a spherically symmetric pressure pulse that moves radially inward.We neglect the rotation of the lead-lithium in this setting. Although thephysical geometry is more complicated than this, a simple model like thisis a necessary step to reaching more complex models. Given the fact thatthe empty region of the cylinder will, in actual operation, become partiallyfilled in by the lead-lithium in compression, the real world system is closerto being spherical than it may seem.We treat the plasma as a spatially homogeneous spherical region, with a504.2. Governing Equationspressure given by the sum of a gas pressure and a magnetic pressure, withthe initial ratio of the gas to magnetic pressure of β0 = 0.1. We assumeadiabatic (reversible) conditions. Initially the pressure of the plasma andthe lead-lithium is 5 MPa, and the pressure pulse hits discontinuously, suchthat the pressure on the outer walls of the lead-lithium jumps to 2 GPa attime t = 0. The motion of the outer boundary between the lead-lithium andpistons, and the inner boundary between the lead-lithium and plasma havevelocities prescribed by the local fluid velocity of the lead-lithium. This isequivalent to the plasma and lead-lithium being immiscible and no lead-lithium leaving the confinement. We additionally neglect radiation lossesdue to material impurities.4.2 Governing EquationsIn this subsection, we present the equations used in our model. The numer-ical values of the parameters are tabulated in the following subsection.Within the lead-lithium, we use the compressible Euler equations forconservation of mass and conservation of momentum, equations (4.1) and(4.2) respectively. The variable ρ denotes mass density, v denotes the localfluid velocity, P denotes pressure, and t is the time. Note that there is atensor product ⊗.ρt +∇ • (ρv) = 0 (4.1)(ρv)t +∇ • (ρv ⊗ v) +∇P = 0 (4.2)In a spherically symmetric coordinate system, equations (4.1) and (4.2)514.2. Governing Equationsreduce to the equations used in our model, (4.3) and (4.4) respectively.ρt + (ρv)r +2r(ρv) = 0 (4.3)(ρv)t + (P + ρv2)r +2r(ρv2) = 0 (4.4)The lead-lithium pressure is given by equation (4.5)P = A(ρ− ρ0)2 +B(ρ− ρ0) + C (4.5)and the total pressure (magnetic plus gas pressure) is governed by equation(4.6)PL(t) = P (rL(t), t) =κ1rL(t)4+κ2rL(t)5. (4.6)The temperature is given asΘ = kT =κ2/rL(t)5n(t)(4.7)with k the Boltzmann constant, for a particle densityn(t) =n0(rL(t)/rL(0))3. (4.8)The pressure of the outer wall is given by equation (4.9)PR(t) = P (rR(t), t) = Patm + (Pmax − Patm)e−t2/t20 (4.9)The position of the lead-lithium and plasma interface rL and the outer wallof the lead-lithium rR are described by equations (4.10) and (4.11). As524.2. Governing EquationsFigure 4.1: Spherically Symmetric Modeldiscussed previously, these equations correspond to the boundaries movingwith the local fluid velocity.drLdt= v(rL(t), t) (4.10)drRdt= v(rR(t), t) (4.11)4.2.1 Parameters, Initializations, and SetupThe baseline parameters used in our study are given in table 4.1. We re-mark that our numerical computations explained in the following sectionwere done in a dimensionless framework. A figure illustrating our simplifiedgeometry is given in figure 4.1.534.3. Numerical ComputationsTable 4.1: A list of the parameters used in the numerical simulations.Parameter ValuePiston impact pressure PR(0) = Pmax 2 GPaInitial total pressure PL(0) 5 MPaInitial plasma particle density n0 3.2× 1016 particles per cm3Initial plasma temperature Θ0 100 eVInitial lead-lithium radius rR(0) 1.5 mInitial plasma radius rL(0) 0.2 mPulse time scale t0 45 µsInitial gas/magnetic pressure ratio β0 0.1Magnetic pressure coefficient κ1 7.273× 103 kg m3 s−2Gas pressure coefficient κ2 145.5 kg m4 s−2Lead-lithium parameter A 922.9 m5 kg−1 s−2Lead-lithium parameter B 4368100 m2 s−2Lead-lithium parameter C = Patm 101325 kg m−1 s−24.3 Numerical Computations4.3.1 Change of Variables to Map Problem to FixedDomainThe equations themselves are challenging to solve numerically due to themoving boundaries, the stiff equation of state, and the geometric sourceterms arising from the spherical geometry [23]. One strategy involves inter-polation and extrapolation while keeping the computational domain in theradial coordinate fixed [38]. We choose another approach here by transform-ing the coordinate system so that there is an underlying, constant compu-tational domain, without using interpolation or extrapolation. We performa change of variables to deal with the moving boundaries, setting τ = t andy = r−rL(t)∆(t) with ∆(t) = rR(t)− rL(t). We will denote Γ(t) by vR(t)− vL(t).Under this transformation, our new equations for mass and momentum544.3. Numerical Computationsconservation are given by equations (4.12) and (4.13)ρτ + {1∆[−(vL + Γy)ρ+ ρv]}y =−2rL + ∆yρv −Γ∆ρ (4.12)(ρv)τ + {1∆[−(vL + Γy)(ρv) + P + ρv2]}y =−2rL + ∆yρv2 −Γ∆ρv (4.13)in the fixed computational domain 0 ≤ y ≤ 1 and τ ≥ 0. Equations (4.12)and (4.13) are expressed in conservation form with time-dependent geomet-ric source terms on the right-hand sides.4.3.2 Finite Volume DiscretizationAfter the coordinate transformation, we use a finite volume discretizationfor the numerical simulations. The fundamental ideas we present here arepresented with excellent clarity and full detail in [23].We discretize the computational domain as follows. First, we choose anN and corresponding spatial mesh size h = 1/N , defining the coordinatesyi = ih with i = 0 and i = N on the boundary. We denote ρji to benumerical approximation to the average cell value∫ yi+h/2yi−h/2ρ(y, tj)dy at timetj . The index j denotes our time-stepping. Similarly we define (ρv)ji .The pseudocode of what is done at each time step is provided below:- Use constant extrapolation to define ρ and ρv at i = −1, 0 and i =N + 1, N + 2.- Compute the eigenvalues λ± =v−(vL+Γy)±√dP/dρ∆ of the linearizedsystems of 4.12 and 4.13.554.3. Numerical Computations- Obtain a time step k = 0.8h/λ∗ where λ∗ is the largest value |λ±|takes.- Using λ±i and λ±i+1, determine the direction of information propagationof the corresponding eigencomponents of the system by taking theiraverage.- Define flux limiters based on how rapidly changing the information isin the different directions.- Solve the homogeneous Riemann problem at each interface using thedirection information from the eigenvalues, but analytically precisenumerical flux.- Compute a second-order flux correction, and use the flux limiters todefine an updated set of values for ρ and ρv, again without the right-hand-side terms.- Advance forward once more in time using the nonzero right-hand sideas update information (a split step).- Update the time, boundary positions, and boundarypressures/densities.In the end, our method is first-order, despite the second-order correctionterm, due to the split stepping and the time-varying source terms. How-ever, this approach does accurately capture many aspects of the solutionsof conservation laws with shocks [23].564.3. Numerical ComputationsFigure 4.2: Profiles of velocity (left) and density (right) at various timevalues compared to the 1/rs growth rate predicted in the acoustic limit.4.3.3 Validation4.3.3.1 Linear Equation of State and Acoustic LimitIn the case of a linear equation of state for the lead-lithium A = 0 in (4.5),and a small-amplitude pressure pulse, the mass and momentum equationscan effectively be reduced to the acoustic equations. Under this reduction,the amplitude of the shock front in both the local fluid velocity and thepressure should grow inversely with the shock position rs. See [20] for details.In figure 4.2, we plot the pressure and velocity profiles at various timesalong with the curves predicted by a 1/rs growth rate. Our numerics arevalidated in this acoustic limit. Qualitative Convergence with N →∞We choose a fixed time of t = 1.243×10−4 s and plot the profiles of velocity,density, and pressure for N = 500 and N = 4000. The plots as displayed infigure 4.3 confirm convergence.574.3. Numerical ComputationsFigure 4.3: Profiles of velocity (left) and density (right) at t = 7.03× 10−4sfor N = 500 and N = 4000. The plots are very similar. First Order ConvergenceWe choose N = 250000 and define the predictions of this model as the“numerically exact” solution. Then we compare the errors of the velocityand density (normalized by their characteristic scales) in the L1−norms forN = 1000, N = 2000, N = 4000, N = 8000 and N = 16000 at the fixedtime t = 1.243× 10−4 s. In order to make this comparison, the N = 250000solution was projected onto the coarser mesh by a linear interpolation. Atthese same N -values, we computed the maximum error in total mass bynumerically integrating∫ rR(t)rL(t)r2ρ(r, t)dr and finding the maximum deviationfrom t = 0 for all times up to the fixed time. The errors are tabulated intable 4.2 and there is a very clear first-order trend. The convergence ratesbased on the slopes of the lines of best fit of the errors and N -values on alog-log plot are indicated in the table.584.4. Model PredictionsTable 4.2: Convergence of Numerical Scheme at t = 1.243× 10−4 sN Velocity Error Density Error Mass Error1000 5.68× 10−4 1.11× 10−3 4.79× 10−42000 2.95× 10−4 5.58× 10−4 2.39× 10−44000 1.39× 10−4 2.77× 10−4 1.20× 10−48000 6.97× 10−5 1.39× 10−4 5.97× 10−516000 3.51× 10−5 6.93× 10−5 2.99× 10−5Rate 1.01 1.00 1.004.4 Model Predictions4.4.1 Plasma Compression and Fusion PossibilityThe figures 4.4, 4.5, and 4.6 depict the position of the inner wall, its veloc-ity, and the pressure of the plasma as a function of time given the baselineparameters in table 4.1. With these baseline parameters, our triple productvalue is 0.52×1015 keV s cm−3, which is smaller than the Lawson triple prod-uct required for yielding more energy output than input. Under a modifiedparameter regime, as given in second last row of table 4.3, we obtain a tripleproduct that exceeds 4× 1015 keV s cm−3. We also make a brief remark onthe plasma survival time as survival times are an important aspect of theengineering. With our simulations, the plasma is compressed over the scaleof hundreds of microseconds, which is within the survival times of plasmassuch as spheromaks [10]. We thus infer that survival time should not be anissue. With these observations, it seems there is promise to yielding fusionenergy through such a design, although we acknowledge the rudimentarynature of the physical modelling considered.594.4. Model PredictionsFigure 4.4: The radius of the plasma vs time in the simulations.Figure 4.5: The radial velocity of the plasma-lead-lithium interface vs timein the simulations.604.4. Model PredictionsFigure 4.6: The total plasma pressure vs time in the simulations.4.4.2 Sensitivity AnalysisHere we analyze what results are obtained in varying some of the physicalparameters of the system. We concern ourselves with five aspects of thereactor performance: firstly, the minimal plasma radius rm; secondly, themaximum total plasma pressure PM; thirdly, the duration over which theplasma is compressed - the time over which the inner wall velocity is negativetc; fourthly, the maximum plasma temperature ΘM; and finally the Lawsontriple product computed asΠL =∫t:vL<0n(t)Θ(t)dt. (4.14)Table 4.3 provides the baseline value of the predictions, and subsequentlyvaries only one of the parameters.While many of the parameters have little effect upon the compression,some parameters cause drastic changes in the reactor behaviour, even whenperturbed by 10%, most notably the impulse pressure of the pistons, the614.4. Model Predictionsinitial size of the plasma, the initial size of the lead-lithium sphere, and thetimescale of the impulse. In interpreting these results, some care is needed,as not all the variables are physically independent. For example, the pressureintensity of the pistons is primarily governed by their impact velocity andthe sound speed of the lead-lithium. Increasing the pressure could come fromincreasing sound speed, perhaps by changing the lead-lithium to a differentmaterial, or by increasing the impact velocity. Observe, though, that purelyby increasing the sound speed at a fixed impulse pressure, the compressionis less effective.With a fixed initial plasma pressure, if the initial plasma radius is smallerthen there is more focusing and the compression is much more effective, justas it is with a higher impulse pressure or a larger lead-lithium initial radius.The time scale is also important; the more drawn out the pulse is, the moreenergy that can be added to the system to drive compression.As a point of interest, we consider the engineering parameters that canmost readily be controlled and the aspects of the model with the largest pos-sible error and perform larger perturbations. We consider the performance-enhancing prospects of doubling the pressure at the outer wall, doubling thepulse time scale, and doubling the initial radius of the lead sphere. Thedoubling of the pressure could in principle be achieved with faster pistonsand doubling the radius of the lead-lithium sphere is a possibility. There isuncertainty as to the precise form of the pressure imparted as a function oftime and doubling the time-scale, which essentially amounts to consideringa symmetric instead of one-sided Gaussian wave, is well within the realm ofmodelling uncertainty.Under optimal conditions, we obtain a Lawson triple product of 16×1015keV s cm−3, giving promise of net energy gain, which occurs when the624.5. Summary and Future WorkTable 4.3: A sensitivity analysis of the reactor performance upon varia-tions in the model parameters. We denote the minimum radius by rm, themaximum total pressure by PM, the compression time by tc, the maximumtemperature by θM, and the Lawson triple product by ΠL. The final threerows of the table illustrate the most promising conditions for fusion.Parameters rm (cm) PM (GPa) tc (ms) ΘM (keV) ΠL (1015 keV s cm−3)Baseline 3.6 6.7 0.43 2.7 0.52β0 ↓ 10% 3.6 6.8 0.43 2.5 0.49β0 ↑ 10% 3.6 6.6 0.43 2.9 0.55c∗ ↓ 10% 3.3 9.5 0.33 3.2 0.56c∗ ↑ 10% 3.8 5.3 0.52 2.5 0.46PR(0) ↓ 10% 4.4 2.8 0.63 1.8 0.37PR(0) ↑ 10% 3.0 15. 0.30 4.0 0.64PL(0) ↓ 10% 3.4 7.5 0.43 2.7 0.53PL(0) ↑ 10% 3.8 6.0 0.43 2.7 0.48rL(0) ↓ 10% 1.9 71. 0.20 8.0 1.4rL(0) ↑ 10% 5.5 1.7 0.85 1.4 0.25rR(0) ↓ 10% 4.3 3.0 0.61 1.9 0.27rR(0) ↑ 10% 3.0 13. 0.31 3.8 0.92t0 ↓ 10% 4.3 3.2 0.60 1.9 0.39t0 ↑ 10% 3.1 13. 0.32 3.7 0.61t0 ↑ 100% 1.3 570 0.15 20. 1.9PR(0) ↑ 100% 1.2 690 0.14 22 16rR(0) ↑ 100% 0.84 4900 0.12 24 2.5impact pressure is doubled. Doubling the pulse time or the lead-lithiuminitial radius also lead to large triple products, which are close to 4 × 1015keV s cm−3.4.5 Summary and Future WorkOur work here has developed the basis of a numerical framework that couldbe useful in fusion modelling, and with this numerical framework, we havemade a preliminary examination into the operational conditions and feasi-bility of magnetized target fusion reactors. Within the model assumptions,current designs for magnetized target fusion reactors are within reach of634.5. Summary and Future Workreaching fusion conditions. There is, however, a high degree of sensitivity ofphysical parameters upon the resultant behaviour. Having a higher impulsepressure, a large confining sphere, and a long impulse are key factors inachieving a highly efficient reactor.There are many avenues for future work. Firstly, some of the modelassumptions such as isentropic conditions could be relaxed, and we couldinclude more material properties of the plasma instead of it being homoge-neous, and consider the possible effects of turbulence. We could also considera geometry that is not spherically symmetric and examine how much theresults differ. Another natural question is how deviations from a sphericalsymmetry behave, and if there is a growth in such deviations, such as whathas been investigated in [37]. It would also be of great interest if an analyticestimate for the plasma compression could be obtained. We also note thatthis model’s predictions about the collapse differ somewhat from previousresults [18], except for some of the more promising results of the sensitivityanalysis. This model represents success in numerically solving conservationlaws with free boundaries with a constant computational domain, but at thesame time we realize that the physics has been greatly simplified. Furtherinvestigation, including tuning the mathematical model to more accuratelyrepresent the underlying physics, is warranted.64Chapter 5Formal Asymptotic AnalysisThis work on asymptotic analysis is taken from [25]. The physical system ispresented in chapter 1.2 and the model is derived in chapter 2. More detailsof the asymptotics including some slight extensions are found in section 6.2;due to page restrictions, they are not included in the paper below.5.1 Context for the AnalysisMathematical models of fluids often deal with nonlinear conservation lawssuch as the Euler or Navier Stokes equations. These equations are oftendifficult to solve or approximate analytically. In our work, we will be dealingwith the hyperbolic conservation laws for mass and momentum conservationgiven by the compressible Euler equations for fluids (see [8] for relevantbackground to these equations). There is an additional obstacle we dealwith, which are the free boundaries that arise naturally in our application.Free boundary, or moving boundary problems, entail a physical boundary(possibly the domain of a PDE or an interface) evolving with time due to theconditions at the boundary itself. Instead of the boundary position beingprescribed in the problem or being a known function of time, it must besolved for within the problem as its evolution depends upon the nature ofthe solution at any given time, where the solution evolves as a PDE in bothspace and time. A common example is that of a Stefan problem describing655.2. Physical Assumptionsthe position of the boundary of an ice cube as it melts: the velocity of theboundary depends on the jump in heat flux at the interface which changeswith time [30]. Another example, which is closely related to the problem ofinterest discussed here, relates to bubble dynamics and the motion of theboundary of a gas bubble in a fluid (see [5] for a thorough introduction tobubble dynamics). This work is novel in that asymptotic analysis allows usto overcome both the obstacle of the nonlinear conservation laws and thedifficulty of the free boundaries to yield insightful results for the application.In our context, the fluid will be lead-lithium between the free boundariesof an inner spherical plasma ball and an outer spherical shell impacted bypistons.Our objective is to predict the minimum radius of the plasma, as this isa key parameter in how efficient such a design could potentially be. The pa-rameter regime for the device winds up being at the boundary of asymptoticvalidity; we will furnish a minimum radius that is qualitatively consistentwith the numerics in its parameter sensitivity, but where the quantitativeaccuracy is a loose estimate. Our work also shows that the asymptotics arequantitatively accurate given a well-ordered parameter regime. Other fac-tors besides the minimum radius, such as confinement time, are not analyzedhere.5.2 Physical AssumptionsWe take the plasma to be spherical and we assume spherical symmetry forall components. Plasma pressures have both a gas pressure and magneticpressure component. In many plasmas, the magnetic pressure is much largerinitially. We model the plasma pressure by a magnetic pressure (scaling665.3. Equations and NondimensionalizationFigure 5.1: The geometry we consider is spherically symmetric, with rL andrR free boundaries.inversely with the fourth power of its radius) [40], neglecting the gas pressurecomponent.The pressure imparted by the pistons is described by a baseline pres-sure plus a very large pressure modulated by a Gaussian with a very shorttime-scale. The pressure of lead-lithium is taken to depend linearly on thedensity. Our model will neglect mixing effects and radiation losses. Thesesimplifications allow us to write more tractable equations and gain quali-tative insight into how the most basic engineering parameters influence thecompression of the plasma.5.3 Equations and NondimensionalizationHere we present the equations relevant to this study and nondimensionalizethe system so that formal asymptotics can be applied. The values andexplanations of the constants used are provided in tables 5.1 and 5.2. Adimensionless diagram of our model is given in figure 5.1.675.3. Equations and NondimensionalizationVariable Description ValueCs Sound speed in lead at 1 atm 2090 m/s%0 Density of lead at 1 atm 11340 kg/m3Patm Atmospheric pressure 101325 PaPplasma,0 Initial plasma pressure 5 MPaPmax Maximum piston pressure 2 GPaΓ Coefficient for magnetic pressure growth 8000 N m2T0 Guassian decay time scale 45 µsRinner,0 Initial plasma radius 0.2 mRouter,0 Initial lead sphere radius 1.5 mTable 5.1: Table of physical constants/parameters.5.3.1 Euler EquationsWithin the lead-lithium, we work with the Euler equations for mass andmomentum conservation. Denoting the mass density by %, the local fluidvelocity by V and the pressure by P in a space-time coordinate system(X, T ), we have:%T +∇X • (%V) = 0 and (%V)T +∇XP +∇X · (%V ⊗V) = 0which with spherical symmetry in space-time coordinates (R, T ), where thevector quantity V is replaced by V , the radial velocity, reduce to%T +(%V )R+2R(%V ) = 0 and (%V )T +PR+(%V2)R+2R(%V 2) = 0. (5.1)From experimental data, the sound speed Cs and density %0 of lead atatmospheric pressure Patm are known [34]. The square of the sound speedis the derivative of P with respect to density % [39] allowing us to express685.3. Equations and NondimensionalizationP (ρ) as an approximate linear function:P (%) = C2s (%− %0) + Patm. (5.2)We approximate the equation of state of lead-lithium by the linearized equa-tion of state for lead and permit negative absolute pressures. Negative ab-solute pressures are observed over short time scales in some materials [13].5.3.2 Boundary ConditionsDue to the impact of the pistons and the interaction of the wave with theplasma boundary, this is a moving boundary problem. We denote RL to bethe radius of the inner wall of the lead-lithium sphere (the plasma radius),and RR to be the radius of the outer wall of the lead-lithium sphere. Thelocal fluid velocity at the inner and outer boundaries matches the velocitiesof these walls respectively:VL(T ) ≡dRLdT= V (RL(T ), T ) and VR(T ) ≡dRRdT= V (RR(T ), T ).(5.3)At any given time, the lead-lithium occupies the spherical shell between RLand RR.At the inner wall, the pressure is given by the plasma magnetic pressure.This givesPL(T ) = P (RL(T ), T ) =ΓRL(T )4=Pplasma,0(RL(T )/RL(−∞))4(5.4)where Γ is a constant chosen so that with PL = ΓR−4L , PL(RL(−∞)) =Pplasma,0. This is an algebraic coupling of the fluid pressure and plasma695.3. Equations and Nondimensionalizationpressure at this wall. At the outer wall, the pressure is given by the pistonpressure which is modelled by a Gaussian pressure in addition to a baselineminimum applied pressurePR(T ) = (Pmax − Pplasma,0)e−T 2/T 20 + Pplasma,0. (5.5)5.3.3 Initial ConditionsThe maximum piston pressure occurs at T = 0, but we allow for negativetimes relative to this. The system begins at T = −∞ whereV (R,−∞) = 0, P (R,−∞) =ΓRL(−∞)4, %(R,−∞) =P (R,−∞)− Pplasma,0C2s+ %0(5.6)withRL(−∞) = Rinner,0, RR(−∞) = Router,0. (5.7)5.3.4 NondimensionalizationTo nondimensionalize, we will look at (5.1) with the change of variables% = %¯ρ, V = V¯ v, P = P¯ p, R = R¯r, T = T¯ t where the bars denotecharacteristic dimensional quantities and ρ, v, p, r, and t are dimensionless.We obtain:%¯T¯ρt+%¯V¯R¯(ρv)r+%¯V¯R¯2r(ρv) = 0 and%¯V¯T¯(ρv)t+P¯R¯pr+%¯V¯ 2R¯(ρv2)r+%¯V¯ 2R¯2r(ρv2) = 0.(5.8)By matching dimensional terms in (5.8)1, we have %¯/T¯ = %¯V¯ /R¯ so thatV¯ = R¯/T¯ . (5.9)705.3. Equations and NondimensionalizationVariable Notation/Equation Characteristic ScaleDensity % %¯ = %0 11, 340 kg m−3Pressure P P¯ = Pmax 2 GPaDistance R R¯ = Router,0 1.5 mTime T T¯ = R¯√%¯/P¯ 3.57 msRadial Velocity V V¯ =√P¯ /%¯ 420. m s−1Table 5.2: Characteristic dimensions for the system. The first three entriesare chosen by the operating conditions. The last two entires can be derivedfrom the first three.By matching dimensional terms in (5.8)2, we have (%¯V¯ )/T¯ = P¯ /R¯ giving%¯ = (P¯ T¯ )/(R¯V¯ ) (5.10)and (%¯V¯ )/T¯ = (%¯V¯ 2)/R¯ giving the same as in (5.9).Equations (5.9) and (5.10) give good guidance as to reasonable scalingsfor the system, and we also have some freedom. We will choose R¯ = 1.5m, the initial outer radius; %¯ = 11340 kg/m3, the density of lead at at-mospheric pressure; P¯ = 2 GPa, the maximum pressure imparted on theouter wall. Substituting (5.9) into 5.10, we have % = (P¯ T¯ 2)/(R¯2) so thatT¯ =√R¯2%¯/P¯ = 3.57 ms. Using this now known value of T¯ in (5.9) givesV¯ = 420. m/s. Table 5.2 summarizes this.5.3.5 Selecting an Asymptotic ParameterThrough the nondimensionalization, various dimensionless parameters ap-pear. Upon considering the physics of the problem, an intense pressureimparted over a very short time-scale (45µs  3.57 ms), we will define = T0/T¯ = 0.0126 such that the Gaussian component of the pressurebecomes e−(tT¯ /T0)2= e−t2/2 . With respect to this value of , the differ-715.3. Equations and Nondimensionalizationent dimensionless parameters have characteristic orders. For example, indimensionless form, p = c2(ρ − 1) + d where c2 = C2s %¯/P¯ ≈ 24.7 andd = Patm/P¯ ≈ 5.07 × 10−5. The density of lead-lithium corresponding tothe maximum pressure is 1.04, which is well described by ρ = 1 + O().If ρ = 1 + O() corresponds to p = O(1) then c2 = O(−1). Thus, it isreasonable to write p(ρ) = b2 (ρ− 1) + a5/2 where b = 0.557 and a = 2.84.In fact, since we are aiming for a leading-order behaviour of the system, wewill neglect a5/2 entirely since it is so small. Also due to its negligible size,O(3/2), we neglect the dimensionless Pplasma,0 term in the pressure condi-tion at the right wall. Discarding these terms can be justified a posteriori bynoting that in all the analysis to come, only pressures that are O() balancein the relevant equations.Proceeding in a similar way with the remaining equations and boundaryconditions, and making use of the fact that knowing either p or ρ is equiv-alent, the final set of equations and conditions we consider are presentedbelow, valid for rL < r < rR, −∞ < t < ∞. In due course, all of theseequations will be used.ρt + (ρv)r +2r(ρv) = 0, (ρv)t +b2ρr + (ρv2)r +2r(ρv2) = 0 (5.11)p =b2(ρ− 1) (5.12)vL(t) =drLdt= v(rL(t), t), vR(t) =drRdt= v(rR(t), t) (5.13)725.4. Formal Asymptotic AnalysisVariable Meaning Formula Value pulse time scale T0R¯√%¯/P¯0.0126b scaled sound speed Cs√√P¯ /%¯0.557γ scaled magnetic pressure coefficientPplasma,0R4inner,0P¯ R¯47/23.52χ scaled initial inner radius Rinner,0R¯√1.19Table 5.3: The dimensionless parameters of the asymptotic model.pL(t) = p(rL(t), t) =γ7/2rL(t)4, pR(t) = p(rR(t), t) = e−t2/2 (5.14)v(r,−∞) = 0, p(r,−∞) =γ7/2rL(−∞)4, ρ(r,−∞) = 1 (5.15)rL(−∞) = χ1/2, rR(−∞) = 1 (5.16)We summarize the constants in the table 5.3. Our main objective is toestimate the smallest value of rL, which we denote r∗.5.4 Formal Asymptotic AnalysisThe analysis has five distinct phases: I, formation; II, focusing; III, reflect-ing; IV, slow collapse; and V, maximum compression. See figure 5.2. Muchof the analysis is done on the linear acoustic equations. Techniques include:Riemann invariants, the velocity potential formulation with a transformationthat turns the linear acoustic equations for a sphere into one-dimensionalwave equations, indirect computations of long-term behaviours by integralsof boundary conditions, energy arguments, and matched asymptotics onthe Rayleigh-Plesset equation. We ultimately compute the plasma radius at735.4. Formal Asymptotic AnalysisFigure 5.2: Relevant space-time scales and boundary motion. In phase I,pulses are formed, moving towards the plasma. In phase II, the pulses moveradially inward, growing in amplitude due to focusing. In phase III, thepulses interact with the plasma with much of the energy reflecting but asmall portion of useful energy remaining. In phase IV, the plasma slowlybegins to collapse and in phase V, the maximum compression is reached.which the velocity of the plasma boundary switches from negative (i.e. com-pressing) to zero at leading order, which we define as the minimum plasmaradius. Effectively I is an inner region that is matched to II; III is an innerregion that can match to II and which matches to IV over a long time; andV is an inner region to IV.5.4.1 Phase I5.4.1.1 Pulse Formation SetupTo analyze the formation of the pulse, we rescale time so that t = τ, andmake the ansatzes that ρ ∼ 1 + ρ1 + `ρ2 and v ∼ m0v0 + m1v1. Seeingas the peak impulse pressure of p = O(1) corresponds to ρ = 1 + O()745.4. Formal Asymptotic Analysisthese expansions are justified. In order to capture the pulse formation, wealso define y = r−1n . This leads to ρv ∼ m0v0 + m1v1 + m0+1ρ1v0 andρv2 ∼ 2m0v20, r = 1 + 1/2y, 2/r ∼ 2 − 21/2y, ∂r = −n∂y and ∂t = −1∂τ .Using this in (5.11), we have to various leading orders:−1(ρ1 + `ρ2)τ + −n(m0v0 + m1v1 + m0+1ρ1v0)y+ 2(m0v0 + m1v1 + m0+1ρ1v0) = 0−1(m0v0 + m1v1 + m0+1ρ1v0)τ + −n−1b2(ρ1 + `ρ2)y + −1/2(2m0v20)y+ 2(2m0v20) = 0.Balancing as many terms as possible to leading order, we choose m0−n = 0and m0 − 1 = −n so that m0 = n = 1/2. This givesρ1,τ + v0,y + `−1ρ2τ + m1−1/2v1,y + 21/2v0 = 0−1/2(v0,τ + b2ρ1,y) + m1−1v1,τ + `−3/2ρ2,y + 1/2(v20)y = 0.We now balance as many terms as possible at the next order with ` − 1 =m1 − 1/2 = 1/2 so that ` = 3/2, and m1 = 1. We thus arrive atρ ∼ 1 + ρ1 + 3/2ρ2, v ∼ 1/2v0 + v1with t = τ and y = r−11/2.To deal with boundary conditions, we first need to find an asymptotic ex-pression for rR ∼ 1 + krR0(τ). To leading order vR =drRdt = O(1/2). Thus,r′R(t) = kr′R0(t) = O(1/2). Using d/dt = −1d/dτ , we have k−1r′R0(τ) =755.4. Formal Asymptotic AnalysisO(1/2) so that k = 3/2.It is worth mentioning that when we write rR0(τ) we really mean rR0(t =τ) and not rR0(t = τ)! Throughout this paper we use the shorthand thata dependent variable evaluated at a rescaled variable should be interpretedby the proper scalings in the original (r, t) variables.The outer wall changes by a distance O(3/2). If we denote yR as theposition of the outer boundary in the y−coordinates, then yR =rR−11/2=rR0(τ) + .... The boundary condition that must be upheld is pR(τ) =b2 (ρ− 1) = e−τ2 so thatρ1(rR(τ), τ) + 1/2ρ2(rR(τ), τ) =1b2e−τ2.We can now perform a Taylor expansion (seeing as the problem is linearand the forcing terms are continuous and differentiable). Keeping just a fewterms, we obtain b2(ρ1(0, τ) + ρ1y(0, τ)(rR0) + 1/2ρ2(0, τ)) = e−τ2whichimpliesρ1(0, τ) =1b2e−τ2and ρ2(0, τ) = 0. (5.17)We obtain the following initial conditions in this scaling:ρ1(y,−∞) = ρ2(y,−∞) = v0(y,−∞) = v1(y,−∞) = 0. (5.18)Figure 5.3 depicts this phase.765.4. Formal Asymptotic AnalysisFigure 5.3: Depiction of phase I. This is the formation phase. In this phase,ρ ∼ 1 + ρ1 + 3/2ρ2 and v ∼ 1/2v0 + v1. Riemann InvariantsBoth the leading order and first correction terms can be written in the formρτvτ+0 1b2 0ρyvy =f(y, τ)g(y, τ) (5.19)and such problems can be solved effectively with a technique known as Rie-mann Invariants [6].The 2 × 2 matrix in (5.19) has eigenvalues ±b with corresponding lefteigenvectors `± = (±b, 1). Upon left multiplying (5.19) by `±, we obtain(±bρ+ v)τ ± b(±bρ+ v)y = ±bf + g. Thus, the system has characteristics.Denoting c± = ±bρ+v, we have that along dy/dτ = ±b, dc±/dτ = ±f+g ≡h±. Given that ρ = v = 0 at τ = −∞, (so c± = 0 at τ = −∞) we cancompute c+(y, τ) = 0+∫ τ−∞ h(y˜(τ˜), τ˜)dτ˜ where y˜(τ˜) = y+b(τ˜−τ) describesthe straight pathway along the rightgoing characteristic starting at τ = −∞and reaching (y, τ). Thus,c+(y, τ) =∫ τ−∞h+(y + b(τ˜ − τ), τ˜)dτ˜ . (5.20)775.4. Formal Asymptotic AnalysisFigure 5.4: Illustration of the right- and left-going characteristic paths andparameterizations.To find c−(y, τ), we need to know its value at the point along the τ−axiswhose leftgoing characteristics eventually reach (y, τ). As c+ is known alongthe τ−axis, and ρ is also specified there, we can find c−. The leftgoingcharacteristic path leading to (y, τ) can be described by (y˜(τ˜), τ˜) = (b(τ −τ˜), τ˜) for τ+y/b ≤ τ˜ ≤ τ. From c+(0, τ+y/b) = bρ(0, τ+y/b)+v(0, τ+y/b)we get c−(0, τ + y/b) = −bρ(0, τ + y/b) + v(0, τ + y/b) = c+(0, τ + y/b) −2bρ(0, τ + y/b). We thus havec−(y, τ) = c+(0, τ +y/b)−2bρ(0, τ +y/b) +∫ ττ+y/bh−(b(τ − τ˜), τ˜)dτ˜ (5.21)Given c±, we can recoverρ =12b(c+ − c−) and v =12(c+ + c−). (5.22)Figure 5.4 illustrates the characteristic curves.785.4. Formal Asymptotic Analysis5.4.1.3 Leading OrderAt leading order, using (5.18) and (5.17)1 we haveρ1,τv0,τ+0 1b2 0ρ1,yv0,y = 0,ρ(0, τ) = 1b2 e−τ2 , ρ(y,−∞) = v(y,−∞) = 0. Thus, c+(y, τ) = 0 by (5.20)since the forcing function h = 0. And c−(y, τ) = −2b 1b2 e−(τ+y/b)2= −2be−(τ+y/b)2 by (5.21). We then use (5.22) to getρ1(y, τ) =1b2e−(τ+y/b)2(5.23)andv0(y, τ) =−1be−(τ+y/b)2. (5.24) CorrectionAt the next order, we obtainρ2τv1τ+0 1b2 0ρ2yv1y =−2v00subject to ρ2(0, τ) = 0, ρ2(y,−∞) = v1(y,−∞) = 0 as per (5.17)2 and(5.18). In this case, h± = ∓2bv0. By (5.20),c+(y, τ) =∫ τ−∞(−2b)v0(y + b(τ˜ − τ), τ˜)dτ˜ =∫ τ−∞2e−(2τ˜+y/b−τ)2dτ˜ .795.4. Formal Asymptotic AnalysisWe can evaluate this integral with a substitution u = 2τ˜ + y/b− τ to getc+ =∫ τ+y/b−∞e−u2du =∫ ∞−τ−y/be−u2du =√pi2erfc(−τ − y/b).Similarly, by (5.21),c−(y, τ) =√pi2erfc(−τ − y/b)− 0 +∫ ττ+y/b(2b)v0(b(τ − τ˜), τ˜)dτ˜=√pi2erfc(−τ − y/b) +2ybe−(τ+y/b)2.We once again use (5.22) to furnishρ2(y, τ) =−yb2e−(τ+y/b)2and v1(y, τ) =√pi2erfc(−τ − y/b) +ybe−(τ+y/b) Phase I ResultsTo leading order, the density and velocity are described by left-going planewaves. Due to the scalings, we have been looking on such a small scalethat the system has not quite “noticed” the spherical nature of the problemand everything appears flat (for y = O(1)). The divergent terms with yprefactors, however, actually result from the spherical geometry as we shallsee in the following section.We haveρ ∼ 1 + 1b2e−(τ+y/b)2− 3/2yb2e−(τ+y/b)2(5.25)andv ∼ 1/2−1be−(τ+y/b)2+ (√pi2erfc(−τ − y/b) +ybe−(τ+y/b)2). (5.26)A plot of these solutions at τ = 1/2 is given in figure 5.5. In the r-805.4. Formal Asymptotic Analysis0.8 0.85 0.9 0.95 1−0.2−0.15−0.1−0.050rVelocity  vnumε0.5 v0ε0.5 v0 + ε v10.8 0.85 0.9 0.95  ρnum1+ε ρ11 + ε ρ1 + ε1.5 ρ20.95 0.96 0.97 0.98 0.99 1−0.03−0.025−0.02−0.015rVelocity  vnumε0.5 v0 ε0.5 v0 + ε v10.95 0.96 0.97 0.98 0.99 11.00051.00061.00071.00081.00091.0011.0011rDensity  ρnum1+ε ρ11+ε ρ1 + ε1.5 ρ2Figure 5.5: Plots depicting the profiles of the velocity and density at avery small value of t to validate the asymptotics of phase I. In the case ofthe second row, the two term asymptotic expansions agree so well with thenumerics the plots cannot be distinguished. Parameters: b = 0.557, χ =1.19, γ = 3.52,  = 0.0126, t = /2 (top row); b = χ = γ = 1,  = 0.001, t =/2 (bottom row).coordinates, the characteristic widths of these pulses are O(√). The plotincludes numerical validation of the results which are discussed in more de-tail in section 5.5.815.4. Formal Asymptotic Analysis5.4.2 Phases II and III5.4.2.1 Spherical Focusing SetupHaving obtained the basic shape of the pulses being formed, we can nowanalyze how these pulses behave as they move towards the centre of thesphere. See figure 5.6. From phase I, we have ρ ∼ 1 + ρ1 + 3/2ρ2 andv ∼ 1/2v0 + v1 and we shall make a similar ansatz here. We will considerr = O(1) and choose a time scale t = kT . With these scalings, (5.11)becomes:−k(ρ1 + 3/2ρ2)T + (1/2v0 + v1)r +2r(1/2v0 + v1) = 0−k(1/2v0 + v1)T + −1b2(ρ1 + 3/2ρ2)r + (v20)r +2r(v20) = 0.Clearly k = 1/2 yields a suitable balance. The equations are:ρ1,T + v0,r +2rv0 = 0 and v0,T + b2ρ1,r = 0where the identical equations are upheld for ρ2 and v1.The initial and boundary conditions are more subtle here. They cannotbe written down clearly because of the different time scales (typically, tomatch at T = 0, we would require τ → ∞ but in this case, the solutionsfound in phase I tend to 0). However, by solving the hyperbolic systemsat different orders, we can find the functional form of the solution, whichallows us to find the solutions by matched asymptotics.825.4. Formal Asymptotic AnalysisFigure 5.6: Depiction of phase II. This is the focusing phase. In this phase,we initially have ρ ∼ 1+ ρ1 + 3/2ρ2 and v ∼ 1/2v0 + v1, and finally obtainρ ∼ 1 + 1/2ρ1 + ρ2 and v ∼ v0 + 1/2v1. Spherical Linear Acoustic LimitFor the two leading orders, the systems of PDEs that need to be solved areof the formρT + vr +2rv = f and vT + b2ρr = g. (5.27)Such equations also appear in linear acoustic problems, often with f = g = 0(and in our case we will not need to explicitly solve such an equation withf, g 6= 0).To solve this, we will make use of a velocity potentialφ(r, T ) =∫ r∞ v(s, T )ds. With this substitution, we haveρT + φrr +2rφr = f and φrT + b2ρr = g. (5.28)Integrating (5.28)2 from r = ∞, assuming ρ = 0 at r = ∞, we haveφT (r, T ) + b2ρ(r, T ) =∫ r∞ g(s, T )ds ≡ G(r, T ) so that φTT = −b2ρT + GT .835.4. Formal Asymptotic AnalysisUsing (5.28)1 here now gives φTT = −b2(−φrr − 2rφr + f) +GT soφTT = b2(φrr +2rφr) +GT − b2f. (5.29)Note that by taking a time derivative of (5.27)1 and a space derivativeof (5.27)2 we haveρTT = −vrT −2rvT + fT and vrT = gr − b2ρrr. (5.30)Equations (5.30)2 and (5.30)1 can then be used in (5.27)1 to obtain ρTT =−(gr − b2ρrr)− 2r (g − b2ρr) + fT so thatρTT = b2(ρrr +2rρr) + fT −2rg − gr. (5.31)A symmetry argument now helps to solve (5.29) and (5.31). If φ = Φ/rand ρ = K/r then φrr + 2rφr =Φrrr and similarly for K. Thus we cantransform these equations into one-dimensional wave equations [39]:ΦTT = b2Φrr + S1(r, T ), KTT = b2Krr + S2(r, T ) (5.32)with S1(r, T ) = r(GT − b2f), and S2(r, T ) = r(fT − 2rg − gr). Leading Order and Correction in Outer RegionThe system of equations for ρ1 and v0 and ρ2 and v1 are the same, and inboth cases the source terms are fortunately zero. With no source terms,and with only leftgoing waves (since there are no sources going right), thesolutions of (5.32) are K = P (r+bT ) and Φ = Q(r+bT ) respectively, whereP here is not the dimensional pressure. Then, being mindful of the fact that845.4. Formal Asymptotic Analysisthe pulses generated have a width of size O(√) and are leftgoing originatingnear r = 1, we can write the solutions to the systems of equations as ρ1 =1rP1(r−1+bT√ ), ρ2 =1rP2(r−1+bT√ ), v0 =1rQ′0(r−1+bT√ ) −√r2 Q0(r−1+bT√ ), andv1 = 1rQ′1(r−1+bT√ )−√r2 Q1(r−1+bT√ ).We will now take these solutions and express them in the inner (y, τ)-coordinates. We will apply the Van-Dyke matching of inner and outer so-lutions [14], and since we have a two-term expansion for the inner solution,we also take two terms in the outer solution. We have t = τ = 1/2T soT = 1/2τ and r = 1 + 1/2y so that r−1+bT√ = y + bτ . This allows us towrite:ρ1 =P˜1(y + bτ)1 + 1/2y∼ P˜1(y + bτ)− 1/2yP˜1(y + bτ)ρ2 =P˜2(y + bτ)1 + 1/2y∼ P˜2(y + bτ)− 1/2yP˜2(y + bτ),v0 =Q˜0′(y + bτ)1 + 1/2y−1/2Q˜0(y + bτ)(1 + 1/2y)2∼ Q˜0′(y+bτ)+1/2(−Q˜0′(y+bτ)y−Q˜0(y+bτ)),v1 =Q˜1′(y + bτ)1 + 1/2y−1/2Q˜1(y + bτ)(1 + 1/2y)2∼ Q˜1′(y+bτ)+1/2(−Q˜1′(y+bτ)y−Q˜1(y+bτ))Now we write(1 + ρ1 + 3/2ρ2)outer = (1 + ρ1 + 3/2ρ2)inner1 + (1b2e−(τ+y/b)2) + 3/2(−yb2e−(τ+y/b)2) = 1 + (P˜1(y + bτ))+3/2(−yP˜1(y + bτ) + P˜2(y + bτ))so that P˜1(x) = 1b2 e−(x/b)2 and P˜2(x) = 0. Doing the same for the velocity,855.4. Formal Asymptotic Analysiswe write:(1/2v0 + v1)outer = (1/2v0 + v1)inner1/2−1be−(τ+y/b)2+ (√pi2erfc(−τ − y/b) +ybe−(τ+y/b)2) = 1/2Q˜0′(y + bτ)+(−yQ˜0′(y + bτ)− Q˜0(y + bτ) + Q˜1′(y + bτ))giving us Q˜0(x) = −√pi2 erfc(x/b) and Q˜1′(x) = 0 so that Q˜1 is constant.Since Pi(x/√) = P˜i(x) and Qi(x/√) = Q˜i(x), i = 0, 1, we have thatduring this initial phase of focusing,ρ ∼ 1 + (1b2re−(r−1+bTb√)2) + o(3/2) (5.33)andv ∼ 1/2(−1bre−(r−1+bTb√)2) + (√pi2r2erfc(−r − 1 + bTb√)) + o(). (5.34) Leading Order and Correction in Inner RegionAs the pulses move inwards, due to the 1/r and 1/r2 terms, the amplitudesgrow and the asymptotic expansions we arrived at above lose their validity.To describe the region where the pulses have grown and arrive at the innerwall, we rescale space and time. We let r = 1/2σ and t − ts = tˆ, wherets is a shifting parameter such that the peak of the leftgoing pressure wavewill have reached the plasma wall at r = χ1/2 at tˆ = 0. We anticipatedifferent scalings for the asymptotic expansions. From r = O(1) to r =O(√), the first perturbation term in the density perturbation given in (5.33)should have grown from O() to O(√). Similarly, from (5.34), the velocityamplitude should have grown from O(√) to O(1) and the velocity term865.4. Formal Asymptotic Analysisthat was O() in that expression should also become O(1) since it grows like1/r2. Therefore, in this next inner region, we write ρ ∼ 1 + 1/2ρ1 + ρ2 andv ∼ v0 + 1/2v1. Using ∂r = −1/2∂σ and ∂t = −1∂tˆ, (5.11) becomes−1(1/2ρ1 + ρ2)tˆ + −1/2(v0 + 1/2(v1 + ρ1v0))σ + −1/2 2σ(v0 + 1/2(v1 + ρ1v0)) = 0−1(v0 +1/2(v1 +ρ1v0))tˆ+−3/2b2(1/2ρ1 +ρ2)σ+−1/2(v20)σ+−1/2 2σv20 = 0which upon looking at balancing terms yield two sets of PDEs:ρ1,tˆ + v0,σ +2σv0 = 0, v0,tˆ + b2ρ1,σ = 0 (5.35)andρ2,tˆ+v1,σ+2σv1 = −(ρ1v0)σ−2σρ1v0, v1,tˆ+b2ρ2,σ = −(ρ1v0)tˆ−(v20)σ−2σv20.(5.36)By virtue of the fact these PDEs can be solved and matched to the outersolution (as done in the following paragraphs), we can be confident in thescalings chosen. We will now consider boundary conditions. For a timescalethat is O(), given a velocity v = O(1), the inner wall can only move adistance O() and thus the inner wall position remains O(√) (recall thatrL(−∞) = χ√). We can find the leading order nonzero contribution to theinner wall motion in the (σ, tˆ)-coordinates. We define σL = rL/√. ThenσL(tˆ) = χ+∫ tˆ−∞√v(σL(sˆ), sˆ)dsˆ. The factor of√ comes from the fact thatdrLdt = v(rL(t), t) =d1/2σLdtˆ= −1/2 dσLdtˆ. Thus,σL(tˆ) = χ+∫ tˆ−∞√v(χ+O(√), sˆ)dsˆ = χ+ 1/2∫ tˆ−∞v(χ, sˆ)dsˆ+O().875.4. Formal Asymptotic AnalysisAs the pressure at the inner wall is given by p = γ7/2/r4L, for rL =O(√), p = O(3/2) = b2 (1/2ρ1(σL, tˆ) + ρ2(σL, tˆ) + ...) so that by a Taylorexpansion, 1/2ρ1(χ, tˆ)+(ρ1σ(χ, tˆ)∫ tˆ−∞ v0(χ, sˆ)dsˆ+ρ2(χ, tˆ))+o() = 0. Thisgives us boundary conditionsρ1(χ, tˆ) = 0 and ρ2(χ, tˆ) = −ρ1σ(χ, tˆ)∫ tˆ−∞v0(χ, sˆ)dsˆ. (5.37)If we again consider (5.35) in potential form with φ0(σ, tˆ) =∫ σ∞ v0(σˆ, tˆ)dσˆthen by integrating the equation from σ =∞ to σ = χ, φ0tˆ(χ, tˆ)+b2ρ1(χ, tˆ) =0 where we used ρ1 = φ0 = 0 at σ =∞. Since ρ1(χ, tˆ) = 0, we must haveφ0tˆ(χ, tˆ) = 0. (5.38)Similarly by integrating (5.36) we have φ1tˆ(χ, tˆ)+b2ρ2(χ, tˆ) =∫ χ∞(−(ρ1v0)tˆ−(v20)σ −2σv20)dσ which with (5.37)2 givesφ1tˆ(χ, tˆ) = b2ρ1σ(χ, tˆ)∫ tˆ−∞v0(χ, sˆ)dsˆ+∫ χ∞(−(ρ1v0)tˆ − (v20)σ −2σv20)dσ.(5.39)We can find the form of the solution to (5.35) using the solutions of(5.32) with S1 = S2 = 0. We note that there are both incoming signals(from σ =∞) and outgoing signals (from the interaction of the pulses withthe plasma) so we have waves propagating in both directions. Also, we areon a spatial scale of O(√) so that it is not necessary to rescale the argumentof the solutions like we did in the previous section. As tˆ → −∞, we willneglect the rightgoing wave and consider only the incoming wave. We writeρ1 =P+1 (σ+btˆ)σ , v0 =Q+0′(σ+btˆ)σ −Q+0 (σ+btˆ)σ2 .Expressing the outer solutions in (r, T )-coordinates in the inner coordi-885.4. Formal Asymptotic Analysisnates and going to leading order, and using the fact that any function ofr−1+bTb√ is a function ofσ−χb + tˆ, we have(1 + 1/2ρ1)inner = (1 + 1/2ρ1)outer1 + 1/2P˜+1 (σ−χb + tˆ)σ= 1 + (1b2σ1/2e−(σ−χb +tˆ)2)giving us a leading order solution in the inner region for the densityρ1 =1b2σe−(σ−χb +tˆ)2(tˆ→ −∞). (5.40)Also,(v0)inner = (1/2v0 + v1)outerQ˜+0′(σ−χb + tˆ)σ−Q˜+0 (σ−χb + tˆ)σ2= 1/2(−1bσ1/2e−(σ−χb +tˆ)2) + (√pi2σ2erfc(σ − χb+ tˆ))so that v0 = −1bσ e−(σ−χb +tˆ)2+√pi2σ2 erfc(−(σ−χb + tˆ)).These solutions are only valid for tˆ → −∞. To get solutions for othertimes, we first find the potential as φ0 =∫ σ∞ v0(σ, tˆ)dσ so thatφ0 = −√pi2σ erfc(−σ−χb − tˆ) (tˆ → −∞). From equation (5.37)1, in order forρ1 to be 0 along σ = χ and to have the incoming solution given in (5.40),the solution isρ1 =1b2σ(e−(σ−χb +tˆ)2− e−(χ−σb +tˆ)2) (5.41)For the potential, we note that φ0(χ,−∞) = −√pi2χ erfc(∞) = 0 so thatby (5.38), φ0(χ, tˆ) = 0 is the boundary condition for all tˆ. Using this, andthe wave solution to φ0, we must have thatφ0 =√piσ(erfc(−χ− σb− tˆ)− erfc(−σ − χb− tˆ))895.4. Formal Asymptotic Analysissov0 =√pi2σ2(erf(σ − χb− tˆ)− erf(−χ− σb− tˆ))−1bσ(e−(σ−χb +tˆ)2+ e−(χ−σb +tˆ)2)(5.42)In writing the above equation, we used the fact that erfc(−x− y)− erfc(x−y) = (1− erf(−x− y))− (1− erf(x− y)) = erf(x− y)− erf(−x− y).We note that the leading order density perturbation and velocity profilesgrow in amplitude like 1/r, or close to such scaling, at their peak values asthe profiles move radially inwards. This is consistent with known resultsfor fluid dynamics with spherical symmetry [20]. Given the equation ofstate (5.12), this also means the peak pressure grows like 1/r. To be moreprecise, the asymptotic profile for ρ1 and p do grow like 1/r, but becausev0 in the inner region also includes an error function term that reduces themagnitude of the Gaussian peak, the scaling may be slightly smaller. Atypical profile of the velocity and pressure, computed numerically [27] isdepicted in figure 5.7. We also observe that v0 and ρ1 tend to 0 as tˆ → ∞and the total distance the plasma wall moves due to these leading orderterms is O(1)×O() = O(). More precisely the leading order displacementsolely due to the v0 term isδ0 = ∫ ∞−∞v0(χ, tˆ)dtˆ = ∫ ∞−∞−2bχe−tˆ2dtˆ = −2√pibχ. (5.43)Any nonnegligible compression will be a result of the correction terms. Equa-tions (5.36), (5.37)2, and (5.39) give insights into what the effects of thesecorrection terms are, which we discuss in the following portion of this paper.905.4. Formal Asymptotic AnalysisFigure 5.7: The pulses are moving to the left. The pressure peak growslike 1/r as the pulses move inwards. The velocity has a peak with negativevalue and one with positive value. The negative value grows roughly like1/r. Parameters: b = χ = γ = 1,  = 0.0025,t = 0.0262. Phase III - Resultant Velocity FieldFigure 5.8 depicts this region of analysis.Given that most of the localized disturbances are reflected, the drivingforce for compression, if a substantial compression is to occur, must comefrom a residual negative radial velocity of the inner wall after the wavescome in. We begin by considering the systems given by (5.36) in a time-independent regime (a long time after the v0 and ρ1 pulses have interactedwith the inner wall). We then find that v1σ + 2σv1 = 0 and ρ2σ = 0 so thatv1 = A/σ2 and ρ2 = B. Note that these A and B values have nothing todo with equation (4.5). Physically, ρ2(σ = ∞) = 0 so that B = 0, butthe constant A cannot be determined. We can, however determine φ1 =∫ σ∞Aσˆ2 dσˆ = −A/σ (tˆ → ∞). This suggests that the correction term leavesa velocity field behind that may do further work to compress the plasma.Strictly speaking, we show that the value of φ1 approaches a constant atσ = χ as tˆ → ∞ but section 5.5.3 validates this residual velocity field915.4. Formal Asymptotic AnalysisFigure 5.8: Depiction of phase III. This is the reflection phase. In this phase,we start with ρ ∼ 1 + 1/2ρ1 + ρ2 and v ∼ v0 + 1/2v1 and finally obtainρ = 1 + o() and v = O(1/2).numerically.Observe that −χφ1(χ,∞) = A andφ1(χ,∞) = φ1(χ,−∞) +∫ ∞−∞φ1tˆ(χ, tˆ)dtˆ.At tˆ = −∞, v1 = 0, and φ1(χ,−∞) = 0. Then by using (5.39), we haveφ(χ,∞) =∫ ∞−∞[b2ρ1σ(χ, tˆ)∫ tˆ−∞v0(χ, sˆ)dsˆ+∫ χ−∞(−(ρ1v0)tˆ − (v20)σ −2σv20)dσ]dtˆ=∫ ∞−∞b2ρ1σ(χ, tˆ)∫ tˆ−∞v0(χ, sˆ)dsˆdtˆ︸ ︷︷ ︸T1−∫ ∞−∞∫ χ∞(ρ1v0)tˆdσdtˆ︸ ︷︷ ︸T2−∫ ∞−∞∫ χ∞(v20)σdσdtˆ︸ ︷︷ ︸T3−∫ ∞−∞∫ χ∞2σv20dσdtˆ︸ ︷︷ ︸T4.925.4. Formal Asymptotic AnalysisWe can simplify these terms and there is cancellation. Indeed,T1 =∫ ∞−∞−4bχe−tˆ2tˆ∫ tˆ−∞−2bχe−sˆ2dsˆdtˆ =∫ ∞−∞4√pib2χ2tˆe−tˆ2(1− erf(−tˆ))dtˆ=4√pib2χ2∫ ∞−∞tˆe−tˆ2erf(tˆ)dtˆ (by parts)=4√pib2χ2(−12e−tˆ2erf(tˆ)|∞−∞ +∫ ∞−∞1√pie−2tˆ2dtˆ) =2√2√pib2χ2T2 =∫ ∞−∞∫ χ∞(ρ1v0)tˆdσdtˆ =∫ χ∞∫ ∞−∞(ρ1v0)tˆdtˆdσ=∫ χ∞0dσ = 0T3 =∫ ∞−∞∫ χ∞(v20)σdσdtˆ =∫ ∞−∞v0(χ, tˆ)2dtˆ=∫ ∞−∞4b2χ2e−2tˆ2dtˆ =2√2√pib2χ2T4 =∫ ∞−∞∫ χ∞2σv20dσdtˆ = −∫ ∞−∞∫ ∞χ2σv20dσdtˆAdding all terms together gives us φ(χ, tˆ) =∫∞−∞∫∞χ2σv20dσdtˆ > 0. Thus,A < 0 and there is a remaining velocity field in the negative radial directionthat can compress the plasma. Obtaining a simple expression for the exactvalue of this integral is not possible, but by taking note of the characteristicshape of the integrand, we can approximate it with high precision. We935.4. Formal Asymptotic Analysis(b, χ) I∗approx I∗numerical(0.557, 1.19) 5.75 5.71(1.11, 1.19) 1.45 1.43(0.288, 1.19) 22.5 22.8(0.557, 2.38) 1.40 1.43(0.557, 0.595) 23.5 23.8Table 5.4: Verifying the numerical validity of the approximation given in(5.44) by doubling/halving the governing parameters and computing I∗ withthe approximation of v20 given in (5.44) and the numerically integrated value.remark thatv20 = [√pi2σ2(erf(σ − χb− tˆ)− erf(−χ− σb− tˆ))−1bσ(e−(σ−χb +tˆ)2+ e−(χ−σb +tˆ)2)]2≈1b2σ2(e−(σ−χb +tˆ)2+ e−(χ−σb +tˆ)2)2≈1b2σ2(e−2(σ−χb +tˆ)2+ e−2(χ−σb +tˆ)2)≈1b2σ2√pi√2(δ(tˆ−σ − χb) + δ(tˆ−χ− σb)) (5.44)where δ denotes the Dirac delta function [41]. We justify the first approxi-mation by noting the terms with erf are modulated by a 1/σ2 which decaysto 0 faster as σ → ∞ than terms with 1/σ and that for σ ≈ χ, the twoerror functions are being subtracted and have similar arguments. The sec-ond approximation comes from the fact that a Gaussian is dominated byits behaviour near its maximum and there is little overlap between the twoGaussians. The final approximation comes from approximating each Gaus-sian by its area modulating a delta function centred at its maximum.945.4. Formal Asymptotic AnalysisFrom the approximation in (5.44),φ1(χ,∞) = I∗ ≈∫ ∞−∞∫ ∞χ2σ1b2σ2√pi√2(δ(tˆ−σ − χb) + δ(tˆ−χ− σb))dσdtˆ=∫ ∞χ∫ ∞−∞√2pib2σ3(δ(tˆ−σ − χb) + δ(tˆ−χ− σb))dtˆdσ=∫ ∞χ2√2pib2σ3dσ =√2pib2χ2.Therefore A ≈√2pib2χ . For the remainder of this paper, we will takethis approximation as the value of I∗. Table 5.4 computes the integral∫∞−∞∫∞χ2σv20dσdtˆ using (5.44) and numerically, and the results show goodagreement. With this A, the long-time velocity profile is v ∼√v1 =√−√2pib2χσ2 (5.45). Phase III ResultsAfter the pulses have interacted with the inner wall, there remains a residualvelocity field which remains as a time-independent solution of the asymptoticequations. We find that ρ ∼ 1 + o() and v ∼ −√2pib2χσ2 1/2 + o(1/2).5.4.3 Phases IV and V5.4.3.1 An Energy ArgumentAt this point we now compute the remaining kinetic energy (all the particleswith a negative radial velocity), assuming there are no reflections at the rightboundary and that all the useful energy is in the vicinity of the inner wall.The kinetic energy density is E = 12ρv2, and to leading order this is E =12(1)(1/2v1)2 = 2v21. The total kinetic energy is 4pi∫ rRrLEr2dr, which in theinner variables at leading order is E = 5/22∫∞χ v1(σ)2σ2dσ = 4pi25/2b4χ2∫∞χdσσ2 .955.4. Formal Asymptotic AnalysisThe total possible compression energy isE =4pi25/2b4χ3. (5.46)The work done in compressing the plasma from rL = χ√ to rL = r∗ theminimum radius −∫PdV is −4pi∫ r∗χ√ P (rL)r2LdrL = −4pi∫ r∗χ√γr2LdrL givingus the total work doneW = 4piγ7/2(1r∗−−1/2χ). (5.47)Equating E = W in (5.46) and (5.47) implies 4pi25/2b4χ3 = 4piγ7/2( 1r∗ −−1/2χ ) so that r∗ = O() withr∗ ∼b4χ3γpi. (5.48)This sort of energy argument was originally used by Lord Rayleigh in consid-ering the compression of a bubble in an incompressible fluid [33]. We presenta formal asymptotic argument for this result in the next two subsections. Outer Region for Motion of Plasma BoundaryThe relevant region in space and time for this analysis is sketched in figure5.9.To proceed with formal asymptotics, we need to consider higher orders ofthe equations of mass and momentum conservation. By taking the equationsof (5.11) to two higher orders than presented in (5.36) with ρ ∼ 1 + 1/2ρ1 +ρ2 + 3/2ρ3 + 2ρ4 and v ∼ v0 + 1/2v1 + v2 + 3/2v3 we obtain the following965.4. Formal Asymptotic AnalysisFigure 5.9: Depiction of phase IV. This is the slow compression phase. Inthis phase, ρ ∼ 1 + O(2) and v = O(1/2). The peak velocity and pressureare growing in this phase.balances:ρ3,tˆ + v2,σ +2σv2 = −(ρ2v0 + ρ1v1)σ −2σ(ρ2v0 + ρ1v1) (5.49)v2,tˆ + b2ρ3,σ = −(ρ1v20 + 2v0v1)σ −2σ(ρ1v20 + 2v0v1)− (ρ2v0 + ρ1v1)tˆ (5.50)ρ4,tˆ + v3,σ +2σv3 = −(ρ3v0 + ρ2v1 + ρ1v2)σ −2σ(ρ3v0 + ρ2v1 + ρ1v2) (5.51)975.4. Formal Asymptotic Analysisv3,tˆ + b2ρ4,σ = −(ρ2v20 + 2ρ1v0v1 + 2v0v2 + v21)σ−2σ(ρ2v20 + 2ρ1v0v1 + 2v0v2 + v21)− (ρ3v0 + ρ2v1 + ρ1v2)tˆ (5.52)From our previous considerations of the behaviour of the solutions as tˆ→∞, we observe that (5.49) and (5.50) can also reach a steady-solution withρ3 = 0 and v2 = constant/σ2 since ρ1, ρ2, v0 all tend to 0 which eliminatesall forcing terms in the equations. Looking at (5.51) and (5.52) similarly, weactually have residual forcing terms that do not vanish as tˆ → ∞, namelythe terms with v21. After a long time, we infer that ρ ∼ 1 + 2ρ4 andv ∼ 1/2v1 + ... By scaling so that t = $, an O(1) time-scale we obtain anew PDE system to evolve, which serves as an outer region for a problemdescribing the motion of the plasma radius. Our balance is:v1,σ +2σv1 = 0, v1,$ + b2ρ4,σ + v21σ +2σv21 = 0. (5.53)Equation (5.53)1 gives an effective incompressibility to the lead-lithium andwe can proceed to solve (5.53) as done to derive the Rayleigh-Plesset equa-tion for bubble dynamics [5]. From (5.53)1 we have v1 = f($)/σ2 so thatupon substituting this into (5.53) we have f′σ2 + b2ρ4σ −2f2σ5 = 0 which canbe integrated from σ = ∞ to σ = σL($) the position of the plasma in-ner radius in the σ−coordinates to get −f′σL+ b2ρ4(σL) +f22σ4L= 0 where weused ρ4(∞, $) = 0. We have ρ4(σL, $) = 0 as well because if r = O(√),p = O(3/2) so that the density at the left boundary is 1 +O(5/2). We then985.4. Formal Asymptotic Analysisgetdfd$=f22σ3L. (5.54)Additionally, the velocity at the inner wall dσLd$ = v1(σL, %) (due to thescaling in this regime there is no extra factor of√ to deal with) sodσLd$=fσ2L. (5.55)Dividing (5.54) and (5.55) we find that dfdσL =f2σLwhich has solutionf = C√σL for a constant C. Thus we can treat v1 as a function of the innerwall position where v1(σL) = C/σ3/2L after using the solution in conjunctionwith (5.55). Givenv1(χ) = −√2pib2χ3(5.56)from equation (5.45), we have C = −√2pib2χ3/2and with respect to these outercoordinatesv1(σL) =−√2pib2χ3/2σ3/2L. (5.57)At this point, the system has not felt the effects of the plasma pressure andthe wall is accelerating inwards. Inner Region for Motion of Plasma BoundaryNumerous balances are possible for the system of (5.11) and it is possi-ble to arrive at the scaling we choose here considering all possible bal-ances and solutions and choose the one that can match to the outer re-gion, but we will consider a physical argument here to obtain the bal-ance. We seek a balance of terms where the pressure of the plasma re-sists the compression (i.e. the pressure gradient balances the momentum995.4. Formal Asymptotic AnalysisFigure 5.10: Depiction of phase V. This is the rapid compression phase wherethe plasma reaches its minimum radius. In this phase, ρ ∼ 1 +O(1/2) andv = O(−1/4).flux). We consider r = O(c). If r = O(c) then p = O(7/2−4c) andρ = 1+O(9/2−4c). Based on the growth predicted by (5.57) this would givea scaling of v = O( 1/2(c/1/2)3/2) = O(5/4−3c/2). The form of the scaling aswritten comes from the velocity growing like 1/r3/2L but being O(√) whenr = O(√). To balance the pressure/density ρ term in the momentum equa-tion with the v2 terms, we require 7/2 − 4c = 5/2 − 3c so that c = 1. Nowwe let r = 1z. Then v ∼ −1/4v−1, ρ ∼ 1 + 1/2ρ1, and to balance the scalesof velocity times time equals distance, we pick t = 5/4$ˆ where $ˆ measuresa time with respect to an arbitrary reference point. See figure 5.10. Theresulting balance is given byv−1,z +2zv−1 = 0 and v−1,$ˆ + b2ρ1,z + (v2−1)z +2zv−1 = 0. (5.58)Using the functional form of v−1 implied by (5.58)1, v−1 = g($ˆ)/z2L, in1005.4. Formal Asymptotic Analysis(5.58)2 and integrating as before, this time up to zL, the position of thewall in the inner coordinates gives us −g′zL+ γz4L+ g22z4L= 0 where we usedb2ρ1(zL, $ˆ) = γ/z4L. It follows then thatdgd$ˆ=γ + g2/2z3LanddzLd$ˆ=gz2Lso that after dividing the equations we get dgdzL =γ+g2/2zLgwhich we separateto get∫ gγ+g2/2dg =∫ dzLzLwith solution γ + g2/2 = DzL for a constant D.We can solve for g, obtaining g = −√2DzL − 2γ (the negative sign givingthe correct direction of the wall) which finally gives usv−1 =−√2DzL − 2γz2L. (5.59)Observe that when zL = z∗L = γ/D the velocity is 0 to leading order andthis will be our leading order estimate to the compression. To find D, weneed to perform matched asymptotics between the inner and outer regionswith respect to the functional dependence of the wall velocity upon the wallposition. Matching to Determine Minimum RadiusWe now put the outer variables of (5.57) in inner variables to match thebehaviours as zL →∞ and σL ↓ 0.−1/4(v−1)inner = 1/2(v1)outer−1/4−√2Dz3/2L= 1/2−√2pib2χ3/2(1/2zL)3/21015.4. Formal Asymptotic Analysisyielding −√2Dz3/2L= −√2pib2χ3/2z3/2Lso that D = pib4χ3 . In the z variable, we have aminimum radius of z∗L = γ/D =b4χ3γpi and we thus find r∗ = b4χ3γpi  consis-tent with (5.48) found with the energy argument. While the energy argu-ment was assumed true for incompressible fluids, it seems in the asymptoticlimit compressibility does not influence the compression to leading order.Note that by having quarter powers of , it may become difficult to distin-guish different asymptotic terms without taking  to be extremely small.When v = O(−1/4), this is very close to v = O(−1/2), the sound speedmagnitude and the wave-like behaviour of the equations may be close tobreaking down. Indeed, when  is too large, the local velocity magnitudecan even exceed the sound speed during the compression phase. Dimensional Expression for Plasma CompressionWe can finally go back to the dimensional parameters in our problem. Usingtables 5.2 and 5.3, we find the minimum radius isR∗L ≈C4sPplasma,0R7inner,0%30piP 4maxR4outer,0T20(5.60)1025.5. Comparison with NumericswhereT0Router,0P 1/2max%1/20 1 (5.61)CsT1/40 P1/4maxR1/2outer,0%1/40= O(1) (5.62)Pplasma,0R4inner,0%7/40P 11/4max R1/2outer,0T7/40= O(1) (5.63)Rinner,0%1/40R1/2outer,0T1/20= O(1). (5.64)Equation (5.60) gives the approximate expression for the minimum ra-dius, which requires that  is small (5.61) and that b, γ, and χ all be O(1)as given in equations (5.62), (5.63), and (5.64).5.5 Comparison with NumericsIn chapter 4, a finite volume numerical framework was developed to solve thenonlinear conservation laws with the moving boundaries. Using this numer-ical framework, we run simulations on the asymptotic problem formulatedin this paper. Our validation is done in multiple stages: we examine theprofiles of the pulses during their formation both numerically and asymp-totically by comparing the profiles of the resultant variables; then, we studythe growths of the amplitudes of the pulses as they move radially inward aspredicted by the asymptotics and numerics; we next proceed to study theexistence of the residual velocity field with the numerics and the growth ofthe inner wall velocity with inner wall position; we observe that the velocityat the inner wall is qualitatively consistent with the numerical results; and1035.5. Comparison with Numericsfinally, we compute the minimal radius in the asymptotic and numericalworkings for various values of b, γ, χ, and .We remark that the asymptotic problem we have considered is primar-ily concerned with a single direction of information propagation: the pulsetravels towards the plasma wall, it is reflected, and then everything elsegoverning the plasma wall takes place locally near the wall and there is noinfluence from leftgoing waves. We need to exercise caution with this nu-merically by removing reflections that occur at the outer wall. This allowsfor better agreement between the asymptotics and numerics but we maketwo remarks. We first note that if these reflections are not removed then areflection occurs at the outer wall and disturbances return to hit the plasmaagain which drastically reduces the compression. Such effects are noticeablefor  . 0.0025. The second point to make is that even in attempting to re-move the reflections at the outer wall, there is likely still a small numericalerror that remains.5.5.1 Pulse FormationWe consider two parameter regimes here. One is the for the parameter setof interest, and the other is for a much smaller  with other dimensionlessparameters set to unity. We fix an end time t = /2 and plot the profiles of vand ρ as computed numerically, and with two different levels of asymptoticaccuracy. From the plots, we are able to verify that the successive terms ofthe asymptotic expansion yield higher accuracies, and that the asymptoticsand numerics are in good agreement. See figure 5.5.1045.5. Comparison with Numerics0 0.2 0.4 0.6 0.8 10246810Maximum PositionMaximum Pressure Value  Numerical PeakAsymptotic Peak0 0.2 0.4 0.6 0.8 1−0.4−0.3−0.2−0.10Minimum PositionMinimum Velocity Value  Numerical PeakAsymptotic PeakFigure 5.11: Numerical and asymptotic descriptions of the growth in peakvalues for the minimum value of velocity and peak value of pressure. Basedon a best fit, the numerically observed scalings are 1/r0.99 and 1/r0.90 forthe left and right plot respectively. Parameters: b = χ = γ = 1,  = FocusingHere, we consider the asymptotically predicted growth rates for the peakamplitudes of the velocity and pressure. From the asymptotic predictions,both of these amplitudes should grow inversely with the position of thesepeaks. This validation is rather delicate: these growths should be upheld inthe limit  ↓ 0; however, there are difficulties in getting the numerics to givehighly accurate results on regions where r  1. What we choose to do ispick a very small value of , 0.001, and consider the growth of the amplitudeson the region r > 0.1 ≈ O(√). We plot the predictions as given by a perfectgrowth of 1/r based on the initial peak pressure and minimum velocity, andthe numerical results for the peak amplitude vs position for this  value.The plots are given in figure 5.11 and we observe strong consistency.1055.5. Comparison with Numerics0 0.1 0.2 0.3−0.3−0.25−0.2−0.15−0.1−0.050Position rVelocity v L  NumericalAsymptoticFigure 5.12: A plot of the numerical and asymptotic velocity profile aftermuch of the pulses have reflected off of the inner wall. The two are nearlyindistinguishable. Parameters: b = χ = γ = 1,  = 0.0025, t = 0.0784.5.5.3 Residual Velocity FieldOur prediction of a steady velocity field resulting from the disturbancesinteracting with the plasma was based primarily on intuition, as we did notformally find the solution for v1 and ρ2 for general tˆ; we simply showed thatthe potential approaches a constant value along the boundary. Numerically,however, we can validate the scaling. After a long time with respect to the tˆtime-scale so that v0 has had its full effect, we plot the profile of the velocityversus the radial position. As the inner wall moves inward, in its wake thereis a velocity profile that scales like 1/r2L. See figure Outer Region Describing Inner Wall VelocityThis validation is rather delicate as there are a number of sources of error.Firstly, the inner wall position at which the boundary condition of (5.56) isbetter approximated by σ = χ −√2√pibχ as per (5.43) but to leading orderwe have taken it as σ = χ. In the case considered in figure 5.13, this already1065.5. Comparison with Numerics0.01 0.015 0.02 0.025 0.03 0.035−1−0.8−0.6−0.4−0.2Inner Wall Position rLInner Wall Velocity v L  NumericalAsymptoticFigure 5.13: A plot of the numerical and asymptotic velocity profile forhow the inner wall velocity implicitly depends on the inner wall position.Parameters: b = χ = γ = 1,  = 0.0025, t ranging from 0.0784 to 0.11.amounts to an error of 17% in estimating the boundary condition whichwill lead to an error in the constants obtained. We additionally know thatv ∼√v1 +v2 so that there is an O() error term in the radial position thatis accrued over a time scale of O(1). Throwing these effects together makesit very hard to get a clean fit between numerical and asymptotic results. Infinding the slope of the log vL vs log rL plot, a least squares fit shows thenumerical scaling is vL ∝ 1/r1.49L which is completely consistent with theasymptotic scaling.5.5.5 Qualitative Agreement of Inner Wall VelocityThe asymptotics predict a series of phenomena at the inner wall: firstly, thewall is stationary until the pressure pulse reaches it whereupon it takes on aGaussian shape then decreases in speed; secondly, a smaller asymptotic termdescribing the wall velocity remains for some time (as the velocity acquiresthe residual profile); thirdly, the wall rapidly speeds up; and finally, the wallis stopped abruptly on a small spatial scale when the plasma pressure finally1075.5. Comparison with NumericsFigure 5.14: A numerical validation that the asymptotics have qualitativelycaptured the behaviour of the inner wall. Parameters: b = χ = γ = 1, = 0.0025. The different qualitative regimes are labelled in the diagram.The r∗-value that we seek occurs at the bottom spike of the curve where themaximum compression label appears.takes over. We observe all of these phenomena in figure Predictions for Plasma CompressionIn this section, we verify that the fundamental asymptotic predictions areconsistent with the numerics. Table 5.5 presents the results. We remarkthat verifying the asymptotic limit is not trivial. As  ↓ 0, the numerics, formodest discretizations, lose accuracy due to the very small radial positionsunder consideration. As a result, we cannot not take  too small. We alsohave to ensure that quantities such as b4χ3γ/pi, etc. remained roughly O(1)for these values of . Picking  too large also leads to problems: for one,the higher order terms may dominate over the desired r∗ = O() behaviour.Another issue is that if  is not small enough, the apparently negligible dis-placement of the inner wall given by (5.43) could be larger than the initialinner radius of χ√. In table 5.5, we compute the minimum radius asymptot-1085.5. Comparison with Numerics r∗asy r∗num Error0.02 0.02000 0.01155 0.008450.01 0.01000 0.00600 0.004000.005 0.00500 0.00444 0.000560.0025 0.00250 0.00242 0.00008Table 5.5: Asymptotic and numerical predictions of minimum radius ofplasma for different values of  with b = 1.05, χ = 0.937 and γ = pi. Basedon a least squares linear regression in a log-log plot of the difference versus, the convergence rate appears to be O(2.3) = o(). The convergence rate isexcellent likely either due to small coefficients in the o() asymptotic seriesterms or fortuitous cancellations of higher order asymptotic error terms.ically and numerically. We verify that the difference in the minimum radiusbetween the asymptotic and exact (numerical) predictions is o() when b, χ,and γ are fixed. We choose b = 1.05, χ = 0.937, γ = pi so that the minimumradius should be 0.99994. Comments on the Model Parameter RegimeWith respect to the parameters of table 5.3, the asymptotic prediction of theminimum radius is r∗ = 0.00229 and the numerical results are larger with avalue of 0.00649. The two predictions are roughly on the same scale, but arenot quantitatively consistent. We speculate this is the result of the fact thatb4χ3γ/pi = 0.18, which is approximately O(√) and not O(1) with respectto  = 0.0126. Somewhere in the sequence of asymptotic regions, the asymp-totic series lost its ordering and this is likely the cause of the discrepancy. Wewould naturally expect these results to be more accurate if  were smaller orb4χ3γ/pi were nearer unity (for example if we use b = 0.8 and  = 0.01 insteadof b = 0.557 and  = 0.0126 so that b4χ3γ = 0.77, the asymptotic predictionis 0.0077 and the numerical result is 0.0088). Despite the discrepancy, thereis qualitative consistency between the numerical and asymptotic modelling.1095.5. Comparison with NumericsIn previous work [27] where the model was very similar, the data show thatperturbations to the inner radius had the greatest impact upon the com-pression. This was followed by a tie between maximum external pressureand maximum outer radius, and then pulse time scale, sound speed, initialplasma pressure. This is reasonably consistent with (5.60), where based onthe respective powers of the parameters, the relative sensitivity of the mini-mal radius with respect to the parameters (i.e., the percentage by which theminimal radius would change for a percent change in the parameter), frommost to least sensitive follows the ordering: initial plasma radius; soundspeed, initial outer radius, and maximum impulse pressure (three-way tie),pulse time scale, and initial plasma pressure. The proper positive/negativecorrelations are also consistent i.e. when the asymptotics show an increasein a parameter decreases the radius, so do the numerics.With the results of the asymptotics, we can make a few application-relevant statements. Based on (5.14)2, (5.24) and (5.46), we note that onlya small fraction of the energy input actually goes towards compressing theplasma. Indeed if we compute the energy input by∫ ∞−∞−4pirR(t)2pR(t)︸ ︷︷ ︸forcevR(t)dt︸ ︷︷ ︸displacement∼4pi3/2b∫ ∞−∞e−2τ2dτ =√8pi3b3/2,the leading order energy is O(3/2) but only O(5/2) ultimately goes towardscompression. A lot of energy is lost in reflection. In our modelling, we ne-glected the gas pressure of the plasma, which would also work to oppose thecompression of the plasma; the degree the plasma is compressed is likely lessthan what we predict asymptotically, i.e., in a model with the gas pressure,we would expect r∗ to be larger than our value given in (5.48) at leading or-1105.6. Conclusions and Future Workder. Also, based on the previous work mentioned in the previous paragraph,we can see that a smaller minimum radius does not necessary yield morepromising fusion conditions (based on the Lawson triple product criterion[21]); the time over which the plasma is compressed is another vital element.While our work here predicts the minimum radius, this is only one of manycomplex components required for the success of magnetized target fusion.5.6 Conclusions and Future WorkIn this chapter, an analytic result for the minimal radius of a plasma has beenobtained in the limit of a very fast impulse time scale. Although the modelis highly simplified, we can describe qualitatively the effects of key designparameters in the magnetized target fusion model in question. Equation(5.60) is qualitatively accurate, given the stipulations outlined in equations(5.61) - (5.64), and it provides a good ballpark quantitative estimate forthe minimum radius. The key parameters that are within control are likelyPmax, Router,0, and T0. In this case, within the limitations of the physicalmodel, the plasma is compressed to a smaller and smaller radius as thepiston pressure, initial outer radius, and pulse time scale all increase. If themedium through which the pressure pulses travel could be modified thendecreasing either its sound speed or density (or both) would be ideal whilesustaining the piston pressure. Also, if the initial plasma pressure or sizecould be decreased, a greater compression can take place.Although this problem contained nontrivial obstacles, including nonlin-ear conservation laws and moving boundaries, through appropriate scalingarguments and suitable solution techniques and estimates, an analytic re-sult that agrees well with numerical simulations has been obtained in the1115.6. Conclusions and Future Worklimit where  ↓ 0. Building upon this asymptotic framework, many new andimportant studies could potentially be done to enhance the physical accu-racy and predictive power of the modelling including: adding more detailedplasma physics, including adding in the gas pressure and allowing for energylosses; considering the plasma and lead-lithium interaction more carefully;incorporating a more sophisticated equation of state for the lead-lithium; oreven studying the effects of localized angular perturbations in the systemgiven that perfect spherical symmetry is not possible in the real device.112Part IVFurther Details andConclusions113Chapter 6Additional InsightsConcerning the Numericaland Asymptotic Modelling ofthe Fusion Reactor6.1 Fusion Reactor Numerical DetailsThis section provides further details on the numerical simulations studiedin chapter 4.Our implementation follows from the discussion found in 3.3.1, but weprovide the background to the coordinate transformation and setup. Adetailed pseudo-code is provided in the appendix A.6.1.1 Transformed Coordinate SystemTo avoid dealing with moving boundaries within each time step, we chose totransform the coordinate system to a fixed computational domain arriving atequations (4.12) and (4.13). The reader may wish to refer back to equations(4.1) through (4.11) for the setting. We begin with the spherically symmetric1146.1. Fusion Reactor Numerical Detailscompressible Euler equations,ρt + (ρv)r +2rρv = 0(ρv)t + (P + ρv2)r +2rρv2 = 0where rL(t) ≤ r ≤ rR(t) with r′L(t) = v(rL(t), t) and r′R(t) = v(rR(t), t). Weconsider a fixed computational domain in space y ∈ [0, 1] where y = r−rLrR−rLand a time parameter τ = t. If we denote ∆(t) = rR(t) − rL(t) and Γ(t) =vR(t)− vL(t) then:∂r =∂y∂r∂y=1∆∂y∂t =∂τ∂t∂τ +∂y∂t∂y= ∂τ +−vL∆− Γ(r − rL)∆2= ∂τ +−vL∆− Γ∆y∆2= ∂τ +−vL − Γy∆With this, we are now able to rewrite the equations for mass and momentumconservation as:ρτ +−vL − Γy∆ρy +1∆(ρv)y +2rL + ∆yρv = 0(ρv)τ +−vL − Γy∆(ρv + P )y +1∆(ρv2)y +2rL + ∆yρv2 = 0.The numerical methods work best when spatial derivative terms are writ-1156.1. Fusion Reactor Numerical Detailsten in a conservative form ut + (f(u))x = 0 as we observe in section 3.3.3.Observe we can rewrite the system in a more conservative form as:ρτ + (−vL − Γy∆ρ+1∆ρv)y +Γρ∆+2rL + ∆yρv = 0(ρv)τ + (−vL − Γy∆(ρv) +1∆P +1∆ρv2)y +Γρv∆+2rL + ∆yρv2 = 0.This is the system given by (4.12) and (4.13).6.1.2 FormulationIn vector form our equations have a flux function−vL−Γy∆ ρ+1∆ρv−vL−Γy∆ (ρv) +1∆P +1∆ρv2with linearized eigenvaluesλ± =v − (vL + Γy)±√P ′(ρ)∆.Our vector that we called u in the discussion of the finite volume method inchapter 3 is taken to be (ρ, ρv)T .Even with this conservative-like form, there are still terms that arey−dependent (and τ−dependent) that appear outside of the fluxes. Aswe cannot avoid these terms, we anticipate our method being reduced fromsecond to first-order [23].1166.2. Asymptotic Commentary6.2 Asymptotic Commentary6.2.1 Different Model than Numerical SimulationsThe asymptotic model in chapter 5 entailed a slightly different formulationthan the numerical modelling done in chapter 4 notably in the fact thatthe asymptotic model worked with a linearized equation of state for lead-lithium, that the plasma only involved a magnetic pressure, and that thepiston impulse was continuous and smooth. Our choice of a more smoothimpulse function was motivated by an attempt to avoid, if possible, nonlinearsystems with shocks and possibly delta-functions (the spatial derivative ofstep functions) appearing in the analysis. We comment, however, that giventhe nature of a Gaussian impulse and its fast time scale, both the continuousand discontinuous idealizations of how the pistons impact the outer wall, dodescribe the same system, and, as we observe in the qualitative analysis, theasymptotics, even with this more smooth impulse function, do predict therelative influence of the different design specifications.6.2.2 Suitability of Linear Lead-Lithium Equation of StateFrom the numerical simulations for the nuclear fusion reactor, we observethat the peak total plasma pressure (and peak pressure within the entiresystem) is approximately 6.7 GPa from table 4.3. With this maximumpressure, we are able to investigate where this falls on the quadratic equationof state curve, and determine how well the linear equation of state used inthe asymptotics describes the behaviour. The red line in figure correspondsto the peak pressure of 6.7 GPa and we see that at this point, the quadraticdoes deviate from the straight line, but the relative error in the densitiescorresponding to this pressure is only around 5%. The simple assumption of1176.2. Asymptotic CommentaryFigure 6.1: Plot of quadratic and linear equation of state models with peakpressure as computed numerically.a linear equation of state appears to be acceptable.In the following subsection, we illustrate how the plasma gas pressurecould be incorporated into the asymptotic analysis.6.2.3 Inclusion of Gas PressureThe asymptotic results illustrated above are valid for a plasma pressuredescribed solely by its magnetic pressure. We can, however, add the gaspressure at the expense of not being able to solve the equations of motionin regions IV and V in 5.4.3 and relying upon the energy argument.If we nondimensionalize equation (4.6) by the parameters listed in table5.2 then we haveP¯PL =κ1R¯4r4+κ2R¯5r51186.2. Asymptotic Commentaryor thatp =γ7/2r4+µ9/2r5where µ = κ2R¯5P¯−9/2. Numerically, µ ≈ 3.380.Even with the inclusion of the gas pressure, the asymptotic order of theinitial plasma pressure remains O(3/2) so that all the analysis is the sameup to the energy argument. Then, with this new pressure, we can computethe work done in compressing the plasma from r = χ√ to r = r∗:W = −∫ r∗χ√4pir2(γ7/2r4+µ9/2r5)dr= 4pi(γ5/2r+µ9/22r2)|r∗χ√= 4pi(−γ2χ+γ5/2r∗−µ22χ2+µ5/22r∗2) (6.1)In equating the remaining kinetic energy equation (5.46) with the workdone in equation (6.1), we find that to leading order, r∗ = s wherepib4χ3=γs+µ2s2.This is a quadratic in 1/s so that the solution is thats = (−γ +√γ2 + 2µpib4χ3µ)−1 =µ−γ +√γ2 + 2µpib4χ3(6.2)where we have taken the + sign in using the quadratic formula because s > 0is physically relevant.At first glance equation (6.2) may not appear consistent with (5.48) in1196.2. Asymptotic Commentarysetting µ = 0, however for µ 1, we can write√γ2 + 2µpib4χ3 as γ√1 + 2µpib4χ3γ2 =γ(1 + µpib4χ3γ2 + o(µ)) so that as µ ↓ 0, equation (6.2) readss =µ−γ + γ + µpib4χ3γ + o(µ)=b4χ3γpi+ o(µ).With the gas pressure included, we estimate the minimum radius to ber∗ =µ−γ +√γ2 + 2µpib4χ3 (6.3)We justify that this claim is still valid in the following section. Numerically,with the value 3.380 and the values of b, χ, γ, and  as in chapter 5, thisevaluates to r∗ ≈ 0.005, which is larger than in chapter 5 where the gaspressure was neglected. Our neglecting the gas pressure in chapter 5 wasnot done with the intent of stating the that it is negligible in the overalldevice performance (indeed, we see that it does play a leading-order role),but rather that for much of the compression, it does not appear, and withthe asymptotic work we sought a basic qualitative description of the system,which requires an analytic description of the compression.Further asymptotic work could possibly be done by including correctionsto the lead-lithium equation of state; however, the asymptotic work hasalready reached the edges of validity.6.2.4 Justification of Equation (6.3)The energy argument attributed to Lord Rayleigh holds for incompressiblefluids, but our system of interest involves the compressible Euler equations.However, if the leading order asymptotic equations maintain the incom-pressible limit as in chapter 5, we would expect (6.3) to hold asymptotic1206.2. Asymptotic Commentaryconsistency.Even with the modified pressure for the plasma, at a scale r = O()and a time t = O(5/4), we still have a balance of terms with ρ ∼ 1 + 1/2ρˆand v ∼ −1/4vˆ. As a result, the leading order equations that describe thecollapse in this inner region, although unsolvable analytically, still imposethatvˆσ +2σvˆ = 0so that the leading order behaviour is still that of an incompressible fluid.This assumes that during the compression, the peak density, after reflection,occurs at the plasma-lead-lithium interface, which is a reasonable assump-tion. If this were not the case, we could imagine that the leading orderdensity perturbation could be larger than O(1/2) somewhere away fromthe boundary and the asymptotics would break down. However, after thereflection, the problem seemingly amounts to a signalling problem with dis-turbance source located at the inner wall, and we would expect the peakvalues to occur there.For a compressible fluid, the results could be quite different becauseto leading order, an asymptotically non-negligible energy could be storedwithin the lead-lithium during compression. Observe that the kinetic energydensity is K = 12ρv2, which, to leading order with ρ ∼ 1 + 1/2ρˆ and v ∼−1/4vˆ implies K = O(−1/2).To compute the potential energy density in our setting, we will considerthe work done in compressing a fluid occupying a volume Vi at constantdensity ρi and pressure p(ρi) down to a volume Vf with constant density ρfand pressure p(ρf ). By conservation of mass, at any intermediate volumeV where the density is ρ, we have that ρiVi = ρV so that V = ρiVi/ρ. The1216.2. Asymptotic Commentarywork done isW = −∫ VfVip(ρ(V ))dV=∫ ρfρiρiViρ2p(ρ)dρ= ρiVi∫ ρfρip(ρ)ρ2dρ.The second equality arose from replacing dV by dVdρ dρ and switching ourvariable of integration from V to ρ.We can use any reference density, and we choose ρi = 1. Then, by notingthat ρiVi = ρfVf we have the work done isW = ρfVf∫ ρf1p(ρ)ρ2dρ.Given this amount of potential energy in a volume Vf , we can find thepotential energy density by dividing by Vf , and if our fluid is at a denistyρf = ρ then the potential energy density becomesP = ρ∫ ρ1p(s)s2ds.1226.2. Asymptotic CommentaryWith p(s) = b2 (s− 1), ρ = 1 + 1/2ρˆ, we findP =1 + 1/2ρˆ∫ 1+1/2ρˆ1b2(s− 1)s2ds=1 + 1/2ρˆb2∫ 1+1/2ρˆ1(1s−1s2)ds=1 + 1/2ρˆb2(log s+1s)|1+1/2ρˆ1=1 + 1/2ρˆb2(log(1 + 1/2ρˆ)− log 1 +11 + 1/2ρˆ− 1)=1 + 1/2ρˆb2(1/2ρˆ−2ρˆ2 + 1− 1/2ρˆ+ ρˆ2 − 1 +O(3/2))=b22ρˆ2 +O(1/2).We see from this that the leading order potential energy contribution isa factor of 1/2 smaller than the leading order kinetic energy contributionand compressibility effects are not asymptotically significant.6.2.5 Numerical ValidationHere we validate equation (6.3) through numerical simulations. The dataof table 6.1 validate our prediction with one particular realization of theconstants b, χ, γ, and µ. As the error appears to be o(), these results supportequation (6.3).1236.2. Asymptotic Commentary r∗asy r∗num Error0.02 0.02883 0.01637 0.012460.01 0.01442 0.00848 0.005940.005 0.00721 0.00603 0.00118Table 6.1: Asymptotic and numerical predictions of minimum radius ofplasma for different values of  with b = 1.05, χ = 0.937, γ = pi, and µ = 4.The asymptotic minimum radius should be 1.4415. The apparent conver-gence rate is O(1.7).124Chapter 7Summary and Future WorkThis chapter provides a summary and considers the future work for theinvestigation into magnetized target fusion.7.1 SummaryIn this thesis we developed a simple model for a magnetized target fusionreactor, and analyzed it from a numerical and asymptotic perspective. Bothmethodologies validated each other in their model predictions. Our numeri-cal work shows us that nuclear fusion energy may be within reach of currentengineering efforts when considered in the light of the sensitivity analysis.The model has a high degree of sensitivity to its parameters, and while thebaseline parameters do not suggest a net energy production, small changesin these parameters could very well push the design into a regime with a netenergy yield.The asymptotic work has allowed for a qualitative description of thereactor operation, illustrating how the pressure pulse works in compressingthe plasma, how much of the input energy is reflected before compressingthe plasma, and how various design parameters should affect the operation,which are consistent with the numerics. We find that an increase in thepiston impact pressure and its time scale, along with an increase in the radiusof the lead-lithium sphere would tend to improve the device performance.1257.2. Future WorkAlso, although counterintuitive, the notion arises that having a transmittingmaterial with a lower sound speed may make the device more effective if theother parameters could remain intact.This work has brought up some interesting ideas, but they should beinterpreted within the realm of our assumptions. Some notable assump-tions were: assuming spherical symmetry, considering an adiabatic plasmawith a highly simplified equation of state, neglecting viscous terms in thelead-lithium, not allowing mixing of the plasma and lead-lithium, and ex-trapolating experimental data for lead to predict the equation of state forlead-lithium down to negative (cavitation) pressures. This work is really justa first step in a much bigger exploration and development of fusion energytechnology.7.2 Future WorkMuch of the work to extend this investigation should focus upon buildinga more accurate model. This could include: a higher dimensional numer-ical framework, more sophisticated equations of state for the plasma andlead-lithium, including energy balances and losses in the system, consider-ing angular instabilities during the implosion, and having a more preciseinteraction between the plasma and its surrounding fluid.126Bibliography[1] Adiabatic invariants in astrophysical plasma. In Plasma Astrophysics,Part I, volume 340 of Astrophysics and Space Science Library, pages103–114. Springer New York, 2006.[2] Avoiding fossil fuel costs with wind energy, 2014.[3] L. Artsimovich. Controlled Thermonuclear Reactions. Hazell Watson& Viney Ltd., 1964.[4] F Bourennani, S Rahnamayan, and GF Naterer. Optimal design meth-ods for hybrid renewable energy systems. International Journal ofGreen Energy, 12(2):148–159, 2015.[5] C. Brennen. Cavitation and Bubble Dynamics. Oxford University Press,Open Access, 1995.[6] C. A. Coulson and A. Jeffrey. Waves: A mathematical approach to thecommon types of wave motion. Longman Inc., New York, 1977.[7] M. J. Edwards, P. K. Patel, J. D. Lindl, L. J. Atherton, S. H. Glenzer,S. W. Haan, J. D. Kilkenny, O. L. Landen, E. I. Moses, A. Nikroo,R. Petrasso, T. C. Sangster, P. T. Springer, S. Batha, R. Benedetti,L. Bernstein, R. Betti, D. L. Bleuel, T. R. Boehly, D. K. Bradley, J. A.Caggiano, D. A. Callahan, P. M. Celliers, C. J. Cerjan, K. C. Chen,127BibliographyD. S. Clark, G. W. Collins, E. L. Dewald, L. Divol, S. Dixit, T. Doepp-ner, D. H. Edgell, J. E. Fair, M. Farrell, R. J. Fortner, J. Frenje, M. G.Gatu Johnson, E. Giraldez, V. Yu. Glebov, G. Grim, B. A. Hammel,A. V. Hamza, D. R. Harding, S. P. Hatchett, N. Hein, H. W. Herrmann,D. Hicks, D. E. Hinkel, M. Hoppe, W. W. Hsing, N. Izumi, B. Jacoby,O. S. Jones, D. Kalantar, R. Kauffman, J. L. Kline, J. P. Knauer, J. A.Koch, B. J. Kozioziemski, G. Kyrala, K. N. LaFortune, S. Le Pape,R. J. Leeper, R. Lerche, T. Ma, B. J. MacGowan, A. J. MacKinnon,A. Macphee, E. R. Mapoles, M. M. Marinak, M. Mauldin, P. W. McK-enty, M. Meezan, P. A. Michel, J. Milovich, J. D. Moody, M. Moran,D. H. Munro, C. L. Olson, K. Opachich, A. E. Pak, T. Parham, H.-S.Park, J. E. Ralph, S. P. Regan, B. Remington, H. Rinderknecht, H. F.Robey, M. Rosen, S. Ross, J. D. Salmonson, J. Sater, D. H. Schneider,F. H. Sguin, S. M. Sepke, D. A. Shaughnessy, V. A. Smalyuk, B. K.Spears, C. Stoeckl, W. Stoeffl, L. Suter, C. A. Thomas, R. Tommasini,R. P. Town, S. V. Weber, P. J. Wegner, K. Widman, M. Wilke, D. C.Wilson, C. B. Yeamans, and A. Zylstra. Progress towards ignitionon the national ignition facilitya). Physics of Plasmas (1994-present),20(7):–, 2013.[8] Tom E Faber. Fluid dynamics for physicists. Cambridge UniversityPress, 1995.[9] Nikolai N Gorelenkov, SD Pinches, and K Toi. Energetic particlephysics in fusion research in preparation for burning plasma experi-ments. Nuclear Fusion, 54(12):125001, 2014.[10] T Gray, MR Brown, CD Cothran, G Marklin, and MJ Schaffer. Stable128Bibliographyspheromak formation by merging in an oblate flux conserver. Physicsof Plasmas (1994-present), 17(3):032510, 2010.[11] R. D. Griffin. Alternative energy, 1992.[12] David Jeffrey Griffiths. Introduction to electrodynamics, volume 3.Prentice Hall Upper Saddle River, NJ, 1999.[13] Eric Herbert, Se´bastien Balibar, and Fre´de´ric Caupin. Cavitation pres-sure in water. Physical Review E, 74(4):041603, 2006.[14] E John Hinch. Perturbation methods. Cambridge university press, 1991.[15] JD Jackson. Catalysis of nuclear reactions between hydrogen isotopesby µ-mesons. Physical Review, 106(2):330, 1957.[16] Minsung Kim, Gilbong Lee, Young-Jin Baik, and Ho-Sang Ra. Per-formance evaluation of geothermal heat pump with direct expansiontype vertical ground heat exchanger. Heat Transfer Engineering, (just-accepted):00–00, 2014.[17] Charles Kittel. Elementary statistical physics. Courier Corporation,2004.[18] Michel Laberge. An acoustically driven magnetized target fusion reac-tor. Journal of Fusion Energy, 27(1-2):65–68, 2008.[19] Michel Laberge. Experimental results for an acoustic driver for mtf.Journal of fusion energy, 28(2):179–182, 2009.[20] LD Landau and EM Lifshitz. Fluid mechanics, vol. 6. Course of The-oretical Physics, pages 227–229, 1987.129Bibliography[21] John D Lawson. Some criteria for a power producing thermonuclearreactor. Proceedings of the Physical Society. Section B, 70(1):6, 1957.[22] L Gary Leal. Advanced transport phenomena: fluid mechanics andconvective transport processes. Cambridge University Press, 2007.[23] Randall LeVeque. Numerical methods for conservation laws, volume132. Springer, 1992.[24] Yonghua Li and Xiaoping Wang. Theoretical calculation and economicassessment of the integrated solar combined cycle power plants. Inter-national Journal of Green Energy, 12(3):192–197, 2015.[25] Michael Lindstrom. Asymptotic analysis of a magnetized target fusionreactor. SIAM Journal on Applied Mathematics, 2015 (under revision).[26] Michael Lindstrom, Sandra Barsky, and Brian Wetton. Investigationinto fusion feasibility of a magnetized target fusion reactor: A prelim-inary numerical framework. Journal of Fusion Energy, 34(1):76–83,2015.[27] Michael Lindstrom, Brian Wetton, and Rob Kiefl. Mathematical mod-elling of the effect of surface roughness on magnetic field profiles in typeii superconductors. Journal of Engineering Mathematics, 85(1):149–177, 2014.[28] DK Mansfield, KW Hill, JD Strachan, MG Bell, SD Scott, R Budny,ES Marmar, JA Snipes, JL Terry, S Batha, et al. Enhancement of toka-mak fusion test reactor performance by lithium conditioning. Physicsof Plasmas (1994-present), 3(5):1892–1897, 1996.130Bibliography[29] S.H. Mohr, J. Wang, G. Ellem, J. Ward, and D. Giurco. Projection ofworld fossil fuels by country. Fuel, 141(0):120 – 135, 2015.[30] John Ockendon, Sam Howison, Andrew Lacey, and Alexander Movchan.Applied partial differential equations. OUP Catalogue, 2003.[31] Abdeen Mustafa Omer. Energy use and environmental impacts: A gen-eral review. Journal of Renewable and Sustainable Energy, 1(5):053101,2009.[32] Kenneth JH Phillips. Guide to the Sun. Cambridge University Press,1995.[33] Lord Rayleigh. Viii. on the pressure developed in a liquid during the col-lapse of a spherical cavity. The London, Edinburgh, and Dublin Philo-sophical Magazine and Journal of Science, 34(200):94–98, 1917.[34] SD Rothman, JP Davis, J Maw, CM Robinson, K Parker, and J Palmer.Measurement of the principal isentropes of lead and lead–antimony al-loy to˜ 400 kbar by quasi-isentropic compression. Journal of PhysicsD: Applied Physics, 38(5):733, 2005.[35] Bimalendu N Roy. Fundamentals of classical and statistical thermody-namics. John Wiley & Sons, 2002.[36] Walter A Strauss. Partial differential equations. an introduction. NewYork, 1992.[37] Victoria Suponitsky, Aaron Froese, and Sandra Barsky. Richtmyer–meshkov instability of a liquid–gas interface driven by a cylindricalimploding pressure wave. Computers & Fluids, 89:1–19, 2014.131Bibliography[38] Sirui Tan and Chi-Wang Shu. Inverse lax-wendroff procedure for nu-merical boundary conditions of conservation laws. Journal of Compu-tational Physics, 229(21):8144–8166, 2010.[39] Gerald Beresford Whitham. Linear and nonlinear waves, volume 42.John Wiley & Sons, 2011.[40] S Woodruff, AID Macnab, and N Mattor. Adiabatic compression ofa doublet field reversed configuration (frc). Journal of Fusion Energy,27(1-2):128–133, 2008.[41] Erich Zauderer. Partial differential equations of applied mathematics,volume 71. John Wiley & Sons, 2011.132Appendix ADetailed Pseudo-Code ofFusion Finite VolumeSchemeThe pseudo-code here provides the details of the program mentioned inchapter 4.After the initializations, our system consists of 2 variables ρ and ρv withP = P (ρ) in the coordinate system (y, τ). The variables are defined on themesh yi = ih, h = 1/N , i ∈ {0, 1, ..., N}. The outline below defines a singleforward step in time of size k. We will denote u = (ρ, ρv)T as our unknown.1. Having found the largest eigenvalue of the linearized system matricesAi about each ui, we determine a suitable time step.2. We extrapolate u and y to indices i = −2,−1 and i = N + 1, N + 2.Extensions of v are defined by ρv/ρ. We are also able to extend thevalues of the sound speed by computing√dP/dρ.3. We compute the difference vector of u, defined at indicesi = −3/2,−1/2, 1/2, ..., N − 1/2, N + 1/2, N + 3/2.4. We interpolate values for u at half-integer indices by averaging, and133Appendix A. Detailed Pseudo-Code of Fusion Finite Volume Schemefrom these interpolated values define interpolated values for v and thesound speed. This also gives interpolated linearized system matrices.5. We average the eigenvalues at adjacent cells in the extrapolated set ofu-values to determine the sign of the eigenvalues at the cell edges (andhence the direction).6. We project the differences onto the eigenvectors of the interpolatedmatrices.7. We obtain limiters at half-integer index positions by projecting the up-wind difference onto the difference and applying the minmod routine.This allows for us to compute the limited high resolution correction.8. We calculate the upwind u-value and corresponding upwind flux athalf-integer indices.9. We neglect the geometric source term and update u with the upwindflux the limited high-resolution correction by taking a step k forwardin time.10. We then step forward by k in time with only the geometric sourceterms and no flux, where the geometric source terms are based on thevalues updated above.11. We update the wall positions based on the left and right velocity beforethe time-step, then update the width of the domain and new mesh.12. We then update the pressure at the left and right boundary based onthe new domain. The corresponding density at the left and right ofthe domain are given from the pressures.134Appendix A. Detailed Pseudo-Code of Fusion Finite Volume Scheme13. We then update velocity by v = ρv/ρ and the relative velocity of theleft and right boundaries.135


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items