Particle-in-cell finite element analysis
at the University of Sydney - Division of Geology & Geophysics

(this page is always under construction -
its author is in the process of learning how to make webpages,
so please excuse its appearance!)

Our geophysics group (headed by Dr. Dietmar Muller) uses the particle-in-cell finite element code "Ellipsis" created by Dr. Louis Moresi (then at CSIRO in Nedlands, WA, now with Monash University's School of Mathematical Sciences).  As part of a nationwide effort to improve/utilize the code, we have extended Ellipsis to 3-dimensions, and are in the benchmarking stages to make sure that it is reliable for geophysical applications.

In a nutshell, the particle-in-cell method combines the strengths of the Lagrangian and Eulerian formulations of mechanics, while bypassing their limitations. For a bit more detail, check out the (~1 page, pdf) article by Dr. Rich Albert in the SGSEG Newsletter (page 4) from June 2002, or consult Dr. Moresi's webpage for some very informative descriptions of the methodology. The computational cost of the PIC procedure is higher than a standard finite element analysis, but oftentimes the additional cost is a small price to pay relative to the advantages that the method provides.

An extensive amount of work has been done with the 2D version of the code, by postgraduate student Craig O'Neill. At some point we will post some animations showing some of Craig's work.

The 3D code is the focal point of by Rich's research, and is to be applied to several problems, including extension of the lithosphere (as per the rifting process and the formation of passive margins), as well as mantle convection and its interaction with Earth's crust. Below are some animations for relatively simple 3D simulations. If you are unable to watch the movies but would like to see them, you can download a viewer at: http://www.apple.com/quicktime
 
 
final image of convection simulation (so far)
MPEG movie #1 (~0.49 MB)
3D animation of the evolving temperature of the two strips
MPEG movie #2 (~0.56 MB)

Above (left) is shown the final stage in a simulation of convection in 3D, where the box is sliced in 3 locations, as shown in this pdf schematic.  This simulation has a temperature-dependent viscosity, and is driven by a thermal perturbation. The green and blue strips are markers of material, with identical properties as the fluid in the box. Of course, the properties could be made to be different, but the main point is to note how the PIC method allows tracking of particles throughout a loading programme. At a later date we will run the problem for longer so that the deformation is more extensive. To see a movie with the 3 slices, click Movie#1.

With a bit more post-processing, we can restrict the animations to, say, just the 2 strips and see what the how temperature of each particle evolves during the simulation. A 3D animation of the two strips in motion (as shown in the above right image) is seen by clicking Movie#2. Animations can be done of iso-variable contours, too, not just for a selected material; it's all just a matter of a few changes in an awk script.
 
 
Extension of a layered lithosphere via 3D PIC finite element analysis
MPEG movie #3 (~1 MB)
A vertically oriented slice through the extended 2-layer system
MPEG movie #4 (~4.0 MB)

Above (left) is pictured the initial and final configurations of a simulation with a 2-layered lithospheric block, although only the upper layer is shown. The colors show the amount of inelastic strain for the top layer. In this simulation, we put two zones of weakness right through the upper layer (seen as the two 'gaps' at the top surface in the initial configuration - although the material in these gaps is present, it has just been removed to facilitate seeing where the zones are located), and subjected the block to 2 separate phases of extension. The first phase occurs in a direction normal to the leftmost face, starting from rest and stopping at a later time. Then after a pause, the second phase of extension begins, this time the direction of extension has components in both horizontal directions. You can see how the block has appeared to drop in elevation, which is just the principal of conservation of mass in action, as the layers thin during stretching. It is apparent that the presence of the staggered weak zones serve to focus futher inelastic deformation, as one would expect. Click Movie#3 for an mpeg animation of this simulation. Another animation (see the image to the above right) of a vertically oriented 2D slice through the 2 layers (upper layer plotted with circles, lower layer with triangles) can be viewed by clicking Movie#4. One of the weak zones can be seen in the image and movie, shown as the red triangles which extend through the upper layer.

Our thanks to the creators of the GMT package, which was used in conjuction with awk and sed to create the images for the animations.

contact us via e-mail for questions/comments: albert@es.usyd.edu.au