Changes between Version 29 and Version 30 of AstroBearProjects/resistiveMHD
- Timestamp:
- 05/27/13 11:06:55 (11 years ago)
Legend:
- Unmodified
- Added
- Removed
- Modified
-
AstroBearProjects/resistiveMHD
v29 v30 2 2 3 3 The resistive equation can be written as:[[BR]] 4 [[latex( \frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B}))]]4 [[latex($\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})$)]] 5 5 The magnetic diffusivity eta is a function of temperature as we all know. [[BR]] 6 6 Now for the first equation, what if we want to include the temperature dependence into the code? From the next section, you can see that when a [[BR]] 7 7 field configuration in equilibrium is subject to strong diffusion, usually heating would occur and surppress the local resistivity and thus the [[BR]] 8 8 diffusivity. By expanding the first equation, we have the form:[[BR]] 9 [[latex( \frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B} + \nabla \eta \times (\nabla \times \textbf{B}))]]9 [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B} + \nabla \eta \times (\nabla \times \textbf{B})$)]] 10 10 Now the second term depends on all of the three components of B. So we end up with equations in which the time variance of Bx, By and Bz depend [[BR]] 11 11 on each other.[[BR]] … … 15 15 with position but ignore its own spacial variation. This would give us a form of resistive MHD similar to that of the thermal diffusion case but [[BR]] 16 16 with temperature dependence built in. [[BR]] 17 [[latex( \frac{\partial \textbf{B}}{\partial t}=\eta (T) \nabla^2 \textbf{B})]]17 [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta (T) \nabla^2 \textbf{B}$)]] 18 18 Unfortunately, this does not work because the diffusion equation itself has to be divergence free. Our roughest approximation, treating the resistivity [[BR]] 19 19 as a constant satisfies the requirement as long as the divergence and the Laplacian are commutable. If we try to do the same thing to the above equation, [[BR]] 20 20 we end up getting:[[BR]] 21 [[latex( \frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B})]]21 [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B}$)]] 22 22 The first term on the right hand side is zero but the second term is not, especially at sharp temperature fronts where grad(T) is large.[[BR]] 23 23 So this approximation only works under '''slowly varying temperature''' situation.[[BR]] … … 25 25 total energy during the evolution, especially for the force-free field case. This phenomenon is explained below. [[BR]] 26 26 In AstoBEAR, we explicitly calculate the resistivity induced current on the cell edges, following equation: [[BR]] 27 [[latex( \textbf{J} = \eta \nabla \times \textbf{B})]] [[BR]]27 [[latex($\textbf{J} = \eta \nabla \times \textbf{B}$)]] [[BR]] 28 28 The stencil for this explicit solver is a 3 by 3 cube surrounding the cell we want to update. The magnetic field is represented by the aux field, which is [[BR]] 29 29 centered on the cell faces. Its curl therefore reside on the cell edges. Here is an example on calculating the diffusive current on the x direction, notice [[BR]] 30 30 that the red arrow is where we are calculating the diffusive current, the green arrows are where the magnetic field originally resides: [[BR]] 31 [[Image(resistive_diagram.png, 20% )]] [[Image(http://www.pas.rochester.edu/~shuleli/res_dia1.png, 30%)]] [[BR]]31 [[Image(resistive_diagram.png, 20%$)]] [[Image(http://www.pas.rochester.edu/~shuleli/res_dia1.png, 30%$)]] [[BR]] 32 32 The actual code looks like the following in 2D (jy are initialized to 1): [[BR]] 33 33 {{{ … … 41 41 }}} 42 42 After calculating and storing jx, jy and jz of the given grid, the change of magnetic field due to diffusion can be calculated as： [[BR]] 43 [[latex( \frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J})]] [[BR]]43 [[latex($\frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J}$)]] [[BR]] 44 44 In the case of AMR, we need to store the emf using function storefixupfluxes. In case of resistivity, it is simply the diffusive current. In 2D case, the emfs [[BR]] 45 45 are only in z. So we need to use f instead. The code looks like: [[BR]] … … 64 64 the fluid is everywhere zero, and no heat is generated by the current. [[BR]] 65 65 If we dot the resistive induction equation with the magnetic field B, we obtain the time evolution equation for magnetic energy:[[BR]][[BR]] 66 [[latex( \partial (B^2)/\partial t + \nabla \cdot \textbf{S} = - \textbf{J}\cdot\textbf{E} = -j^2/\eta)]][[BR]][[BR]]67 where [[latex( \textbf{S}=\textbf{J} \times \textbf{B})]] is the magnetic energy flux caused by resistive diffusion and [[latex(j=|\textbf{J}|)]] is the magnitude [[BR]]66 [[latex($\partial (B^2)/\partial t + \nabla \cdot \textbf{S} = - \textbf{J}\cdot\textbf{E} = -j^2/\eta$)]][[BR]][[BR]] 67 where [[latex($\textbf{S}=\textbf{J} \times \textbf{B}$)]] is the magnetic energy flux caused by resistive diffusion and [[latex($j=|\textbf{J}|$)]] is the magnitude [[BR]] 68 68 of the diffusive current.[[BR]] 69 69 In this equation, the S term accounts for the redistribution of magnetic energy (and thus the redistribution of total energy), and the last term accounts for [[BR]] 70 70 the loss of magnetic energy due to reconnection.[[BR]] 71 71 The total energy change for the resistive step is therefore:[[BR]][[BR]] 72 [[latex( \partial \epsilon/\partial t + \nabla \cdot \textbf{S} = 0)]][[BR]][[BR]]73 Here the [[latex( j^2/\eta)]] dissipation term is absent because the dissipation of magnetic energy does not change the total energy: the loss of magnetic energy [[BR]]72 [[latex($\partial \epsilon/\partial t + \nabla \cdot \textbf{S} = 0$)]][[BR]][[BR]] 73 Here the [[latex($j^2/\eta$)]] dissipation term is absent because the dissipation of magnetic energy does not change the total energy: the loss of magnetic energy [[BR]] 74 74 is converted into thermal energy. [[BR]] 75 75 In the code, the energy flux as a result of magnetic diffusion needs to be calculated explicitly using: [[BR]] 76 [[latex( \textbf{S}=\textbf{J} \times \textbf{B})]] [[BR]]76 [[latex($\textbf{S}=\textbf{J} \times \textbf{B}$)]] [[BR]] 77 77 The energy fluxes reside on the cell faces while the diffusive current reside on the cell edges. We therefore need to compute a face average of the diffusive [[BR]] 78 78 current as well as the magnetic field components which are not normal to the face using the surrounding edges. In the last diagram, the blue arrows connected by [[BR]] 79 79 dashed lines are what used to compute the energy flux. [[BR]] 80 80 The magnetic field can be updated from the diffusive currents by: [[BR]] 81 [[latex( \frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J})]] [[BR]]81 [[latex($\frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J}$)]] [[BR]] 82 82 The updating of magnetic field and the calculation of the energy fluxes can be done at the same time, as in the code: [[BR]] 83 83 {{{ … … 95 95 Finally, the energy is updated using the divergence of the energy fluxes. This finishes the diffusive process.[[BR]] 96 96 The resistive time scale is explicitly calculated by: [[BR]] 97 [[latex( dt_{resistive} = \frac{c dx}{\eta})]] [[BR]]97 [[latex($dt_{resistive} = \frac{c dx}{\eta}$)]] [[BR]] 98 98 where c is the cfl number (set to be 0.5 in the code).[[BR]][[BR]] 99 99 … … 104 104 105 105 To investigate the Sweet-Parker problem, we introduce the following magnetic field configuration called "sheer pinch": [[BR]] 106 [[latex( B_y(x) = b_0 tanh(x/a))]][[BR]][[BR]]106 [[latex($B_y(x) = b_0 tanh(x/a)$)]][[BR]][[BR]] 107 107 108 where in computational units we choose [[latex( b_0 = 1)]] and [[latex(a = 0.5)]].[[BR]]108 where in computational units we choose [[latex($b_0 = 1$)]] and [[latex($a = 0.5$)]].[[BR]] 109 109 The density profile is chosen so that the pressure equilibrium can be maintained with constant temperature.[[BR]][[BR]] 110 [[latex( \rho(x) = \rho_0/cosh^2(x/a)+\rho_c)]][[BR]][[BR]]111 where [[latex( \rho_0 = 1)]] and [[latex(a = 0.5)]]. The temperature is set to be constant as 0.5. [[BR]]110 [[latex($\rho(x) = \rho_0/cosh^2(x/a)+\rho_c$)]][[BR]][[BR]] 111 where [[latex($\rho_0 = 1$)]] and [[latex($a = 0.5$)]]. The temperature is set to be constant as 0.5. [[BR]] 112 112 The domain is set to be -6.4 < x < 6.4 and -12.8< y < 12.8, with fixed resolution 480 * 960. The boundaries are all open. The initial profile is plotted below: [[BR]] 113 [[Image(http://www.pas.rochester.edu/~shuleli/hhc_plot.png, 30% )]][[BR]][[BR]]113 [[Image(http://www.pas.rochester.edu/~shuleli/hhc_plot.png, 30%$)]][[BR]][[BR]] 114 114 The initial state is in pressure equilibrium though unstable. There are two ways to generate instabilities. The first way is to artificially increase the resistivity [[BR]] 115 115 at the center of the domain. This increase will result in a higher reconnectivity, which will eventually bend magnetic field. This creates an X point where field [[BR]] … … 117 117 to the direction of the sheer pinch. The box surrounding the X point where the outflows (Petschek shock) come out of is called the "Sweet-Parker Box". The following [[BR]] 118 118 diagrams show how increased resistivity at the center of a sheer pinch drives Petschek shock.[[BR]] 119 [[Image(http://www.pas.rochester.edu/~shuleli/multi_test/resistive_instab.png,30% )]] [[Image(http://www.pas.rochester.edu/~shuleli/hhc_flowpattern.png, 40%)]] [[BR]][[BR]]119 [[Image(http://www.pas.rochester.edu/~shuleli/multi_test/resistive_instab.png,30%$)]] [[Image(http://www.pas.rochester.edu/~shuleli/hhc_flowpattern.png, 40%$)]] [[BR]][[BR]] 120 120 The Mach number in 2D case is plotted in pseudo-color: [[BR]] 121 [[Image(http://www.pas.rochester.edu/~shuleli/hhc_mach.png, 30% )]][[BR]][[BR]]121 [[Image(http://www.pas.rochester.edu/~shuleli/hhc_mach.png, 30%$)]][[BR]][[BR]] 122 122 123 123 To test Sweet-Parker problem in AstroBEAR, we construct the sheer pinch using the above setup, and modify . [[BR]][[BR]] … … 125 125 The following figure shows the Petschek shock from a reconnection spot at the center. Colored variable is the kinetic energy in log scale, magnetic field is illustrated [[BR]] 126 126 by white lines. [[BR]] 127 [[Image(http://www.pas.rochester.edu/~shuleli/multi_test/sp1_0020.png, 30% )]]127 [[Image(http://www.pas.rochester.edu/~shuleli/multi_test/sp1_0020.png, 30%$)]] 128 128 [[BR]] 129 129 A movie with 2-level AMR:[[BR]] … … 131 131 132 132 The second way to perturb the sheer pinch is to put in a sinusoidal perturabtion on magnetic field: [[BR]][[BR]] 133 [[latex( B_x = B_p sin(kx)cos(ky))]][[BR]]134 [[latex( B_y = -B_p cos(kx)sin(ky))]][[BR]][[BR]]133 [[latex($B_x = B_p sin(kx)cos(ky)$)]][[BR]] 134 [[latex($B_y = -B_p cos(kx)sin(ky)$)]][[BR]][[BR]] 135 135 where Bp is the perturbation amplitude and the wave number:[[BR]] 136 [[latex( k = 2\pi/L_y)]].[[BR]]136 [[latex($k = 2\pi/L_y$)]].[[BR]] 137 137 This perturbation creates periodical X points and O points, which leads to bright and dark spots of magnetic pressure. Material flow into the dark spots as a result [[BR]] 138 138 of pressure imbalance, which creates periodical dense "islands".The growth rate and the size of the "islands" depend on resistivity and the strength of the perturbation.[[BR]][[BR]] 139 139 140 [[Image(http://www.pas.rochester.edu/~shuleli/multi_test/mihr_0180.png,40% )]][[BR]][[BR]]140 [[Image(http://www.pas.rochester.edu/~shuleli/multi_test/mihr_0180.png,40%$)]][[BR]][[BR]] 141 141 142 142 To watch a full movie, [http://www.pas.rochester.edu/~shuleli/multi_test/mihr.gif click here][[BR]][[BR]] … … 149 149 In cgs units, the real magnetic diffusivity can be written as:[[BR]] 150 150 151 [[latex( \eta = c^2/4\pi \sigma)]]151 [[latex($\eta = c^2/4\pi \sigma$)]] 152 152 153 153 where sigma is the material conductivity which depends on the electron-ion collision rate mu which in term depends on the electron temperature to the three halves. [[BR]] 154 154 The cross field diffusivity can thus be expressed as function of electron temperature if we assume the plasma is fully ionized:[[BR]] 155 155 156 [[latex( \eta = 8.2243 \times 10^5 Z_{eff} F(Z_{eff}) \ln \Lambda T_{e}^{-3/2})]]156 [[latex($\eta = 8.2243 \times 10^5 Z_{eff} F(Z_{eff}) \ln \Lambda T_{e}^{-3/2}$)]] 157 157 158 158 where T is the electron temperature in eV, Zeff is the effective ion charge.The Coulomb logarithm for interested number density is plotted below: [[BR]] 159 [[Image(http://www.pas.rochester.edu/~shuleli/CoulombLog.png, 30% )]][[BR]]159 [[Image(http://www.pas.rochester.edu/~shuleli/CoulombLog.png, 30%$)]][[BR]] 160 160 [[BR]][[BR]] 161 161 162 162 Here, T is the same as before, n is the electron number density in the units of per cubic centimeter. F is a function of effective ion charge:[[BR]][[BR]] 163 163 164 [[latex( F(Z) = \frac{1+1.198Z+0.222Z^2}{1+2.996Z+0.753Z^2})]][[BR]][[BR]]164 [[latex($F(Z) = \frac{1+1.198Z+0.222Z^2}{1+2.996Z+0.753Z^2}$)]][[BR]][[BR]] 165 165 166 166 Since the electron's maxwellian will be distorted when there is a strong outer magnetic field applied, the realistic resistivity can depend on the local field orientation. [[BR]] 167 167 The anisotropicity is around 2 for the parallel and cross field resistivity. This requires the electron mean free path to be much longer than the gyroradius under such a field: [[BR]] 168 [[latex( R_B << \lambda_{mfp})]][[BR]][[BR]]168 [[latex($R_B << \lambda_{mfp}$)]][[BR]][[BR]] 169 169 or, the gyro frequency to be much greater than the mean electron ion collision frequency. [[BR]] 170 170 This effect is not considered currently. [[BR]][[BR]]
