Solar Numerical Simulations - A 3 Week Research Project (And Blog Hiatus)

In the previous blog post, I laid out the basis for how numerical simulations can be done by the process of "discretizing" partial differential equations (which afford no analytical solutions).  This technique also goes by the name of "finite differences" and incorporates an array of methods - such as the FTCS, or 'forward time-centered scheme - which I applied to an aurora using basic assumptions, of which the most important was the change from a particulate to a fluid dynamic treatment.

I also noted that the case of solar applications is generally more difficult. Readers can get an excellent view (graphically) of the nature of the basic shear angle and magnetic shear stress problem (in the vicinity of the solar photosphere) by reference to this earlier blog link:

http://brane-space.blogspot.com/2010/09/predicting-next-giant-solar-flares.html

Note, the upper image of a massive solar active region and the x-y axes, and the direction of the magnetic neutral line (NL) relative to the x -axis especially. This neutral line (also called the "magnetic inversion line" ) separates regions of (+) and (-) polarity.  In general, the degree of magnetic stress, indicating how much free magnetic energy is stored in the region, is given by what we call the "shear angle", denoted φ.   A rough estimate of the shear angle magnitude in the image can be obtained by drawing a tangent line from the origin of the superposed x-y coordinate system to the upper 'hump' of the neutral line, just below the letters, 'NL'. In general, the approach of a flare is signified by angles   φ > 84 degrees.

In solar numerical approaches to magnetic shear, the central partial DE which is attacked is given in terms of the vector potential A:

DIV^2 A  + d/dA {½ B(y)^2} = 0

where B(y) is the magnetic induction or magnetic component in the y direction. Historically, two problems have been defined: 1) B(y) = f(A) for whcih the functional form of the field is pre-prescribed (obviously simpler) and 2) d = d (x) or the displacement d(x) of the footpoints (say in an arch model) from the x-axis are imposed. One then follows the evolution of the field through a series of force-free equilibria.

The first numerical efforts at solving problem (2) were made by Sturrock and Woodbury and Barnes and Sturrock (1972), but they only found one solution.

My goal, in a research project over the next 2 1/2 weeks, is to find additional solutions, numerically. I also seek to have these solutions consistent also with those solutions of  the shear velocity field (briefly described in the link), thus:

v (φ,W) =  a sin (φ ) exp [1 - 2[φ ]/ W

where a defines a normalization constant peculiar to the region, and generated from a mean of shearing scales by comparing regions in the same cycle. Then, [ φ] is the absolute value of the shearing angle (defined from the netural line orientation to the x-y axis) and 'W' is the assigned cell size of the shearing region in appropriate units.

Selection and discretization of the appropriate equations will be no easy task.  For example, the discreteness of the numerical algorithm always introduces modifications to the equations which are to be solved.   In addition, explicit finite difference methods for hyperbolic equations (see end of previous blog)are subject to a stability limit for the time step:


(delta t)  < O(delta x/ v (ty))

known as the Courant condition (or the CFL, i.e. Courant-Friederichs-Levi condition).  Explicit differencing implies that the update is obtained with quantities known in the vicinity of the particular node i, i.e., usually requires nodes i-1, i, and i+1.(Example - in the 2nd diagram shown in the link a1 might be node i-1, and a2 at i+1 in an exaggerated sense.)  Assume information is carried in the simulation with a velocity v(ty) implying that the information is carried a distance L = v(tp) *(delta t ) in a single time step delta t. However, if L is much larger than the grid separation delta x then it travels more than a single grid spacing in one time step.

What one gets is a delicate balancing act and almost continuous cross-checking with each test check of a discretizaton equation in a computing run.  For example, the maximum velocity that is resolved by an explicit scheme is given by v_max = delta x/ delta t yielding the condition: delta  t < v(max)/delta x as in the Courant condition above. In effect, if this tested v_max conflicts with the magnitude of  my shear velocity for a given region, v (φ,W), there is the warning of a numerical instablity. If care isn't taken then, the errors grow exponentially and one can end up with a meaningless result.

The long and short of it is that this will be my focus the next two and a half weeks, and hence no further blogs until about May 22nd. I do hope readers will understand, and perhaps avail themselves of some past blogs they may only have scanned over before.   This is a project I've been eyeing completing for some time, and I believe the time is now right.  To regular readers, and followers, I regret any inconvenience!