| 1 | | = ======================================= = |
| 2 | | = The 2.5-D Magnetized Cloud with Resistivity = |
| 3 | | = ======================================= = |
| | 1 | '''Resistivity''' [[BR]][[BR]] |
| | 2 | |
| | 3 | The resistive equation can be written as:[[BR]] |
| | 4 | |
| | 5 | [[latex(\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B}))]] |
| | 6 | |
| | 7 | The magnetic diffusivity eta is a function of temperature as we all know. [[BR]] |
| | 8 | |
| | 9 | 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]] |
| | 10 | |
| | 11 | [[latex(\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B} + \nabla \eta \times (\nabla \times \textbf{B}))]] |
| | 12 | |
| | 13 | 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]] |
| | 14 | |
| | 15 | 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]] |
| | 16 | |
| | 17 | 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]] |
| | 18 | |
| | 19 | [[latex(\frac{\partial \textbf{B}}{\partial t}=\eta (T) \nabla^2 \textbf{B})]] |
| | 20 | |
| | 21 | 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]] |
| | 22 | |
| | 23 | [[latex(\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B})]] |
| | 24 | |
| | 25 | 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]] |
| | 26 | So this approximation only works under '''slowly varying temperature''' situation.[[BR]] |
| | 27 | |
| | 28 | 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]] |
| | 29 | |
| | 30 | In AstoBEAR, we explicitly calculate the resistivity induced current on the cell edges, following equation: [[BR]] |
| | 31 | [[latex(\textbf{J} = \eta \nabla \times \textbf{B})]] [[BR]] |
| | 32 | The magnetic field is represented by the aux field, which is 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 that the red arrow is where we are calculating the diffusive current, the green arrows are where the magnetic field originally resides: [[BR]] |
| | 33 | |
| | 34 | << figure. [[BR]] |
| 21 | | [[Image(http://www.pas.rochester.edu/~shuleli/mceq10.png, 6%)]] |
| | 45 | In the case of AMR, we need to store the emf using function storefixupfluxes. In case of resistivity, it is simply the diffusive current. The code looks like: [[BR]] |
| | 46 | |
| | 47 | << code. [[BR]] |
| | 48 | |
| | 49 | 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]] |
| | 50 | If we dot the resistive induction equation with the magnetic field B, we obtain the time evolution equation for magnetic energy:[[BR]][[BR]] |
| | 51 | [[latex(\partial (B^2)/\partial t + \nabla \cdot \textbf{S} = - \textbf{J}\cdot\textbf{E} = -j^2/\eta)]][[BR]][[BR]] |
| | 52 | 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 of the diffusive current.[[BR]][[BR]] |
| | 53 | In this equation, the S term accounts for the redistribution of magnetic energy (and thus the redistribution of total energy), and the [[latex(j^2/\eta)]] term accounts for the loss of magnetic energy due to reconnection.[[BR]] |
| | 54 | The total energy change for the resistive step is therefore:[[BR]][[BR]] |
| | 55 | [[latex(\partial \epsilon/\partial t + \nabla \cdot \textbf{S} = 0)]][[BR]][[BR]] |
| | 56 | 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 is converted into thermal energy. [[BR]] |
| | 57 | In the code, the energy flux as a result of magnetic diffusion needs to be calculated explicitly using: [[BR]] |
| | 58 | |
| | 59 | [[latex(\textbf{S}=\textbf{J} \times \textbf{B})]] [[BR]] |
| | 60 | |
| | 61 | 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 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 dashed lines are what used to compute the energy flux. [[BR]] |
| | 62 | The magnetic field can be updated from the diffusive currents by: [[BR]] |
| | 63 | |
| | 64 | [[latex(\frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J})]] [[BR]] |
| | 65 | |
| | 66 | The updating of magnetic field and the calculation of the energy fluxes can be done at the same time, as in the code: [[BR]] |
| | 67 | |
| | 68 | << code. [[BR]] |
| | 69 | |
| | 70 | We then call storefixupfluxes for the three energy flux components (as xflux, yflux and zflux in the code).[[BR]] |
| | 71 | Finally, the energy is updated using the divergence of the energy fluxes. This finishes the diffusive process.[[BR]] |
| | 72 | [[BR]] |
| | 73 | |
| | 74 | The resistive time scale is explicitly calculated by: [[BR]] |
| | 75 | |
| | 76 | [[latex(dt_{resistive} = \frac{c dx}{\eta})]] [[BR]] |
| | 77 | |
| | 78 | where c is the cfl number (set to be 0.5 in the code). To change the cfl number, the following line of code can be modified: [[BR]] |
| | 79 | |
| | 80 | where 2d0 can be replaced by any value above 1 (for example, putting this number to 5d0 will reduce the cfl to 0.2). |
| | 81 | |
| | 82 | [[BR]][[BR]] |
| | 83 | |
| | 84 | '''Resistivity Test Problems'''[[BR]][[BR]] |
| 44 | | In the next section we examine a more complicated scenario: the force-free configuration.[[BR]][[BR]][[BR]] |
| 45 | | |
| 46 | | = =================================== = |
| 47 | | = The Force Free Field Subject to Resistivity = |
| 48 | | = =================================== = |
| 49 | | |
| 50 | | We still use the 2.5-D model for the force-free field, where the field has a B_z, B_phi component that only depends on Bessel functions of r. The following vector graph shows how the field looks like initially.[[BR]] |
| 51 | | |
| 52 | | [[Image(http://www.pas.rochester.edu/~shuleli/3dfield/mcfieldvector.png, 40%)]] |
| 53 | | |
| 54 | | In the infinite conduction case where the field is frozen in, the force-free field configuration can sit literally forever because the Lorentz force is self-canceling everywhere. Now with the the resistivity, we would expect the field to drift and break the force-free status. We conduct simulation with strong resistivity and high beta to let the field diffusion play a prominent role in the evolution. The setup has a beta of beta = 0.2 and with '''vr > va > cs'''. One could say that this setup resembles figure 1, the major differences are that (1) the field is force-free; (2) there is no cavity initially. The density and temperature is uniform distributed. The following graph shows the field strength, velocity, density, temperature and the B_x component, respectively on each row.[[BR]] |
| 55 | | |
| 56 | | [[Image(http://www.pas.rochester.edu/~shuleli/resff1.png, 80%)]] |
| 57 | | |
| 58 | | Let us look at the field diffusion first. One of the drastic difference between the FF and non-FF case is that in the FF case although the field dissipates and reduce in strength, it does not spread in space very much. One of the reasons is that in the FF configuration, the field strength at the edge are largely composed by the B_x and B_y field because they are first order Bessel comparing to B_z which is zero order. B_z contributes mostly to the central area around r = 0. The B_phi field is reconnecting itself during the diffusion so that it reduces the strength in a way that it does not have to spread over. Let us look at row 4 for B_x evolution. The blue and red lobes in the B_x field correspond to opposite directions they point to. And the field diffusion happens so that the positive lobe would spread into the negative lobe's territory and vice versa. If the positive field and the negative field move into one zone, they inevitably cancel each other which at the macroscopic level, appears as field reconnection. The reconnection rate thus connects to how fast the colored lobes expand which in turn relates to how strong the resistivity is. The B_x evolution on row 4 shows that the positive and negative lobes are diminishing due to the reconnection and the zone with green (denotes zero field) increases. This dissipation happens without significant expansion since it happens at the inside and edge regions at the same rate, which is slightly different from traditional diffusion where edge always diffuse faster because it has a larger gradient. [[BR]] |
| 59 | | |
| 60 | | From the density plot, we see a dense ring is formed at the edge of the total field strength plot, where the pressure gradient is the greatest. And the material flow is the fastest at the ring. [[BR]] |
| 61 | | The material flow is interesting in itself that it reverts its direction during the evolution: during early phase of the evolution, the material is flowing into the strong field region, which is identical to figure 1: the inflow is there to compensate the pressure loss. Then, however, at time t = 0.1, it reverts its direction to point outwards. The material flow in this situation is subtle since it depends on which region is losing pressure faster. In the force-free configuration, the field cancellation caused by the reconnection is happening not just inside the strong field region but also at the edge. The reversed flow may indicate that the region outside the ring is losing pressure faster comparing to the region inside the ring at time t = 0.1. The temperature plot has no surprise that we also get heating where field pressure is dropping.[[BR]] |
| 62 | | |
| 63 | | We next compare the former case to a weak field case: beta = 0.2 and beta = 20. In the beta = 20 situation, the velocity relation is '''vr > cs > va'''. The following plot shows comparison of field distribution, temperature and density respectively. Obviously the high beta case does not have a dense ring and the density does not redistribute much during the evolution. Another difference is that the central heating is not obvious in the high beta case, which is similar to the comparison we did before: the pressure loss in the high beta case is not significant enough to trigger the heating. [[BR]] |
| 64 | | |
| 65 | | [[Image(http://www.pas.rochester.edu/~shuleli/resff2.png, 50%)]] |
| 66 | | |
| 67 | | Finally we compare the first simulation to a case where the sound speed is comparable to the diffusion speed: '''vr = cs > va'''. This is achieved by raising the uniform temperature initially but keep everything else the same. The following graph shows the result. We can see that there is little difference between the hot case and the high beta case in terms of ring formation and central heating: both cases does not have ring formation and significant central heating. But in the hot case, the field is dissipating faster comparing to the original situation (See row 1). [[BR]] |
| 68 | | |
| 69 | | [[Image(http://www.pas.rochester.edu/~shuleli/resff3.png, 50%)]] |
| 70 | | |
| 71 | | [[BR]][[BR]][[BR]][[BR]] |
| 72 | | = =================================== = |
| 73 | | = 3-D Tangled Field Cloud Evolution = |
| 74 | | = =================================== = |
| 75 | | |
| 76 | | In the 3-D tangled field case, the field will exert force on the cloud material because it is not force-free, and at the same time, it will diminish due to magnetic reconnection as a result of field diffusion. [[BR]] |
| 77 | | [[Image(http://www.pas.rochester.edu/~shuleli/clumpfieldlines.png, 50%)]][[BR]][[BR]][[BR]][[BR]] |
| 78 | | |
| 79 | | = =================================== = |
| 80 | | = Viscosity in AstroBEAR = |
| 81 | | = =================================== = |
| 82 | | |
| 83 | | The viscosity term in the Navier-Stokes Equation can be separated as:[[BR]][[BR]] |
| 84 | | [[latex($dv/dt = \nabla \cdot T$)]][[BR]] |
| 85 | | For the applications we are considering, the fluid is isotropic. Under this assumption, the tensor can be expressed as the following combination of two scalar dynamic viscosities [[latex($\mu$)]] and [[latex($\mu'$)]]:[[BR]][[BR]] |
| 86 | | [[latex($T = 2 \mu \nabla v + \mu' \nabla \cdot v I$)]][[BR]][[BR]] |
| 87 | | where I is the unit tensor.[[BR]] |
| 88 | | It can be proven that the stress tensor should have zero trace, which leads to the following relation between [[latex($\mu$)]] and [[latex($\mu'$)]]:[[BR]][[BR]] |
| 89 | | [[latex($\mu' = -2/3\mu$)]][[BR]][[BR]] |
| 90 | | Combining these relations, we can write down the final form of the separated viscous equation:[[BR]][[BR]] |
| 91 | | [[latex($dv/dt=\nabla \cdot (2 \mu \nabla v) - 2/3 \nabla (\mu \nabla \cdot v)$)]][[BR]][[BR]] |
| 92 | | Notice that this equation already counts the effect of all the elements in a viscosity tensor, assuming the fluid itself is isotropic.[[BR]] |
| 93 | | [[BR]][[BR]] |
| 94 | | |
| 95 | | = =================================== = |
| 96 | | = Viscosity Implementation = |
| 97 | | = =================================== = |
| 98 | | The viscous source term is calculated rigorously by first finding the tensor elements at the surrounding face centers of a cell center considered: |
| 99 | | The [[latex($\mu$)]] sources [[latex($2\mu \nabla v$)]] are naturally located at the cell centers, the [[latex($\mu'$)]] sources however needs an averaging to obtain the cell centered velocity divergence.[[BR]] |
| 100 | | Therefore the scheme requires a 20 point stencil in 3-D configuration, which makes the implicit scheme a bit complicated. The following graph shows such a stencil in 2-D setting:[[BR]][[BR]] |
| 101 | | [[Image(http://www.pas.rochester.edu/~shuleli/vis_stencil.png, 35%)]][[BR]][[BR]] |
| 102 | | Using this stencil setting, one can get face centered velocity divergence for the cell considered, and then calculate the gradient components located at the cell center.[[BR]] |
| 103 | | In the viscous plasma context, the physical viscosity can be calculated rigorously. The plasma viscosity is proportional to the product of ion mean free path and the ion thermal energy:[[BR]][[BR]] |
| 104 | | [[latex($\mu \sim \lambda_{ifp}v_{thermal}$)]][[BR]][[BR]] |
| 105 | | In our code, we offer the option to calculate the realistic viscosity using the following equation:[[BR]][[BR]] |
| 106 | | [[latex($\mu = \frac{m_{i}}{3 \sqrt{2} \pi d_{i}^2} (\frac{8kT}{\pi m_{i}})^{\frac{1}{2}}$)]][[BR]][[BR]] |
| 107 | | The ion mass and local temperature should be readily available, but a good estimate of ion diameter [[latex($d_{i}$)]] usually requires look up a ionic radius table. [[BR]] |
| 108 | | As an estimate, one can use that of hydrogen, about 0.1~0.2 nm.[[BR]] |
| 109 | | One can see that the viscosity itself is anisotropic and depends on the local temperature. [[BR]] |
| 110 | | That is why we have both the [[latex($\mu$)]] and [[latex($\mu'$)]] terms in our viscosity equation, and have to calculate face centered source terms before applying the final divergence.[[BR]][[BR]] |
| 111 | | The viscosity interface is as follows:[[BR]] |
| 112 | | The viscosity source term can be turned on in astrobear by putting the viscosity flag to true:[[BR]] |
| 113 | | '''lViscous = .true.'''[[BR]] |
| 114 | | The second parameter '''VisType''' determines which type of viscosity one requires: 1 for constant viscosity; 2 for real physical viscosity.[[BR]] |
| 115 | | The third parameter '''viscosity''' is the user defined constant viscosity.[[BR]] |
| 116 | | The viscosity module is tested to be working compatibly with implicit thermal conduction module and uniform gravity module.[[BR]][[BR]] |
| 117 | | |
| 118 | | = ======================= = |
| 119 | | = Code Tests = |
| 120 | | = ======================= = |
| 121 | | The 2-D Harris current sheet test can be set up as follows:[[BR]][[BR]] |
| 122 | | [[latex($B_y(x) = b_0 tanh(x/a)$)]][[BR]][[BR]] |
| 123 | | where in computational units we choose [[latex($b_0 = 1$)]] and [[latex($a = 0.5$)]].[[BR]] |
| 124 | | The density profile is chosen so that the pressure equilibrium can be maintained with constant temperature.[[BR]][[BR]] |
| 125 | | [[latex($\rho(x) = \rho_0/cosh^2(x/a)+\rho_c$)]][[BR]][[BR]] |
| 126 | | where [[latex($\rho_0 = 1$)]] and [[latex($a = 0.5$)]]. The temperature is set to be constant as 0.5. [[BR]] |
| 127 | | 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.[[BR]] |
| 128 | | The initial state is in pressure equilibrium though unstable. We put in a perturbation: [[BR]][[BR]] |
| 129 | | [[latex($B_x = B_p sin(kx)cos(ky)$)]][[BR]] |
| 130 | | [[latex($B_y = -B_p cos(kx)sin(ky)$)]][[BR]][[BR]] |
| 131 | | where [[latex($\B_p$)]] is the perturbation amplitude and [[latex($k = 2\pi/L_y$)]]. |
| 132 | | The magnetic field and the density profiles are plotted below:[[BR]][[BR]] |
| 133 | | [[Image(http://www.pas.rochester.edu/~shuleli/hhc_plot.png, 50%)]][[BR]][[BR]] |
| 134 | | The magnetic energy is diminished at the X point which will convert to bulk ram pressure that drives the y direction outflow around the x = 0 line. [[BR]] |
| 135 | | For the outflow pattern, the leading edge is bounded by intermediate shocks and the trailing edge is bounded by slow shocks(Petschek-type). There are totally four slow mode shocks stemming from the Sweet-Parker region, as shown in the following figure.[[BR]][[BR]] |
| 136 | | [[Image(http://www.pas.rochester.edu/~shuleli/petschek_diagram.png, 50%)]][[BR]][[BR]] |
| 137 | | The time evolution model can be given by the Semenov solution, which is obtained by assuming a large resistivity box at the stagnation point. See the following figure.[[BR]][[BR]] |
| 138 | | [[Image(http://www.pas.rochester.edu/~shuleli/semenov_evolution.png, 50%)]][[BR]][[BR]] |
| 139 | | The flow pattern for our setup is shown in the following plot:[[BR]][[BR]] |
| 140 | | [[Image(http://www.pas.rochester.edu/~shuleli/hhc_flowpattern.png, 50%)]][[BR]][[BR]] |
| 141 | | The 2-D mach number plot is shown below:[[BR]][[BR]] |
| 142 | | [[Image(http://www.pas.rochester.edu/~shuleli/hhc_mach.png, 50%)]][[BR]][[BR]] |
| 143 | | |
| 144 | | The following plots show the AstroBEAR test of various settings in the current sheet outflow problem. [[BR]][[BR]] |
| 145 | | [[Image(http://www.pas.rochester.edu/~shuleli/cs_outflow.png, 50%)]][[BR]][[BR]] |
| 146 | | The Workman setting has central symmetry while the reconnection sites locate at the two edges. The positive radius of curvature of field lines at the reconnection site results in an inflow to the X point.[[BR]] |
| 147 | | The first panel shows its kinetic energy distribution. [[BR]] |
| 148 | | The Semenov setting is achieved by applying a resistivity function that varies with position. This resistivity gradient results in bending of magnetic field lines with negative radius of curvature. A vertical outflow from the X point is thus created. The resulted Sweet-Parker region is expanding throughout the simulation till the flow speed at its boundary reaches the Alfven speed. Panels 2 and 3 shows the kinetic energy and vertical flow momentum. Panel 4 shows the kinetic energy of Semenov outflow with an enhanced resistivity contrary. [[BR]] |
| 149 | | |
| 150 | | |
| 151 | | Other possible tests:[[BR]] |
| 152 | | [[BR]][[BR]] |
| 153 | | (1) Flux Explusion Problem (Analytic 2D, resistivity only) [[BR]] |
| 154 | | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_test1.png, 50%)]][[BR]][[BR]] |
| 155 | | This problem studies the field redistribution when there is a rotating rigid cylinder with finite conductivity is placed in the flow. Can be found in text books. |
| 156 | | |
| 157 | | (2) Hartman Flow (Analytic 2D, involves viscosity and resistivity) [[BR]] |
| 158 | | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_test2.png, 50%)]][[BR]][[BR]] |
| 159 | | This problem solves the state of a stationary flow of viscose and resistive fluid between two plates with superimposed transverse magnetic field. Can be found in many MHD text books and papers. |
| 160 | | See for instance: |
| 161 | | [[https://www.pas.rochester.edu/~shuleli/hartmann_apj08.pdf]][[BR]][[BR]] |
| 162 | | |
| 163 | | |
| 164 | | (3) Hydromagnetic Rayleigh-Bernard Convection (2D or 3D. Involves thermal diffusion, viscosity, resistivity) [[BR]] |
| 165 | | [[Image(http://www.pas.rochester.edu/~shuleli/rmhd_test3.png, 50%)]][[BR]][[BR]] |
| 166 | | The stabilizing effect of vertical magnetic field on a classical Rayleigh-Bernard convection problem. This problem involves find the critical Chandrasekhar number with a given Rayleigh number. For a thermal convecting flow at fixed Rayleigh number, the applied vertical non-dimensional magnetic field exceeds the critical Chandrasekhar number. For details, see for instance, Chandrasekhar's 1961 book. |
| 167 | | Also experiments. see for instance: |
| 168 | | [[https://www.pas.rochester.edu/~shuleli/RB_pre00.pdf]][[BR]][[BR]] |
| 169 | | |
| 170 | | The following LLNL journal contains useful tests on resistive MHD solver: [[BR]][[BR]] |
| 171 | | [[https://e-reports-ext.llnl.gov/pdf/374198.pdf]][[BR]][[BR]] |
| 172 | | |
| 173 | | Tests |
| 174 | | The stabilization effect of RT by constant viscosity and constant thermal conduction.[[BR]][[BR]] |
| 175 | | [[Image(http://www.pas.rochester.edu/~shuleli/rt_graph.png, 50%)]][[BR]][[BR]] |
| 176 | | |
| 177 | | The Harris Sheet Test Simulation:[[BR]][[BR]] |
| 178 | | Click the following links for animations of kinetic energy and magnetic mach.[[BR]] |
| 179 | | [http://www.pas.rochester.edu/~shuleli/currentsheet/csheet_kin.gif GIF:Kinetic Energy][http://www.pas.rochester.edu/~shuleli/currentsheet/csheet_mach.gif GIF:Magnetic Mach][[BR]] |
| 180 | | [[Image(http://www.pas.rochester.edu/~shuleli/currentsheet/csheet_test.png, 40%)]][[BR]][[BR]] |
| 181 | | |
| 182 | | |
| | 125 | Notice that in our island formation problem, the perturbation is anti-symmetric about the center. So the density pattern is also anti-symmetric. It is trivial to use a symmetric perturbation in order to obtain a symmetric island pattern. [[BR]] |
