##### Document Text Contents

Page 1

Saas-Fee Advanced Course 27

Lecture Notes 1997

Page 2

Springer

Berlin

Heidelberg

New York

Barcelona

Budapest

Hong Kong

London

Milan

Paris

Singapore

Tokyo

Page 261

9. Numerical Radiation Hydrodynamics 247

9.2 Computational Strategy

Overall Procedure. It is costly to try to solve directly the full angle-frequency

dependent radiation equations coupled to the equations of hydrodynamics be-

cause the dimensionality of the problem (= 4) is too large. So we look for

an approximate technique. Notice that we need only the frequency-integrated

moments to solve the material energy and momentum equations. Therefore

we are led to try a splitting technique in which the spacetime evolution is han-

dled separately from the angle-frequency information. The method described

below is called the multifrequency/grey method. For ease of exposition, as-

sume we use explicit hydro.

1. Suppose we are given the solution at t^. Update velocities to t^~^^ (ma-

terial momentum equation), hence radii and densities to t^^^.

2. If we can assume that the spectral distributions (Ejy/E) and {F^IF) are

known at t̂ "*"̂ (actually they are not!) we can calculate K^ and XF?

hence A:E = f^E/f^p and fcp = X F / X R - Similarly, if we can assume that

the angular distribution of the radiation field is known at t'^'^^ (again, it

isn't), then we have /^ and g,,, hence / and q.

3. At t^'^^ we can now solve the gas energy equation plus the radiation en-

ergy and momentum equations. We thus obtain the run of T^-^^^E^^^,

and F^"*" .̂ This system is nonlinear because the material properties de-

pend nonlinearly on T""^^; to solve it we linearize and iterate to conver-

gence, as described earlier.

In principle, these steps give us the solution at the advanced time. But ac-

tually they don't because we assumed we already knew (Ej^/E), {Fty/F),fi^,

and qiy at P"^^, when in fact we didn't. Therefore at the advanced time level

p + i we should now:

4. Solve the angle-dependent transfer equation for the geometric factors

ft, and qj, using source-sink terms evaluated at t^^^ [i.e., KI,{T^'^^),

7?.(r»+i) ].

5. Then solve the frequency-dependent moment equations for the spectral

profiles {Ey/E) and (Fj^/F). Use these updated profiles to update the

information in ACE,XFJ^E, ^uid A:F at t^'^^.

6. Go to step 3, and iterate to consistency.

— In step 4, we ignore frequency-coupling in order to reduce the dimension-

ality of the problem to (r,/x) at a given (i^^t). From this step we derive the

geometric factors /^ and ^^ which depend mainly on angular information

about the radiation field.

— In step 5, we dispose of angular information by using moments, and attempt

to improve the frequency-dependent information contained in {Ej,/E) and

{Fu/F) by accounting for the (i/,r) coupling at a given ,̂ assuming that

the geometric factors {fi,,qi,) are known.

Page 262

248 D. Mihalas

- This procedure works fairly well because fy is only a ratio which is mainly

geometry dependent, whereas {E^jE) and (Fiy/F) are only spectral pro-

files, which are insensitive to geometry.

— Very few calculations have been done to this level of consistency. Most ei-

ther ignore the frequency-coupling in the moments, and lag the Eddington

factors, or they handle the frequency-couphng but assume /»/ = | (multi-

group diffusion).

Spatial Differencing and Unresolved Fronts. It is straightforward to write

difference formulae for the equations written above. We won't produce them

here; the reader can easily work them out. But consider one very difficult

problem: how can we calculate a good value for XF at an interface? A seem-

ingly reasonable value giving the correct optical depth increment between cell

centers is

QJd ^ " ^ d - l / 2 + ^ ^ d - f 1/2

But trouble arises when we fail to resolve a front in the flow (a shock, ion-

ization front, ...). The problem is that the opacity can vary as T" with

a ~ 1 2 — 1 4 i n regions where H and/or He get excited and/or ionize. There-

fore even small changes in T between zones imply a huge change in x- This

rapid variation in x makes it hard to choose the correct value at the inter-

face, and the errors made can clobber the energy transport at the interface.

Generally the result is disasterous. How can we fix this problem? R. Christy

[6] found that we can get "reasonable" results using an energy-weighted har-

monic mean:

(^''-/•^)^(0.-./. + (^^-/^)^(f). Q\ _ V-^/d-l/2 VX/d+i/2

This recipe has been very popular, but it is only ad hoc. The real solution

to the problem is to use an adaptive grid. The problem can't be solved by

brute force with tiny Lagrangian zones. For example, the ionization zone in

Cepheids and RR Lyrae envelopes is about 0.001 to 0.01 scale heights thick,

and moves back and forth over 3-5 scale heights during a pulsation cycle.

A prohibitive number of zones (and timesteps!) would be required to resolve

the front. In contrast, an adaptive grid moves with the front and resolves it.

This approach is the preferred technique.

Method of Solution. For the explicit hydro scheme sketched above, at t^^^ we

get equations coupling j'^^+i^ ^^+i ^ ^^^+1 at d = 1, ... , D. Linearize the equa-

tions around the current estimate of the solution. We get a block tridiagonal

system which can be solved by Gaussian elimination, and iterated to con-

vergence (Newton-Raphson procedure). In Unearizing opacities we can write

/̂ E = k^Kp and XF = ^ F X R ^^^ assume {k^.k^) remain fixed. Or (better)

Page 521

Index 507

- SNI 354,356,357

- SN la 356, 358, 359, 408, 412-415,

417,418

- SNIb 354-360,393

- SNIc 354-360,393

- SN II 354, 358-360, 376, 393

- SN II-L 354, 357, 360

- SNII-P 354,357,360

- SNIIb 357

- SNIIn 354

- SNIIP 376

- spectrum 356

- Type II 251

- Type I 353

- Type II 353,372

Synchrotron radiation transfer 462

TAMA 392

Tangent space 185,188

Test problems

- advection 70

- compositional step-function 351

- detonation 434, 435, 437

- grid adaption 297

- radially-symmetric 103

- reaction network 432

- shock tube see Shock tube problem

- simple combustion 90

- simple EOS 432

- Sod's 2,16, 63

- SPH 464,471-473

- step test 352

Theorem of Harten 73

Thermalization length 221

Thermodynamics, first law of 197, 213

Thermonuclear reactions 212

Three body reaction 419, 420

Three plus one formalism 143-145

TIGA 393

Time centering 279

Time scale

- acoustic vs dynamic 130

- burning 407

- convective 408

- core bounce 367

- core collapse 360, 387

- diffusion of composition 407

- diffusion of heat 407

- dynamical 264

- dynamical vs collapse 360

- dynamical, SN 251

- electromagnetic 119

- flow 427

- formation of jets 345

- free fall 327

- grid 296

- grid response vs flow 259

- hot bubble convection 381

- human life 345

- ignition 407

- neutrino diffusion 363

- nuclear burning 214, 264, 427

- nuclear vs convective 408

- radiative vs dynamic 13, 95

- reaction 91

- resLction vs flow 50, 90, 91, 414

- RT instability 375, 378

- simulation vs dynamic 348

- simulation vs time step 47

- sound travel 408, 450

- strong and em. reactions 361

- sweep-up 314

- thermal 264

- variety of 47

Time slicing 145,146

- geodesic 146,147

Time step restriction 312, 430

Time step, individual 463

Total variation 68

Transverse derivative terms 104

Triple umbilic point 125,126

Triple-a reaction 432

Turbulent combustion 413-415, 417,

428

TVD region 73

- second-order methods 74

Two body reaction 420

Two-temperature description 228

Two-temperature means 228

Ultra-relativistic limit 273

Umbilic point 126

Unphysical solution

- artificial stellar pulsation 319

- combustion 90

- detonation 92, 93, 410, 434-436

- from discrepancy EOS vs opaxrity

274

- magnetic monopoles 130

- negative pressure 112, 432

- shock speed 55

- species mass fraction 429

- with linearized Riemann solver 112

- with nonrelativistic Riemann solver

139

- with SPH 463, 467

Page 522

508 Index

Unstable, unconditionally 236

Upstream direction 43

Upwind direction 43

Vanishing viscosity approach 23

Vanishing viscosity solution 24

Variable smoothing length 463, 470,

471

Variables

- dependent 199

- independent 199

- physical 199

- primitive, recovery of 139

- state 8

Vectorization 431,474

Vibrational excitation 273

VIRGO 392

Viscosity

- artificial see Artificial viscosity

- bulk 225

- numerical see Artificial viscosity

- pseudo see Pseudo viscosity

- radiation 220, 223, 224

- vanishing see Vanishing viscosity

Viscous

- energy dissipation 285

- length scale 286

- momentum transfer 285

- pressure 238

- pressure tensor 285

- profile 112

- shock layer 1

Vortex shedding 456

Wall heating 115

Wave

- acoustic 100,130, 234, 237

- - spurious 114

- Alfven see Alfven wave

compressional 122

- deflagration 91

- entropy

- - spurious 115

- intermediate 123, 126,127

- magnetosonic 122

- - fast 122,123

- - slow 122,123

- rarefaction see Rarefaction wave

- shear 122,141

- Taylor 410

- thermal 222

- transverse 122

Weak solution 24

Working surface 446, 447, 450

X-rays 170,372,405

Young stellar objects (YSO) 437

Zeldovich-von Neumann-Doering

model (ZND) 411

ZND structure 91-93

Saas-Fee Advanced Course 27

Lecture Notes 1997

Page 2

Springer

Berlin

Heidelberg

New York

Barcelona

Budapest

Hong Kong

London

Milan

Paris

Singapore

Tokyo

Page 261

9. Numerical Radiation Hydrodynamics 247

9.2 Computational Strategy

Overall Procedure. It is costly to try to solve directly the full angle-frequency

dependent radiation equations coupled to the equations of hydrodynamics be-

cause the dimensionality of the problem (= 4) is too large. So we look for

an approximate technique. Notice that we need only the frequency-integrated

moments to solve the material energy and momentum equations. Therefore

we are led to try a splitting technique in which the spacetime evolution is han-

dled separately from the angle-frequency information. The method described

below is called the multifrequency/grey method. For ease of exposition, as-

sume we use explicit hydro.

1. Suppose we are given the solution at t^. Update velocities to t^~^^ (ma-

terial momentum equation), hence radii and densities to t^^^.

2. If we can assume that the spectral distributions (Ejy/E) and {F^IF) are

known at t̂ "*"̂ (actually they are not!) we can calculate K^ and XF?

hence A:E = f^E/f^p and fcp = X F / X R - Similarly, if we can assume that

the angular distribution of the radiation field is known at t'^'^^ (again, it

isn't), then we have /^ and g,,, hence / and q.

3. At t^'^^ we can now solve the gas energy equation plus the radiation en-

ergy and momentum equations. We thus obtain the run of T^-^^^E^^^,

and F^"*" .̂ This system is nonlinear because the material properties de-

pend nonlinearly on T""^^; to solve it we linearize and iterate to conver-

gence, as described earlier.

In principle, these steps give us the solution at the advanced time. But ac-

tually they don't because we assumed we already knew (Ej^/E), {Fty/F),fi^,

and qiy at P"^^, when in fact we didn't. Therefore at the advanced time level

p + i we should now:

4. Solve the angle-dependent transfer equation for the geometric factors

ft, and qj, using source-sink terms evaluated at t^^^ [i.e., KI,{T^'^^),

7?.(r»+i) ].

5. Then solve the frequency-dependent moment equations for the spectral

profiles {Ey/E) and (Fj^/F). Use these updated profiles to update the

information in ACE,XFJ^E, ^uid A:F at t^'^^.

6. Go to step 3, and iterate to consistency.

— In step 4, we ignore frequency-coupling in order to reduce the dimension-

ality of the problem to (r,/x) at a given (i^^t). From this step we derive the

geometric factors /^ and ^^ which depend mainly on angular information

about the radiation field.

— In step 5, we dispose of angular information by using moments, and attempt

to improve the frequency-dependent information contained in {Ej,/E) and

{Fu/F) by accounting for the (i/,r) coupling at a given ,̂ assuming that

the geometric factors {fi,,qi,) are known.

Page 262

248 D. Mihalas

- This procedure works fairly well because fy is only a ratio which is mainly

geometry dependent, whereas {E^jE) and (Fiy/F) are only spectral pro-

files, which are insensitive to geometry.

— Very few calculations have been done to this level of consistency. Most ei-

ther ignore the frequency-coupling in the moments, and lag the Eddington

factors, or they handle the frequency-couphng but assume /»/ = | (multi-

group diffusion).

Spatial Differencing and Unresolved Fronts. It is straightforward to write

difference formulae for the equations written above. We won't produce them

here; the reader can easily work them out. But consider one very difficult

problem: how can we calculate a good value for XF at an interface? A seem-

ingly reasonable value giving the correct optical depth increment between cell

centers is

QJd ^ " ^ d - l / 2 + ^ ^ d - f 1/2

But trouble arises when we fail to resolve a front in the flow (a shock, ion-

ization front, ...). The problem is that the opacity can vary as T" with

a ~ 1 2 — 1 4 i n regions where H and/or He get excited and/or ionize. There-

fore even small changes in T between zones imply a huge change in x- This

rapid variation in x makes it hard to choose the correct value at the inter-

face, and the errors made can clobber the energy transport at the interface.

Generally the result is disasterous. How can we fix this problem? R. Christy

[6] found that we can get "reasonable" results using an energy-weighted har-

monic mean:

(^''-/•^)^(0.-./. + (^^-/^)^(f). Q\ _ V-^/d-l/2 VX/d+i/2

This recipe has been very popular, but it is only ad hoc. The real solution

to the problem is to use an adaptive grid. The problem can't be solved by

brute force with tiny Lagrangian zones. For example, the ionization zone in

Cepheids and RR Lyrae envelopes is about 0.001 to 0.01 scale heights thick,

and moves back and forth over 3-5 scale heights during a pulsation cycle.

A prohibitive number of zones (and timesteps!) would be required to resolve

the front. In contrast, an adaptive grid moves with the front and resolves it.

This approach is the preferred technique.

Method of Solution. For the explicit hydro scheme sketched above, at t^^^ we

get equations coupling j'^^+i^ ^^+i ^ ^^^+1 at d = 1, ... , D. Linearize the equa-

tions around the current estimate of the solution. We get a block tridiagonal

system which can be solved by Gaussian elimination, and iterated to con-

vergence (Newton-Raphson procedure). In Unearizing opacities we can write

/̂ E = k^Kp and XF = ^ F X R ^^^ assume {k^.k^) remain fixed. Or (better)

Page 521

Index 507

- SNI 354,356,357

- SN la 356, 358, 359, 408, 412-415,

417,418

- SNIb 354-360,393

- SNIc 354-360,393

- SN II 354, 358-360, 376, 393

- SN II-L 354, 357, 360

- SNII-P 354,357,360

- SNIIb 357

- SNIIn 354

- SNIIP 376

- spectrum 356

- Type II 251

- Type I 353

- Type II 353,372

Synchrotron radiation transfer 462

TAMA 392

Tangent space 185,188

Test problems

- advection 70

- compositional step-function 351

- detonation 434, 435, 437

- grid adaption 297

- radially-symmetric 103

- reaction network 432

- shock tube see Shock tube problem

- simple combustion 90

- simple EOS 432

- Sod's 2,16, 63

- SPH 464,471-473

- step test 352

Theorem of Harten 73

Thermalization length 221

Thermodynamics, first law of 197, 213

Thermonuclear reactions 212

Three body reaction 419, 420

Three plus one formalism 143-145

TIGA 393

Time centering 279

Time scale

- acoustic vs dynamic 130

- burning 407

- convective 408

- core bounce 367

- core collapse 360, 387

- diffusion of composition 407

- diffusion of heat 407

- dynamical 264

- dynamical vs collapse 360

- dynamical, SN 251

- electromagnetic 119

- flow 427

- formation of jets 345

- free fall 327

- grid 296

- grid response vs flow 259

- hot bubble convection 381

- human life 345

- ignition 407

- neutrino diffusion 363

- nuclear burning 214, 264, 427

- nuclear vs convective 408

- radiative vs dynamic 13, 95

- reaction 91

- resLction vs flow 50, 90, 91, 414

- RT instability 375, 378

- simulation vs dynamic 348

- simulation vs time step 47

- sound travel 408, 450

- strong and em. reactions 361

- sweep-up 314

- thermal 264

- variety of 47

Time slicing 145,146

- geodesic 146,147

Time step restriction 312, 430

Time step, individual 463

Total variation 68

Transverse derivative terms 104

Triple umbilic point 125,126

Triple-a reaction 432

Turbulent combustion 413-415, 417,

428

TVD region 73

- second-order methods 74

Two body reaction 420

Two-temperature description 228

Two-temperature means 228

Ultra-relativistic limit 273

Umbilic point 126

Unphysical solution

- artificial stellar pulsation 319

- combustion 90

- detonation 92, 93, 410, 434-436

- from discrepancy EOS vs opaxrity

274

- magnetic monopoles 130

- negative pressure 112, 432

- shock speed 55

- species mass fraction 429

- with linearized Riemann solver 112

- with nonrelativistic Riemann solver

139

- with SPH 463, 467

Page 522

508 Index

Unstable, unconditionally 236

Upstream direction 43

Upwind direction 43

Vanishing viscosity approach 23

Vanishing viscosity solution 24

Variable smoothing length 463, 470,

471

Variables

- dependent 199

- independent 199

- physical 199

- primitive, recovery of 139

- state 8

Vectorization 431,474

Vibrational excitation 273

VIRGO 392

Viscosity

- artificial see Artificial viscosity

- bulk 225

- numerical see Artificial viscosity

- pseudo see Pseudo viscosity

- radiation 220, 223, 224

- vanishing see Vanishing viscosity

Viscous

- energy dissipation 285

- length scale 286

- momentum transfer 285

- pressure 238

- pressure tensor 285

- profile 112

- shock layer 1

Vortex shedding 456

Wall heating 115

Wave

- acoustic 100,130, 234, 237

- - spurious 114

- Alfven see Alfven wave

compressional 122

- deflagration 91

- entropy

- - spurious 115

- intermediate 123, 126,127

- magnetosonic 122

- - fast 122,123

- - slow 122,123

- rarefaction see Rarefaction wave

- shear 122,141

- Taylor 410

- thermal 222

- transverse 122

Weak solution 24

Working surface 446, 447, 450

X-rays 170,372,405

Young stellar objects (YSO) 437

Zeldovich-von Neumann-Doering

model (ZND) 411

ZND structure 91-93