| 224 | | |
| 225 | | == 1D Radiative Shock Simulations == |
| 226 | | These simulations follow Ch. 4 of Delamarter '01 and his treatment of the 1D steady radiative shock problem. The purpose of this problem is to check that the cooling source terms are being handled correctly. Here are some useful links: |
| 227 | | |
| 228 | | [http://www.pas.rochester.edu/~bearclaw/delamarterplots.html AstroBEAR 1.0 test page] |
| 229 | | |
| 230 | | [http://www.pas.rochester.edu/~bearclaw/paper_library/Delamarter01.pdf Delamarter '01 pdf] |
| 231 | | |
| 232 | | === Equations === |
| 233 | | The shock jump equations for a stationary shock are used to solve for the initial post-shock values. The post-shock velocity v,,2,, can be written as: |
| 234 | | |
| 235 | | [[latex($v_2 = v_1 \frac{(\gamma - 1)M^2 + 2}{(\gamma + 1)M^2} $)]] |
| 236 | | |
| 237 | | where v,,1,, is the ambient velocity, and M is the ambient mach number. Remember that the mach number M = v,,1,,/c where c is the ambient sound speed, and [[latex($c = \sqrt{\frac{\gamma P_{1}}{\rho_{1}}}$)]] where P,,1,, is ambient pressure, and [[latex($\rho_{1}$)]] is ambient density. The post-shock density and pressure ([[latex($\rho_{2}$)]] and P,,2,, respectively) can be found by using mass flux and momentum flux conservation across the shock: |
| 238 | | |
| 239 | | [[latex($\rho_1 v_1 = \rho_2 v_2$)]] |
| 240 | | |
| 241 | | [[latex($\rho_1 v_1 ^2 + P_1 = \rho_2 v_2 ^2 + P_2$)]] |
| 242 | | |
| 243 | | These post-shock values become the boundary conditions for the fluid equations in the cooling region: |
| 244 | | |
| 245 | | [[latex($\rho v = \rho_2 v_2$)]] |
| 246 | | |
| 247 | | [[latex($\rho v ^2 + P = \rho_2 v_2 ^2 + P_2$)]] |
| 248 | | |
| 249 | | [[latex($\frac{d}{dx}(\frac{1}{2}\rho v ^3 + \frac{\gamma}{\gamma - 1} P v) = -\Lambda$)]] |
| 250 | | |
| 251 | | where [[latex($\rho$)]], v, and P represent the density, velocity, and pressure in the cooling region as functions of x, and [[latex($\Lambda$)]] is the cooling rate. |
| 252 | | |
| 253 | | [[BR]] |
| 254 | | === Initial Parameters === |
| 255 | | n,,1,, = 60 particles/cc |
| 256 | | |
| 257 | | v,,1,, = 10^7^ cm/s |
| 258 | | |
| 259 | | T,,1,, = 10^4^ K |
| 260 | | |
| 261 | | For analytic cooling, |
| 262 | | |
| 263 | | [[latex($\Lambda = n^{2} \alpha T^{\beta}$)]] |
| 264 | | |
| 265 | | where n and T and the number density and temperature in the post-shock region respectively. |
| 266 | | |
| 267 | | For these simulations, |
| 268 | | |
| 269 | | [[latex($\beta$)]] = 2 |
| 270 | | |
| 271 | | [[latex($\alpha$)]] = 1.23786 * 10^-34^ erg*cm^3^/s/K^2^ |
| 272 | | |
| 273 | | cell length = 2.5 x 10^15^ cm |
| 274 | | |
| 275 | | problem domain = 400 cells |
| 276 | | |
| 277 | | final time ~ 4000 years |
| 278 | | |
| 279 | | [[BR]] |
| 280 | | [[CollapsibleStart(Results)]] |
| 281 | | === Results === |
| 282 | | A run with no cooling was done to make sure the post-shock values were correct. With no cooling (adiabatic), the hydrodynamic quantities should jump discontinuously at the shock, and then remain constant while the shock itself remains stationary. Only the first 100 cells of the simulation are shown here. |
| 283 | | |
| 284 | | || Adiabatic || |
| 285 | | || [[Image(ad_press0000.png, width=400)]] || [[Image(ad_temp0000.png, width=400)]] || |
| 286 | | ||= [attachment:ad_press.gif Pressure Movie] =||= [attachment:ad_temp.gif Temperature Movie] |
| 287 | | || [[Image(ad_rho0000.png, width=400)]] || [[Image(ad_vel0000.png, width=400)]] || |
| 288 | | ||= [attachment:ad_rho.gif Density Movie] =||= [attachment:ad_vel.gif Velocity Movie] =|| |
| 289 | | |
| 290 | | Then, a run with the aforementioned analytic cooling parameters was done. The hydrodynamic quantities will change discontinuously at the shock as before, but then change continuously according to the cooling rate. The ambient temperature is set to be the floor temperature, so that no cooling occurs at or below this temperature. Once the post-shock temperature reaches this floor temperature, the hydrodynamic quantities will remain constant. This is called the quiescent region. Only the first 100 cells of the simulation are shown here. |
| 291 | | |
| 292 | | || Analytic Cooling || |
| 293 | | || [[Image(an_press0000.png, width=400)]] || [[Image(an_temp0000.png, width=400)]] || |
| 294 | | ||= [attachment:an_press.gif Pressure Movie] =||= [attachment:an_temp.gif Temperature Movie] |
| 295 | | || [[Image(an_rho0000.png, width=400)]] || [[Image(an_vel0000.png, width=400)]] || |
| 296 | | ||= [attachment:an_rho.gif Density Movie] =||= [attachment:an_vel.gif Velocity Movie] =|| |
| 297 | | |
| 298 | | For Dalgarno-McCray (DM) cooling, some of the initial parameters had to be changed. The initial velocity used in the previous simulations would give rise to a post-shock temperature in the unstable region of the DM curve. This led to the radiative instability, and thus unsteady shocks. To avoid this, the initial velocity was lowered to 80 km/s. Also, the cooling rate from the DM curve is much higher than that of the analytic form. In order to see the cooling region, the grid had to be shrunk, resulting in a cell length of 10^13^ cm. |
| 299 | | |
| 300 | | || DM Cooling || |
| 301 | | || [[Image(DM_temp0010.png, width=400)]] || |
| 302 | | ||= [attachment:DM_temp.gif Temperature Movie] =|| |
| 303 | | |
| 304 | | For Non-equilibrium (NEQ) cooling, the length scale had to be adjusted. The shock velocity is the same as the DM case, but the cell length is now 7.75e11 cm. |
| 305 | | |
| 306 | | || NEQ Cooling || |
| 307 | | || [[Image(NEQ_temp0010.png, width=400)]] || |
| 308 | | ||= [attachment:NEQ_temp.gif Temperature Movie] =|| |
| 309 | | |
| 310 | | Z Cooling has the same initial parameters and scaling as the NEQ case. |
| 311 | | |
| 312 | | || Z Cooling || |
| 313 | | || [[Image(Z_temp0010.png, width=400)]] || |
| 314 | | ||= [attachment:Z_temp.gif Temperature Movie] =|| |
| 315 | | |
| 316 | | Here is a fun plot that shows the temperature profiles for DM, NEQ, and Z cooling. This gives you a good idea of the temperature ranges that each type of cooling works in. Also, you can see the differences in the cooling length which is why I would normally change the length scale depending on the type of cooling. |
| 317 | | |
| 318 | | [[Image(allcool_temp.png,width=400)]] |
| 319 | | |
| 320 | | [[CollapsibleEnd]] |
| 321 | | |
| 322 | | [[BR]] |
| 323 | | === Including MHD === |
| 324 | | Making these simulations work with MHD was no simple task. Having a magnetic field adds another jump condition, and it makes the equations within the cooling region more complicated. Here are what the equations look like now: |
