Changes between Version 25 and Version 26 of AstroBearProjects/resistiveMHD


Ignore:
Timestamp:
01/28/13 00:26:02 (12 years ago)
Author:
Shule Li
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • AstroBearProjects/resistiveMHD

    v25 v26  
    22
    33The resistive equation can be written as:[[BR]]
    4 
    54[[latex(\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B}))]]
    6 
    75The 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 
     6Now 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]]
     7field configuration in equilibrium is subject to strong diffusion, usually heating would occur and surppress the local resistivity and thus the  [[BR]]
     8diffusivity. By expanding the first equation, we have the form:[[BR]]
    119[[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 
     10Now 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]]
     11on each other.[[BR]]
     12So 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  [[BR]]
     13explicit solver to treat the problem instead.[[BR]]
     14One 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  [[BR]]
     15with 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]]
     16with temperature dependence built in. [[BR]]
    1917[[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 
     18Unfortunately, this does not work because the diffusion equation itself has to be divergence free. Our roughest approximation, treating the resistivity  [[BR]]
     19as 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]]
     20we end up getting:[[BR]]
    2321[[latex(\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B})]]
    24 
    2522The 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]]
    2623So 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 
     24Upon 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  [[BR]]
     25total energy during the evolution, especially for the force-free field case. This phenomenon is explained below. [[BR]]
    3026In AstoBEAR, we explicitly calculate the resistivity induced current on the cell edges, following equation: [[BR]]
    3127[[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 
     28The 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  [[BR]]
     29calculating the diffusive current on the x direction, notice that the red arrow is where we are calculating the diffusive current, the green arrows are  [[BR]]
     30where the magnetic field originally resides: [[BR]]
    3431[[Image(resistive_diagram.png, 30%)]] [[BR]]
    35 
    36 
    37 The actual code looks like: [[BR]]
    38 
    39 << code. [[BR]]
    40 
     32The actual code looks like the following in 2D (jy are initialized to 1): [[BR]]
     33{{{
     34    DO i=mjy(1,1),mjy(1,2); DO j=mjy(2,1),mjy(2,2); DO k=mjy(3,1),mjy(3,2)
     35       IF(nDim==2)THEN
     36          jy(i,j,k) = -jy(i,j,k)*(Info%q(i,j,k,iBz)-Info%q(i-1,j,k,iBz))/dx
     37       ELSE
     38          jy(i,j,k) = -jy(i,j,k)*((Info%aux(i,j,k,3)-Info%aux(i-1,j,k,3))+(Info%aux(i,j,k-1,1)-Info%aux(i,j,k,1)))/dx
     39       END IF
     40    END DO; END DO; END DO
     41}}}
    4142After calculating and storing jx, jy and jz of the given grid, the change of magnetic field due to diffusion can be calculated as： [[BR]]
    42 
    4343[[latex(\frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J})]] [[BR]]
    44 
    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]]
     44In 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]]
     45are only in z. So we need to use f instead. The code looks like: [[BR]]
     46{{{
     47    IF(nDim==2)THEN
     48       ALLOCATE(fjx(mjy(1,1):mjy(1,2),mjy(2,1):mjy(2,2),mjy(3,1):mjy(3,2),1))
     49       ALLOCATE(fjy(mjx(1,1):mjx(1,2),mjx(2,1):mjx(2,2),mjx(3,1):mjx(3,2),1))
     50       ALLOCATE(bflux(1))
     51       bflux(:)=iBz
     52       fjx(:,:,:,1)=jy(:,:,:)*dt/dx
     53       fjy(:,:,:,1)=-jx(:,:,:)*dt/dx
     54       CALL storefixupfluxes(Info,mjy,1,fjx,bflux)
     55       CALL storefixupfluxes(Info,mjx,2,fjy,bflux)
     56       DEALLOCATE(bflux)
     57       DEALLOCATE(fjx); DEALLOCATE(fjy)
     58    ELSE
     59       CALL StoreEmfs(Info,mjx,1,jx*dt/dx)
     60       CALL StoreEmfs(Info,mjy,2,jy*dt/dx)
     61    END IF
     62}}}
     63In 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  [[BR]]
     64the fluid is everywhere zero, and no heat is generated by the current. [[BR]]
    5065If we dot the resistive induction equation with the magnetic field B, we obtain the time evolution equation for magnetic energy:[[BR]][[BR]]
    5166[[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]]
     67where [[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]]
     68of the diffusive current.[[BR]]
     69In 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]]
     70the loss of magnetic energy due to reconnection.[[BR]]
    5471The total energy change for the resistive step is therefore:[[BR]][[BR]]
    5572[[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]]
     73Here 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]]
     74is converted into thermal energy. [[BR]]
    5775In the code, the energy flux as a result of magnetic diffusion needs to be calculated explicitly using: [[BR]]
    58 
    5976[[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]]
     77The 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]]
     78current 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]]
     79dashed lines are what used to compute the energy flux. [[BR]]
    6280The magnetic field can be updated from the diffusive currents by: [[BR]]
    63 
    6481[[latex(\frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J})]] [[BR]]
    65 
    6682The 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 
     83{{{
     84    DO i=msx(1,1),msx(1,2); DO j=msx(2,1),msx(2,2); DO k=msx(3,1),msx(3,2)
     85       IF(nDim==2)THEN
     86          Info%aux(i,j,k,1) = Info%aux(i,j,k,1)+(jz(i,j,k)-jz(i,j+1,k))*dt/dx
     87          xflux(i,j,k,1)=-0.25*(jz(i,j+1,k)+jz(i,j,k))*(Info%q(i-1,j,k,iBy)+Info%q(i,j,k,iBy))
     88       ELSE
     89          Info%aux(i,j,k,1) = Info%aux(i,j,k,1)+((jy(i,j,k+1)-jy(i,j,k))+(jz(i,j,k)-jz(i,j+1,k)))*dt/dx
     90          xflux(i,j,k,1)=0.25*(jy(i,j,k+1)+jy(i,j,k))*(Info%q(i-1,j,k,iBz)+Info%q(i,j,k,iBz))-0.25*(jz(i,j+1,k)+jz(i,j,k))*(Info%q(i-1,j,k,iBy)+Info%q(i,j,k,iBy))
     91       END IF
     92    END DO; END DO; END DO
     93}}}
    7094We then call storefixupfluxes for the three energy flux components (as xflux, yflux and zflux in the code).[[BR]]
    7195Finally, the energy is updated using the divergence of the energy fluxes. This finishes the diffusive process.[[BR]]
    72 [[BR]]
    73 
    7496The resistive time scale is explicitly calculated by: [[BR]]
    75 
    7697[[latex(dt_{resistive} = \frac{c dx}{\eta})]] [[BR]]
    77 
    7898where 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]]
     99{{{
     100    explicit_maxspeed(level) = max(explicit_maxspeed(level),2d0*resistivity/dx)
     101}}}
     102where 2d0 can be replaced by any value above 1 (for example, putting this number to 5d0 will reduce the cfl to 0.2).[[BR]]
    83103
    84104'''Resistivity Test Problems'''[[BR]][[BR]]
    … …  
    96116The 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]]
    97117[[Image(http://www.pas.rochester.edu/~shuleli/hhc_plot.png, 30%)]][[BR]][[BR]]
    98 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 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 lines continue to come in and annihilate because of the lower field pressure at the center. The reconnection heat will drive outflows out of the X point, parallel 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 diagrams show how increased resistivity at the center of a sheer pinch drives Petschek shock.[[BR]]
     118The 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]]
     119at 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]]
     120lines continue to come in and annihilate because of the lower field pressure at the center. The reconnection heat will drive outflows out of the X point, parallel [[BR]]
     121to 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]]
     122diagrams show how increased resistivity at the center of a sheer pinch drives Petschek shock.[[BR]]
    99123[[Image(http://www.pas.rochester.edu/~shuleli/multi_test/resistive_instab.png,30%)]][[BR]][[BR]]
    100124The flow pattern for our setup is shown in the following plot:[[BR]][[BR]]
    … …  
    103127[[Image(http://www.pas.rochester.edu/~shuleli/hhc_mach.png, 30%)]][[BR]][[BR]]
    104128
    105 To test Sweet-Parker problem in AstroBEAR, we . [[BR]][[BR]]
     129To test Sweet-Parker problem in AstroBEAR, we construct the sheer pinch using the above setup, and modify . [[BR]][[BR]]
    106130
    107 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 by white lines. [[BR]]
     131The 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]]
     132by white lines. [[BR]]
    108133[[Image(http://www.pas.rochester.edu/~shuleli/multi_test/sp1_0020.png, 30%)]]
    109134[[BR]]
    … …  
    116141where Bp is the perturbation amplitude and the wave number:[[BR]]
    117142[[latex(k = 2\pi/L_y)]].[[BR]]
    118 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 of pressure imbalance, which creates periodical dense "islands".
    119 The growth rate and the size of the "islands" depend on resistivity and the strength of the perturbation.[[BR]][[BR]]
     143This 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]]
     144of 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]]
    120145
    121146[[Image(http://www.pas.rochester.edu/~shuleli/multi_test/mihr_0180.png,40%)]][[BR]][[BR]]
    … …  
    123148To watch a full movie, [http://www.pas.rochester.edu/~shuleli/multi_test/mihr.gif click here][[BR]][[BR]]
    124149
    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]]
     150Notice 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 [[BR]]
     151a symmetric perturbation in order to obtain a symmetric island pattern. [[BR]]