Computational modeling and simulation of nonlinear electromagnetic and multiphysics problems

In this dissertation, nonlinear electromagnetic and multiphysics problems are modeled and simulated using various three-dimensional full-wave methods in the time domain. The problems under consideration fall into two categories. One is nonlinear electromagnetic problems with the nonlinearity embedded in either the permeability or the conductivity of the material's constitutive properties. The other is multiphysics problems that involve interactions between electromagnetic and other physical phenomena. A numerical solution of nonlinear magnetic problems is formulated using the three-dimensional time-domain finite element method (TDFEM) combined with the inverse Jiles-Atherton vector hysteresis model. A second-order nonlinear partial differential equation (PDE) that governs the nonlinear magnetic problem is constructed through the magnetic vector potential in the time domain, which is solved by applying the Newton-Raphson method. To solve the ordinary differential equation (ODE) representing the magnetic hysteresis accurately and efficiently, several ODE solvers are specifically designed and investigated. To improve the computational efficiency of the Newton-Raphson method, the multi-dimensional secant methods are incorporated in the nonlinear TDFEM solver. A nonuniform time-stepping scheme is also developed using the weighted residual approach to remove the requirement of a uniform time-step size during the simulation. Breakdown phenomena during high-power microwave (HPM) operation are investigated using different physical and mathematical models. During the breakdown process, the bound charges in solid dielectrics and air molecules break free and are pushed to move by the Lorentz force produced by the electromagnetic fields. The motion of free electrons produces plasma currents, which generate secondary electromagnetic fields that couple back to the externally applied fields and interact with the free electrons. When the incident field intensity is high enough, this will lead to an exponential increase of the charged particles known as breakdown. Such a process is first described by a nonlinear conductivity of the solid dielectric as a function of the electric field to model the dielectric breakdown phenomenon. The air breakdown problem encountered with HPM operation is then simulated with the plasma current modeled by a simplified plasma fluid equation. Both the dielectric and air breakdown problems are solved with the TDFEM together with a Newton's method, where the dielectric breakdown is treated as a pure nonlinear electromagnetic problem, while the air breakdown is treated as a multiphysics problem. To describe the plasma behavior more accurately, the plasma density and velocity are modeled by the equations of diffusion and motion, respectively. This results in a multiphysics and multiscale system depicted by the nonlinearly coupled full-wave Maxwell and plasma fluid equations, which are solved by a nodal discontinuous Galerkin time-domain (DGTD) method in three dimensions. The air breakdown during the HPM operation and the resulting plasma formation and shielding are modeled and simulated. Several important numerical issues in the simulation of nonlinear electromagnetic and multiphysics problems have been investigated and discussed. A continuity-preserving and divergence-cleaning scheme for electromagnetic problems involving inhomogeneous materials has been proposed based on the purely and damped hyperbolic Maxwell equations. A divergence-cleaning method is presented to enforce Gauss's laws and normal flux continuity by introducing auxiliary variables and damping terms into the original Maxwell's equations, which result in artificial propagation and dissipation of the numerical errors. Based on the DGTD method, dynamic h- and p-adaptation algorithms are developed for a full-wave analysis of electromagnetic and multiphysics problems. The dynamic h-adaptation algorithm can dynamically refine the mesh to resolve the local variation of the fields during the wave propagation, while the dynamic p-adaptation algorithm can determine and adjust the basis order in real time during the simulation. Both algorithms developed and investigated in this dissertation are highly flexible and efficient, and are powerful simulation tools in the solution of nonlinear electromagnetic and multiphysics problems

Fourth-order finite volume algorithm with adaptive mesh refinement in space and time for multi-fluid plasma modeling, A

2022 Spring.Includes bibliographical references.Improving our fundamental understanding of plasma physics using numerical methods is pivotal to the advancement of science and the continual development of cutting-edge technologies such as nuclear fusion reactions for energy production or the manufacturing of microelectronic devices. An elaborate and accurate approach to modeling plasmas using computational fluid dynamics (CFD) is the multi-fluid method, where the full set of fluid mechanics equations are solved for each species in the plasma simultaneously with Maxwell's equations in a coupled fashion. Nevertheless, multi-fluid plasma modeling is inherently multiscale and multiphysics, presenting significant numerical and mathematical stiffness. This research aims to develop an efficient and accurate multi-fluid plasma model using higher-order, finite-volume, solution-adaptive numerical methods. The algorithm developed herein is verified to be fourth-order accurate for electromagnetic simulations as well as those involving fully-coupled, multi-fluid plasma physics. The solutions to common plasma test problems obtained by the algorithm are validated against exact solutions and results from literature. The algorithm is shown to be robust and stable in the presence of complex solution topology and discontinuities, such as shocks and steep gradients. The optimizations in spatial discretization provided by the fourth-order algorithm and adaptive mesh refinement are demonstrated to improve the solution time by a factor of 10 compared to lower-order methods on fixed-grid meshes. This research produces an advanced, multi-fluid plasma modeling framework which allows for studying complex, realistic plasmas involving collisions and practical geometries

Numerical investigation of MPD thrusters using a density-based method with semi-discrete central-upwind schemes for MHD equations

The magnetohydrodynamic (MHD) equations which combines the Navier-Stokes equations with the Maxwell equations are essential for the investigations of many research areas as earth's core modelling, metal casting, fusion devices and electrical and aerospace devices. In the present work, the central-upwind schemes proposed by Kurganov, Noelle and Petrova for hydrodynamics are extended and combined with the divergence cleaning method of Dedner in order to investigate the performance of the self-field and applied magneto-plasma dynamic thrusters which still involving some outstanding problems. This new algorithm is developed for the single temperature, ideal and resistive MHD equations in a finite volume discretization framework with Gaussian integration. The electrical conductivity is predicted according to the Spitzer-Härm formulation and the real gas ratio of specific heats proposed by Sankaran is implemented for higher discharge current. To improve the quality of the solution, the limiter function of first and second order interpolation scheme is used. The accuracy and the robustness of the obtained solver are demonstrated through numerical simulations of ideal MHD benchmark problems. first, the ability of the developed code to handle shocks, rarefactions and contact discontinuities is tested with the Brio-Wu shock-tube problem. The Minmod and the Van Albada limiter functions has been found to perform better than the other limiter used and the obtained results agree well with both the analytic and the simulations results of previous work. Secondly, the complex and multiple shock interactions and the transition from smooth to turbulent flow involved in the Orszag-Tang vortex problem is well described by the present code and the comparison with the WENO-5 scheme of Shen shows good agreement. Lastly, The ability to described the interaction of an denser cloud with a MHD shock is tested by simulating the 2D cloud-shock interaction problem. The main phases of the interaction are well captured by the solver and the temporal progression of the density contour is in accordance with those obtained by Xisto. The ability of the developed resistive solver to deal with plasma flow acceleration is tested by simulating the well experimental investigated thrusters: The full scale benchmark thruster and the extended anode thruster of Princeton. The results show good agreement with the experimental and simulations results of previous work for discharge current less than the critical current just before the beginning of the onset phenomenon. Simulations are also conducted on the Villani-H thruster to determine the effect of geometric changes over the thruster performance and a first designing attempts is proposed according to the stability analysis. Confident with the results obtained with ideal and resistive MHD problems, the present code is extended to applied-field MPD thrusters. The purpose is to achieve a high thrust level required for space missions with less input power than with self-field MPD thrusters and thus avoid the onset instabilities. For the verification of the code, the NASA Lewis Research Center's (NASALeRC) MPD thruster is chosen because of its wide range of experimental data bank. The method presented reproduce the theory of thrust production and plasma acceleration. Some difficulties as the limitation of the maximum rotational speed and the depletion of the plasma density on the anode surface have been captured. Moreover, the present density-based method compares very well with experimental data of Myers and simulations of Mikellides

Numerical resistive relativistic magnetohydrodynamics

La presente tesis se desarrollada dentro del marco de la Magnetohidrodinámica Resistiva Relativistica (RRMHD; por sus siglas en inglés) y uno de sus principales objetivos es el caracterizar las condiciones físicas que optimizan la disipación de campos magnéticos en plasmas relativistas, especialmente en aquellos que son de interés astrofísico. Para alcanzar este objetivo, realizamos el estudio de los denominados modos de ruptura dobles, bajo condiciones ideales (IDTMs; por sus siglas en inglés) que maximizan sus tasas de crecimiento. Se demuestra que en el régimen relativista los IDTMs pueden crecer en escalas de tiempo de unos pocos tiempos característicos de Alfvén, desarrollando regímenes explosivos, incluso si las condiciones en las que se desarrollan no son estrictamente ideales. Ello nos permite concluir que los IDTMs relativistas pueden utilizarse para explicar los fenómenos de reconexión astrofísicos más violentos, como por ejemplo los que se cree que acontecen en la magnetosferas de estrellas de neutrones. Para lograr este objetivo, piedra angular de la tesis, se construyó un nuevo código RRMHD, apto para ser usado en el estudio de plasmas en el contexto astrofísico. Este nuevo código denominado CUEVA , se basa en una formulación conservativa de volúmenes finitos de las ecuaciones de RRMHD. La evolución de un estado inicial dado se realiza mediante la técnica conocida como método de líneas. Dado que en el régimen ideal el sistema de ecuaciones de la RRMHD es matemáticamente rígido, la integración temporal se lleva a cabo con métodos parcialmente implícitos. En C UEVA se implementan dos familias principales de integradores de tiempo: los métodos denominados RKIMEX y MIRK. Como un subproducto de esta tesis, hemos desarrollado un resolvedor aproximado del tipo HLLC. El nuevo resolvedor captura de forma exacta discontinuidades de contacto estacionarias. En combinación con técnicas de reconstrucción espacial de orden ultra-alto, en esta tesis caracterizamos la resistividad y la viscosidad numérica de CUEVA de manera exhaustiva. Ello nos ha permitido delimitar con claridad que el desarrollo de los IDTMs es de origen físico y no un artefacto numérico.The main goal of this thesis is the study of magnetic reconnection in relativistic plasma of astrophysical interest. We pay special attention to the dynamical effects that the resistive dissipation of magnetic fields may have in such process. Thus, our approach to the study is numerical, i.e. we developed numerical models that mimic as closely as possible the physical conditionsunder which reconnection happens in relativistic astrophysical plasma. The goals of this thesis can be grouped in two sets: computational and physical. Considering that we aim to obtain physical results (ideally) indepen-dent of the numerical methods that we employ, we pursue the development of a new multidimensional RRMHD code for astrophysical applications. The code include, different numerical algorithms for the time-evolution, for the solution of the Riemann problem and for the intercell reconstruction. Only in this way we may calibrate the impact of the numerical methods on the results. Therefore, concerning the physics of reconnection in astro-physical plasma, our first and foremost goal is to understand the dependenceof the growth rate of resistive instabilities (such as TM instabilities) on the physical properties of the plasma. We shall consider specifically the case of relativistic ideal double tearing modes (RIDTMs), where the adjective ideal refers to the fact that the instability develops at a timescale of the order of the MHD timescale (i.e. Alfvén crossing time) and its growth rate is independent of the resistivity (and thus, it does not diverge in the ideal RMHD regime). This is an unexplored and interesting setup where two parallel current sheets may interact as TMs develop yielding an explosive reconnection episode, where the reconnection time scales are so small that may be of interest for explaining a number of astrophysical objects

Discretisations and Preconditioners for Magnetohydrodynamics Models

The magnetohydrodynamics (MHD) equations are generally known to be difficult to solve numerically, due to their highly nonlinear structure and the strong coupling between the electromagnetic and hydrodynamic variables, especially for high Reynolds and coupling numbers. In the first part of this work, we present a scalable augmented Lagrangian preconditioner for a finite element discretisation of the $\mathbf{B}$-$\mathbf{E}$ formulation of the incompressible viscoresistive MHD equations. For stationary problems, our solver achieves robust performance with respect to the Reynolds and coupling numbers in two dimensions and good results in three dimensions. Our approach relies on specialised parameter-robust multigrid methods for the hydrodynamic and electromagnetic blocks. The scheme ensures exactly divergence-free approximations of both the velocity and the magnetic field up to solver tolerances. In the second part, we focus on incompressible, resistive Hall MHD models and derive structure-preserving finite element methods for these equations. We present a variational formulation of Hall MHD that enforces the magnetic Gauss's law precisely (up to solver tolerances) and prove the well-posedness of a Picard linearisation. For the transient problem, we present time discretisations that preserve the energy and magnetic and hybrid helicity precisely in the ideal limit for two types of boundary conditions. In the third part, we investigate anisothermal MHD models. We start by performing a bifurcation analysis for a magnetic Rayleigh--B\'enard problem at a high coupling number $S=1{,}000$ by choosing the Rayleigh number in the range between 0 and $100{,}000$ as the bifurcation parameter. We study the effect of the coupling number on the bifurcation diagram and outline how we create initial guesses to obtain complex solution patterns and disconnected branches for high coupling numbers.Comment: Doctoral thesis, Mathematical Institute, University of Oxford. 174 page

Generalized averaged Gaussian quadrature and applications

A simple numerical method for constructing the optimal generalized averaged Gaussian quadrature formulas will be presented. These formulas exist in many cases in which real positive GaussKronrod formulas do not exist, and can be used as an adequate alternative in order to estimate the error of a Gaussian rule. We also investigate the conditions under which the optimal averaged Gaussian quadrature formulas and their truncated variants are internal

MS FT-2-2 7 Orthogonal polynomials and quadrature: Theory, computation, and applications

Quadrature rules find many applications in science and engineering. Their analysis is a classical area of applied mathematics and continues to attract considerable attention. This seminar brings together speakers with expertise in a large variety of quadrature rules. It is the aim of the seminar to provide an overview of recent developments in the analysis of quadrature rules. The computation of error estimates and novel applications also are described

High-Order Particle Integration for Particle-In-Cell Schemes using Boris with Spectral Deferred Corrections

The study of plasmas plays an important role in both science and technology. Plasma dynamics can be found wherever charged particles or materials interact with and generate electromagnetic fields, covering more orders of magnitude in scale and density than any other type of matter. Plasma phenomena dominate the dynamics of the sun, stars and space between them and are important to a variety of technologies, from fusion reactors to spacecraft propulsion. As plasma behaviour and associated mathematical relations are naturally complex, numerical methods and computer simulation play a crucial role in furthering the field. Particle-in-Cell (PIC) is a class of numerical scheme currently used in the simulation of hot diffuse plasmas, or denser plasmas at small scales. One crucial part of such schemes is the particle integrator, which solves the particle equations of motion, typically via time discretisation of the Newton-Lorentz force. For nearly forty years, the dominant algorithm for charged particle tracking has been leapfrog integration using Boris' algorithm. The combined scheme is often referred to simply as the classic Boris integrator and provides a directly computable centre-difference discretisation for the implicit system. As the Boris algorithm is intrinsically second order accurate, a tunable order algorithm based on Boris and spectral deferred corrections (Boris-SDC) was recently proposed and demonstrated to exhibit high order time convergence for a single-particle Penning trap with exactly known electromagnetic fields. The faster reduction in error as time-step size is decreased allowed Boris-SDC to be more computationally efficient than classic Boris, but whether the advantageous characteristics of Boris-SDC would extend to PIC and approximated fields was not investigated. This thesis contributes the implementation and performance testing of Boris-SDC within PIC schemes and generalises Boris-SDC to the relativistic regime. This relativistic extension to Boris-SDC is shown to retain higher order time convergence and improved computational performance when compared to classic Boris even in highly relativistic regimes ($>99\%$ speed of light). The relativistic Boris-SDC integrator is shown to produce less unphysical drift than classic Boris in the force-free scenario where electric and magnetic forces cancel. The algorithmic modifications required to implement Boris-SDC within PIC are highlighted first and the impact of spatial electric field approximation on particle integrator performance is demonstrated for the electrostatic case (ESPIC). The relativistic Boris-SDC integrator is then derived and implemented within the open-source PIC code Runko, demonstrating capability of Boris-SDC to work with existing codes. Finally, performance tests are conducted in the form of work-precision comparisons to classic Boris for two electrostatic benchmarks, the two-stream instability and Landau damping. The spatial field approximation inherent to ESPIC imposes an error saturation on the time convergence of the particle integrator, inversely proportional to the spatial resolution, which limits the achievable global error. The limited accuracy is found to erode the computational performance of Boris-SDC, as the low level of error required to offset the added computational cost of SDC cannot be reached. Above the spatial saturation point however, Boris-SDC is found to retain high order time convergence and higher accuracy for a fixed time-step size than classic Boris. As a final note, suggestions are given for further work, including use of higher order spatial methods, investigation on the significance of momentum error vs. spatial error as well as PIC applications wherein Boris-SDC might be useful despite the lack of a clear performance gain