Basic principles of numerical radiation hydrodynamics and the RALEF code
Sponsored Links
This presentation is the property of its rightful owner.
1 / 87

Basic principles of numerical radiation hydrodynamics and the RALEF code PowerPoint PPT Presentation

  • Uploaded on
  • Presentation posted in: General

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. Main constituents of the RALEF-2D package.

Download Presentation

Basic principles of numerical radiation hydrodynamics and the RALEF code

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript

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

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


Thermal conduction

Radiation transport

EOS and opacities

Laser absorption

– 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

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


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”.

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!

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.

Computational mesh: block structure

Every block can be subdivided into np1np2 topologically rectangular parts;

different parts can be composed of different materials.

RALEF: the ALE technique in action

Ncyc = 1

Ncyc = 100

Limitations of the ALE technique

Example: laser irradiation of a Cu foil

t = 1.6 ns

fixed-pressure boundary

Cu foil


Simulation stops because a singularity develops at the plasma-vacuum boundary !

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.

Laser irradiation of a thin Sn disk

free-outflow boundary

reflective (symmetry) boundary

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.



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 !

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.

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.

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 !

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

Thermal conduction

(a test bed for radiation transport)

– 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):

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.

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.

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.

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.

Test problems: non-linear wave into a cold wall


Square grid 100x100

Initial condition:

Boundary condition:


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

Standard grids for test problems

Only the Kershaw grid is strongly distorted in the sense that some of the coefficients (1)(1) become negative.

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).

Radiation transport

– 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):

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 !

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.

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.

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 .


  • 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.


  • a significant amount of numerical diffusion in space.

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 !

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:

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!


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.


Our numerical scheme

  • is not conservative (violates the finite-difference version of the Gauss theorem), and

  • hasno strict convergence in the diffusion limit.

Radiation transport:

test problems

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.

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%.

Diffusion limit in a slab

No convergence, the error level depends on the level of mesh distortion.

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

Fireball expansion in xy-geometry

By t = 0.002102926 the exact position of the front must be at r = 1.5.

Final state: the conduction solution

By t = 0.002102926 the exact position of the front must be at r = 1.5.

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.

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% !

SN method: the “ray effect”

Consider radiation transport from a central “hot” rod across a vacuum cavity.


Convergence of the Sn method


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)

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:

Illustrative problems

Importance of radiation transport for laser plasmas


Sn foil:Δzfoil = 0.42 µm, Rfoil = 75 μm, Rmesh = 120 μm.

CO2 laser

Peak absorbed laser flux: F00= 1.41010 W/cm2


Sn foil

Ring-Gaussian spatial profile


Temporal profile (normalized)

No radiative transfer:

thermal conduction only

Radiation transport:

S8with 8 spectral groups

Irradiation by a short (20ps) YAG ( = 1.064 μm)pulse

Sn microsphere: RSn = 17 μm, Rmesh = 150 μm.


Nominal laser energy: Elas-absorbed = 1 mJ

(reflection ignored)


Absorbed laser energy: Elas-absorbed = 0.4 mJ

 0.6 mJ missed the target

Gaussian temporal profile (normalized)

Gaussian spatial profile (normalized)

Sn droplet irradiated by CO2-laser ( = 10.6 μm)

Sn microsphere: RSn = 17 μm, Rmesh = 150 μm.


Absorbed laser energy:Elas-absorbed = 10 mJ


Peak absorbed laser flux: F00= 6.2109W/cm2

Gaussian temporal profile (normalized)

Gaussian spatial profile (normalized)

Spectral Sn opacities

In band: 2 -groups

Click to see the movie: Time is in nanoseconds

Evolution after the laser is turned off: equilibrium EOS

The effect of metastable EOS

Fully equilibrium EOS

(Maxwell construction in the two-phase region)

Metastable EOS

(van der Waals loops in the two-phase region)


  • Our agenda for future work:

  • laser deposition with refraction and reflection,

  • hydrodynamics with metastable → stable phase transitions,

  • material mixing,

  • non-LTE radiation transport.

  • Login