| | 54 | '''Implicit Implementation of Resistivity''' [[BR]][[BR]] |
| | 55 | |
| | 56 | The resistive equation can be written as:[[BR]] |
| | 57 | |
| | 58 | [[latex($\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})$)]] |
| | 59 | |
| | 60 | The magnetic diffusivity eta is a function of temperature as we all know. [[BR]] |
| | 61 | |
| | 62 | 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 field configuration in equilibrium is subject to strong diffusion, usually heating would occur and surppress the local resistivity and thus the diffusivity. By expanding the first equation, we have the form:[[BR]] |
| | 63 | |
| | 64 | [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B} + \nabla \eta \times (\nabla \times \textbf{B})$)]] |
| | 65 | |
| | 66 | 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 on each other.[[BR]] |
| | 67 | |
| | 68 | So the coefficient array is thus no longer a tri-diagonal matrix. For cases where resistive speed is slow (which is usually the case), we can use explicit solver to treat the problem instead.[[BR]] |
| | 69 | |
| | 70 | One may wonder if it is possible to throw away the second term on the right hand side of the diffusion equation to just let the diffusivity to vary 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 with temperature dependence built in. [[BR]] |
| | 71 | |
| | 72 | [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta (T) \nabla^2 \textbf{B}$)]] |
| | 73 | |
| | 74 | Unfortunately, this does not work because the diffusion equation itself has to be divergence free. Our roughest approximation, treating the resistivity 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, we end up getting:[[BR]] |
| | 75 | |
| | 76 | [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B}$)]] |
| | 77 | |
| | 78 | 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]] |
| | 79 | So this approximation only works under '''slowly varying temperature''' situation.[[BR]] |
| | 80 | |
| | 81 | Upon a closer look at the simulations conducted below, we can observe that the energy inside the domain does not conserve. There is a small increase of total energy during the evolution, especially for the force-free field case. This phenomenon is explained below. [[BR]] |
| | 82 | In the case of resistive MHD, the energy can be dissipated in the form of Joule heat, comparing to the infinite conductivity case, where the voltage inside the fluid is everywhere zero, and no heat is generated by the current. [[BR]] |
| | 83 | If we dot the resistive induction equation with the magnetic field B, we obtain the time evolution equation for magnetic energy:[[BR]][[BR]] |
| | 84 | [[latex($\partial (B^2)/\partial t + \nabla \cdot S = - \textbf{J}\cdot\textbf{E}$)]][[BR]][[BR]] |
| | 85 | Since we know that:[[BR]][[BR]] |
| | 86 | [[latex($\textbf{J}=\eta\textbf{E}$)]][[BR]] |
| | 87 | and also:[[BR]][[BR]] |
| | 88 | [[latex($\textbf{E}=\nabla \times \textbf{B}$)]][[BR]][[BR]] |
| | 89 | This leads to an extra energy source term:[[BR]][[BR]] |
| | 90 | [[latex($\rho\frac{\partial e}{\partial t}=\eta (\nabla \times \textbf{B})^2$)]][[BR]][[BR]] |
| | 91 | This term can be large in the fast-varying field region, and can be a major mechanism of the dissipation of energy released by the magnetic reconnection. [[BR]][[BR]][[BR]] |
| | 92 | |
| | 93 | Here I overview two approaches to solve the resistive MHD equations. [[BR]][[BR]] |
| | 94 | (1)Solving the decomposed format with divergence cleaning (Implicit) [[BR]] |
| | 95 | This is the simplest method in terms of solving scheme. This method is based on the simplification that in some cases eta is a constant in the region of interest. This will simplify the above equation into the following form:[[BR]] |
| | 96 | |
| | 97 | [[latex($\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B}$)]] |
| | 98 | |
| | 99 | This equation holds true for each separate field component, and the form resembles the linear isotropic thermal conduction equation. We can solve the elliptic equation for three times or we can write a expanded linear system so that the solver is only called once: [[BR]] |
| | 100 | |
| | 101 | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_matrix1.png, 20%)]][[BR]] |
| | 102 | |
| | 103 | In the constant diffusivity case, the matrix M are identical for the Bx, By and Bz components. Define C=[-6k,k,k,k,k,k,k], where [[latex($k = \eta dt/dx^2$)]]. One can immediately recognize that k is the diffusion CFL number. Then the expanded matrix is still "tridiagonal" in a 3-D sense and can be blocked into three identical parts: [[latex($M_i=C-[1,0,0,0,0,0,0]$)]] on each row. The RHS vector are source terms that can be easily obtained by the following the CR scheme:[[BR]] |
| | 104 | [[BR]] |
| | 105 | [[latex($S_{xi}=(-C-[1,0,0,0,0,0,0]) \cdot B_{xi}$)]][[BR]] |
| | 106 | [[BR]] |
| | 107 | |
| | 108 | The matrix equation can be solved using linear solvers at hand in parallel.[[BR]] |
| | 109 | The field obtained above resides at cell centers. We can do interpolations to retrieve the new magnetic field at each cell faces. But since we are treating the x, y and z fields independently and applying sources that is not necessarily divergence free to them, we are running the risk of introducing divergence to the face centered magnetic field components which should be divergence free. Experiments show that this decomposing method can introduce large divergence to the face centered field. Thus a divergence cleaning scheme is required once we finish the diffusion update. Here we introduce the divergence cleaning scheme implemented in AstroBEAR. Considering a cell with face centered field components:[[BR]] |
| | 110 | |
| | 111 | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_dia1.png, 20%)]][[BR]][[BR]] |
| | 112 | |
| | 113 | We can find the divergence source term: [[latex($d=\nabla \cdot B$)]]. It is obvious that if B is face centered, then its divergence should be cell centered.[[BR]] |
| | 114 | Now we solve the following equation system: [[latex($\nabla^2 \phi = d$)]] by using d as a source. This equation system is again well defined, and can be constructed using the same C vector described before with k = 1. Once we obtained phi, we can go back to clean the face centered field: [[BR]] |
| | 115 | |
| | 116 | [[latex($B'=B-\nabla \phi$)]] [[BR]] |
| | 117 | |
| | 118 | We can easily see that the new field is divergence free: |
| | 119 | |
| | 120 | [[latex($\nabla \cdot B' = \nabla \cdot B - \nabla \cdot \nabla \phi = \nabla \cdot B - \nabla^2 \phi = 0 $)]] [[BR]] |
| | 121 | |
| | 122 | Since phi is cell centered, its gradient is face centered which matches the face centered field components. So this method is exact. [[BR]] |
| | 123 | |
| | 124 | The advantages and disadvantages are listed below: [[BR]] |
| | 125 | |
| | 126 | (A1) Simple scheme, works fine with slowly varying temperature situations. [[BR]] |
| | 127 | |
| | 128 | (A2) The divergence cleaning scheme is simple and exact. [[BR]] |
| | 129 | |
| | 130 | (D1) The scheme requires constructing cell centered field components and then reconstruct the face centered field components from the updated cell centered field. The cell centered field components are not divergence free to start with, which is diffused using a decomposed scheme. The process may introduce divergence that is large enough that the cleaning mechanism could not work well enough. [[BR]] |
| | 131 | |
| | 132 | (D2) The divergence cleaning scheme has a flexible boundary condition. Normally in the code, we take the Neumann boundary condition. This condition impose the requirement that the face centered field components at the domain boundaries do not need to be adjusted. The condition in itself may not be optimistic, and is likely to be depending on the specific problem. [[BR]] |
| | 133 | |
| | 134 | (D3) Cannot treat fast varying temperature situation where the temperature length scale is comparable or smaller comparing to the field length scale. [[BR]] |
| | 135 | |
| | 136 | (2)Solving the coupled format with varying temperature (Explicit) [[BR]] |
| | 137 | |
| | 138 | Here we consider the full equation:[[BR]] |
| | 139 | |
| | 140 | [[latex($\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})$)]] |
| | 141 | |
| | 142 | Here we treat the face centered field components directly by looking at the following diagram: [[BR]] |
| | 143 | |
| | 144 | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_dia2.png, 30%)]][[BR]]]] |
| | 145 | |
| | 146 | We use Bz at the top face of the cell as an example for demonstration. Here, the changing of Bz does not only depend on the Bz alone, but also the field components of Bx and By. The second order accurate scheme requires 12 closest points surrounding the face center. In the diagram, F denotes face centers, X denotes edge centers. The face centers F1 to F8 have distance of [[latex($sqrt{1/2}$)]], and the face centers F9 to F12 have distance of 1. A larger lattice cell can be considered using further face centers when a higher accuracy is required. The correlation vector C in this case has 13 elements which relates the face center we are considering F0 to the 12 adjacent face centers F1 to F12. Its elements can be written as: [[BR]] |
| | 147 | |
| | 148 | [[latex($C_0 = -(k_{X0} + ... + k_{X4})$)]][[BR]] |
| | 149 | [[latex($C_1 = -k_{X0}$)]][[BR]] |
| | 150 | [[latex($C_2 = k_{X0}$)]][[BR]] |
| | 151 | [[latex($C_3 = k_{X2}$)]][[BR]] |
| | 152 | [[latex($C_4 = -k_{X2}$)]][[BR]] |
| | 153 | [[latex($C_5 = k_{X3}$)]][[BR]] |
| | 154 | [[latex($C_6 = -k_{X3}$)]][[BR]] |
| | 155 | [[latex($C_7 = -k_{X4}$)]][[BR]] |
| | 156 | [[latex($C_8 = k_{X4}$)]][[BR]] |
| | 157 | [[latex($C_9 = k_{X1}$)]][[BR]] |
| | 158 | [[latex($C_{10} = k_{X2}$)]][[BR]] |
| | 159 | [[latex($C_{11} = k_{X3}$)]][[BR]] |
| | 160 | [[latex($C_{12} = k_{X4}$)]][[BR]] |
| | 161 | |
| | 162 | where k are the CFL number at the edge centers. Since the CFL are defined at the cell centers initially, we need to interpolate the edge centered CFL using the four nearest cell center temperatures. The relation between temperature and the diffusivity is described in the next section. Remember that because the x, y, z field components now reside on different face centers, the value of C vector for a certain cell will be different for different field components. [[BR]] |
| | 163 | |
| | 164 | The matrix equation has the form of: [[BR]] |
| | 165 | |
| | 166 | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_matrix2.png, 20%)]][[BR]]]] |
| | 167 | |
| | 168 | where M_x, M_y and M_z are the diagonal contributions from the face centers F9 to F12: M = [C_0, C_9, C_10, C_11, C_12, ....]-[1,0,0,0...]; and the blocked M is different for x, y and z components as stated before. There are also nondiagonal blocks which connect Bx block to By and etc. These are marked by D, which comes from the C_1 to C_8 contributions. The source vector is obtained by construct the matrix with -C_0-1; C1; C2; ... C12 and dot it with the field vector: [Bx1 ... Bxn; By1 ... Byn; Bz1 ... Bzn]. The scheme involves looking at each face centers and find out it C vector. Once the C vector is obtained, it can be used to construct the matrix and the source vector. [[BR]] |
| | 169 | This matrix equation is still solvable using linear solver package in parallel, but the construction of the linear system is much more complicated comparing to approach (1).[[BR]] |
| | 170 | |
| | 171 | (A1) Overall more sophisticated than (1), can treat rapid changing temperature. [[BR]] |
| | 172 | (A2) Method can be extended to treating viscosity. [[BR]] |
| | 173 | (A3) Using a curl field as the source thus preventing unwanted field divergence from generating. [[BR]] |
| | 174 | (D1) Computationally much more intense while still requires subcycling. The CFL is instead given by [[latex($\eta dt/l^2$)]][[BR]][[BR]] |
| | 175 | |
| | 176 | '''Microphysical Resistivity as Explicit Source Term''' |
| | 177 | |
