Basic principles of numerical radiation hydrodynamics and the RALEF code

1. Basic principles of numerical radiation hydrodynamics and the RALEF code M. M. Basko Keldysh Institute of Applied Mathematics, Moscow, Russia ELI, Prague 10 July 2014

2. Main constituents of the RALEF-2D package The RALEF-2D code has been developed with the primary goal to simulate high-temperature laser plasmas. Its principal constituent blocks are Hydrodynamics Thermal conduction Radiation transport EOS and opacities Laser absorption

3. – energy deposition by thermal conduction (local), – energy deposition by radiation (non-local), – eventual external heat sources. Equations of hydrodynamics The RALEF-2Dcode is based on a one-fluid, one-temperature hydrodynamics model in two spatial dimensions (either x,y, or r,z): ideal hydrodynamics

4. Selection of principal dependent variables Illustration with the 1D planar example: Divergent form: ρ, ρu, ρE Non-divergent form: ρ, u, e   Naturally leads to a conservative numerical scheme Easier to avoid unphysical values of the internal energy e and temperature T RALEF

5. Selection of principal independent variables Divergent Eulerian form: t, x Non-divergent Lagrangian form: t, m Eulerian simulation: spatial mesh in x,y,zis fixed in time. Lagrangian simulation: Lagrangian mesh in m=∫ρdx is “attached” to the fluid; typically, the simulation time is severely limited by mesh “collapse”.

6. The ALE technique (adaptive mesh) The Arbitrary Lagrangian-Eulerian approach allows free motion of the x,y,zmesh – independent of the motion of the fluid! Implementation in RALEF: every time step consists of 2 substeps: a purely Lagrangian step – the mesh is “attached” to the fluid, followed by a rezoning step, where a completely new x,y,z mesh is constructed, and the principal dynamic variables are remapped to the new mesh. At substep 2, a new mesh is constructed by solving the Dirichlet problem for a separate Poisson-like differential equation with a user-defined weight function. In RALEF, there are quite a number of user-accessible parameters for control over the mesh dynamics – which is often the most tricky part of a successful simulation!

7. Mesh topology in RALEF Computational mesh consists of blocks with common borders. Each block is topologically equivalent to a rectangle. Common block faces have equal number of cells. Mesh geometry: • Cartesian (x,y) • cylindrical (r,z) Mesh library: Different mesh options, distinguished by the value of variable igeom, are combined into a mesh library, which is being continuously expanded.

8. Computational mesh: block structure Every block can be subdivided into np1np2 topologically rectangular parts; different parts can be composed of different materials.

9. RALEF: the ALE technique in action Ncyc = 1 Ncyc = 100

10. Limitations of the ALE technique Example: laser irradiation of a Cu foil t = 1.6 ns fixed-pressure boundary Cu foil laser Simulation stops because a singularity develops at the plasma-vacuum boundary !

11. Boundary conditions • To obtain a particular solution of the hydrodynamics equations, we need • to define the initial state, and • to impose certain boundary conditions. The initial state is conceptually simple and fully free to be set by the user: either equilibrium or non-equilibrium. Setting up a physically consistent and computationally efficient combination of boundary conditions is a considerably more complex task (requires basic understanding of hydrodynamics and some physical intuition). • Basic types of boundary conditions in RALEF: • reflective (symmetry) boundary (fixed in space → Eulerian) • prescribed pressure (moves in space → Lagrangian) • prescribed velocity (moves in space → Lagrangian) • free outflow (fixed in space → Eulerian) • prescribed inflow (fixed in space → Eulerian) Additional degree of freedom: the boundary conditions and the ALE mode can be changed at will at any time during the simulation.

12. Laser irradiation of a thin Sn disk free-outflow boundary reflective (symmetry) boundary

13. Different materials By finding an appropriate combination of boundary conditions and ALE options, one can adequately simulate practically any 2D problem with a single material. Multiple materials pose an additional challenge: In RALEF, mixing of different materials within a single mesh cell is not allowed  hence, any material interface must be treated as a Lagrangian curve, which usually tends to get tangled: as a result, the simulation stops. Al Sn

14. Equation of state For ideal hydrodynamics (without thermal conduction, viscosity, etc) we need an equation of state in the form RALEF can easily accommodate any thermodynamically stable EOS. However, because the principal variable is E=e+u2/2 rather than e (or T), the numerical scheme is not positive with respect to T sporadic appearance of negative temperatures! • An additional conceptual difficulty arises for a van-der-Waals type of EOS in the region of liquid-vapor phase coexistence – i.e. below the binodal, where we have two different branches of the same EOS: • a metastable branch, and • an absolutely stable equilibrium branch, obtained by means of Maxwell construction. This latter problem has not been resolved yet !

15. Implicit and explicit numerical hydrodynamics Implicit scheme Explicit scheme Implicit schemes are numerically stable and allow large time steps, but are difficult to implement. Numerical stability of explicit schemes requires the Courant–Friedrichs–Lewy condition (CFL condition) to be fulfilled which for certain problems may incur prohibitively long computing times. RALEF has explicit hydrodynamics, i.e. is oriented on simulating fast processes.

16. Numerical scheme for hydrodynamics The numerical scheme for the 2D hydrodynamics is built upon the CAVEAT-2D (LANL, 1990) hydrodynamics package and has the following properties: • it uses cell-centered principal variables on a multi-block structured quadrilateral mesh (either in the x-y or r-z geometry); • is fully conservative and belongs to the class of second-order Godunov schemes (no artificial viscosity is needed); • the mesh is adapted to the hydrodynamic flow by applying the ALE (arbitrary Lagrangian-Eulerian) technique; • the numerical method is based on a fast non-iterative Riemann solver (J.K.Dukowicz, 1985), well adapted for handling arbitrary equations of state.

17. Principal variables: Equation of state: i, ui, Ei The Godunov numerical method ― all assigned to the cell centers ! Lagrangian phase: the Riemann problem is solved at each cell face  fluxes F of momentum uand total energy E  new ui(t+dt)and Ei(t+dt); mass is conserved. ← 1st order A fast non-iterative Riemann solver by J.K.Dukowicz (1985) is used. A “snag”: node velocities uvi! No artificial viscosity is needed !

18. Importance of the 2-nd order + ALE Non-linear stage of the Rayleigh-Taylor instability of a laser-irradiated thin carbon foil RALEF: 1-st order RALEF: 2-nd order

19. Thermal conduction (a test bed for radiation transport)

20. – energy deposition by thermal conduction (local), – energy deposition by radiation (non-local), – eventual external heat sources. Equations of hydrodynamics The newly developed RALEF-2D code is based on a one-fluid, one-temperature hydrodynamics model in two spatial dimensions (either x,y, or r,z):

21. Implicit versus explicit algorithms for conduction Fully implicit scheme Explicit scheme An explicit scheme is stable under the condition , which practically always incurs prohibitively long computing times. A fully implicit scheme is always stable but requires iterative solution – which can be implemented for thermal conduction, but becomes absolutely unrealistic for radiation transport. Symmetric semi-implicit (SSI) scheme RALEF is based on the SSI algorithm for both the thermal conduction and radiation transport.

22. The SSI method for thermal conduction The key ingredient to the RALEF-2D code is the SSI (symmetric semi-implicit) method of E.Livne & A.Glasner (1985), used to incorporate thermal conduction and radiation transport into the 2D Godunov method (in order to avoid costly matrix inversion required by fully implicit methods). The numerical scheme for thermal conduction (M.Basko, J.Maruhn & A.Tauschwitz, J.Com.Phys., 228, 2175, 2009) has the following features: • it uses cell-centered temperatures from the FVD (finite volume discretization) hydrodynamics on distorted quadrilateral grids, • is fully conservative (based on intercellular fluxes with an SSI energy correction for the next time step), • (almost) unconditionally stable, • space second-order accurate on all grids for smooth , • symmetric on a local 9-point stencil, • computationally efficient.

23. Time step control in the SSI method Constraints on t are based on usual approximation-accuracy considerations, but here we need two separate criteria with two control parameters: Ts is a problem-specific sensitivity threshold for temperature variations.

24. Test problems: linear steady-state solutions Time convergence to the steady state The linear solution is reproduced exactly on all grids, unless special constraints to ensure positiveness are imposed on stronglydistorted grids. For sufficiently large time steps t, the SSI method does not converge to the steady-state solution  an indication of its conditional stability (neutral for large t). Piecewise linear solutions with -jumps are reproduced exactly only with harmonic mean f,ij , and only on rectangular grids.

25. Test problems: non-linear wave into a cold wall Parameters: Square grid 100x100 Initial condition: Boundary condition: Solution: Test results: This solution cannot be simulated with the harmonic-mean f , whereas excellent results are obtained with the arithmetic-meanf : the front position is reproduced with an error of 0.1% for  . Temperature in all cells

26. Standard grids for test problems Only the Kershaw grid is strongly distorted in the sense that some of the coefficients (1)(1) become negative.

27. Test problems: non-linear steady-state solution Here we test against a steady-state solution of the form with the source term . Our scheme clearly demonstrates the 2nd order convergence rate on all grids, and is no less accurate than the best published schemes of Morel et al. (1992), Shashkov et al. (1996), Breil & Maire (2007).

28. Radiation transport

29. – energy deposition by thermal conduction (local), – energy deposition by radiation (non-local), – eventual external heat sources. Equations of hydrodynamics The newly developed RALEF-2D code is based on a one-fluid, one-temperature hydrodynamics model in two spatial dimensions (either x,y, or r,z):

30. Radiation transport Transfer equation for radiation intensity I in the quasi-static approximation (the limit of c → ): Quasi-static approximation: radiation transports energy infinitely fast (compared to the fluid motion)  the energy residing in radiation field at any given time is infinitely small ! In the present version, the absorption coefficient k and the source function B= B(T) are calculated in the LTE approximation. Coupling with the fluid energy equation: Radiation transport adds 3 extra dimensions (two angles and the photon frequency) the 2D hydrodynamics becomes a 5D radiation hydrodynamics !

31. Not to be mixed up with spectral multi-group diffusion In many cases the term “radiation hydrodynamics” (RH) is applied to hydrodynamic equations augmented with the multi-frequency diffusion equation for the spectral radiation energy density ; the coupling term to the fluid is Here the information about the angular dependence of the radiation field is lost; one simply has to solve some 30 – 100 additional mutually independent diffusion equations. Not much of a challenge for computational physics: there already exist numerically stable, positive and conservative numerical schemes on distorted (non-rectangular) grids. In the RALEF code we have such a scheme implemented for the thermal conduction.

32. Challenges for computational RH Ideally, the numerical scheme for solving the radiation transfer equation must possess the following important properties: • it must be positive in the sense that non-negative boundary values of Iν guarantee Iν≥ 0at all collocation points ― provided that Bν ≥ 0 ; • it must be conservative in the sense that both locally and globally the finite-difference equivalent of the Gauss theorem is fulfilled, • it must reproduce the diffusion limit on arbitrary quadrilateral (or triangular) grids in optically thick regions ― especially when the optical thickness of individual grid cells becomes >> 1  a high-order accuracy scheme is needed. The difficulty of constructing such a scheme has been recognized early in the theory of neutron transport: K.D.Lathrop, J.Comp.Phys., 4 (1969) 475. To the best of our knowledge, no such scheme is known in the open literature.

33. Calculation of Iν: the method of short characteristics • angular directions are discretized by using the Sn method with n(n+2) fixed photon propagation directions over the 4π solid angle; • for each angular direction and frequency , the radiation field I is found by solving the transfer equation with the method of short characteristics (A.Dedner, P.Vollmöller, JCP, 178, 263, 2002). Mesh nodes are chosen as collocation points for the radiation intensity I . Advantages: • even on strongly distorted meshes, it is guaranteed that light rays pass through each mesh cell; • the algorithm is generally computationally more efficient than that with long characteristics. Disadvantages: • a significant amount of numerical diffusion in space.

34. Conservativeness: flux conservation in vacuum Consider a scheme where the radiation intensity is assigned to mesh vertices. Consider one cell with a negligibly small absorption (k 0): In our example on the left we have Conclusion: a numerical scheme based on nodal collocation points is generally non-conservative ! Conservativeness can be restored by assigning the intensity I to cell faces (R.Ramis, 1992) ― equivalent to using the fluxes H as independent variables But then the high-order of finite-difference approximation, needed for the diffusion limit, is lost !

35. The diffusion limit The diffusion limit is the limit of k>> | lnT|, | ln|, | lnk | In this limit we can apply the asymptotic expansion: In the quasi-static LTE case the diffusion limit is equivalent to the approximation of radiative thermal conduction:

36. Consider two opposite beamlets propagating along an infinitely thin column h: Heating by two opposite beamlets in the diffusion limit In the diffusion limit we have Having combined the two opposite beams, we get Conclusion: our finite-difference scheme must correctly reproduce the 2nd spatial derivatives of the unknown function Iν– i.e. to be of high order accuracy! The latter is a major challenge for non-rectangular grids in 2 and 3 dimensions – and especially so when the cell optical thickness Δτ >> 1!

37. Achievements We have developed an original numerical scheme for radiation transport which • is strictly positive, • from the practical point of view ― reproduces the diffusion limit with a fair degree of accuracy on not too strongly distorted meshes, • works robustly in both (x,y) and (r,z) geometries. Failures Our numerical scheme • is not conservative (violates the finite-difference version of the Gauss theorem), and • hasno strict convergence in the diffusion limit.

38. Radiation transport: test problems

39. Numerical diffusion: a searchlight beam in vacuum The short-characteristic method produces a significant amount of numerical diffusion for light beams with sharp edges. For thermal radiation a certain amount of numerical diffusion may be more an advantage than a drawback.

40. Test problems: radiative cooling of a slab Computational mesh: Exact solution: Convergence of the Sn method Typical accuracy in problems with optically thin mesh cells is 12%.

41. Diffusion limit in a slab No convergence, the error level depends on the level of mesh distortion.

42. Test “fireball” (the limit of radiative heat conduction) Consider an initially hot unit sphere with T0 = 1 and k = 100/T. Propagation of a radial heat wave is governed by the conduction equation which admits an analytical self-similar solution with the front position

43. Fireball expansion in xy-geometry By t = 0.002102926 the exact position of the front must be at r = 1.5.

44. Final state: the conduction solution By t = 0.002102926 the exact position of the front must be at r = 1.5.

45. Final state: the transport solution with k0 = 100 From practical point of view, in this particular case we observe a slow convergence to the exact diffusion solution.

46. Final state: the transport solution with cell = 10 In this example we observe no convergence to the exact diffusion solution. Numerical energy “leak” is  40% !

47. SN method: the “ray effect” Consider radiation transport from a central “hot” rod across a vacuum cavity. S6 Convergence of the Sn method

48. Opacities

49. Opacity options in RALEF-2D Here we profit from many years of a highly qualified work at KIAM (Moscow) in the group of Nikiforov-Uvarov-Novikov (the THERMOS code based on the Hartree-Fock-Slater atomic modeling). Opacity options: • power law, • ad hoc analytical, • inverse bremsstrahlung (analytical), . . . . . . . • GLT tables (source opacities from Novikov)

50. Laser absorption in RALEF-2D • An arbitrary number of monochromatic cylindrical laser beams is allowed: for every beam the transfer equation is solved by the same method of short characteristics as for thermal radiation. • No refraction, no reflection, the inverse bremsstrahlung absorption coefficient, artificially enhanced absorption beyond the critical surface. Present model: