My code in PyLith is now successfully running some simulations that are responding to changing boundary conditions. a bit of analysis that Brad and I did showed that the ~5cm displacement was only about 15-20% off of the estimated displacement for the end of the ice shelf. Additionally, the grounding zone is showing a beautiful curve towards the flotation height.
After several weeks of slow and frustrating progress, we are finally nearing real data!
This week was a bit of a bust. I didn’t really end up making much progress at all. I am still not able to get PyLith to work properly. I am having a lot of trouble getting gravity to work properly. At Brad’s advice, I rebuilt the simulation starting from a similar example provided in the code, changing just one element at a time. Unfortunately, I am still having big problems with gravity.
PyLith does not currently allow the user to specify any type of body force, only pressure, displacement, or point force conditions on “boundaries” which can be any collection of points in the mesh. But, in order to “turn on” gravity, the user must use the provided gravity field that is capable of applying body forces.
Unfortunately, the gravity is causing strange happenings in my simulations. It does not appear to be spatially uniform and is causing impossible curves. But, the σ-yy component of stress (the compressive stress in the y diretion – so basically how much any little section of ice is pressing against the ice above and below it) seems correct in distribution and magnitude.
The other problem that I have run in to is that the shape of my models that don’t include gravity but only a pressure condition don’t have the correct shape either.
Now I am going to go into the actual equations that we are dealing with so if you don’t want to deal with the math you can skip down a bit, but I promise it’s pretty simple.
So, because the code already builds in the way a material responds to forces, all the user needs to do is input the forces and material properties. One of the first things I did was figure out what kind of boundary conditions I need. This included the more simple conditions of fixing one end of the beam in place to simulate the part of the ice that is stuck on the ground, and having no normal or shear forces on the top free surface. This means that the air is not exerting a pressure directly into the ice or along the surface (friction is a good example of a shear force). I am neglecting air pressure because it is so small and also ends up cancelling out for the most part.
So the trickiest condition is the bottom boundary pressure. Since I only care about what is happening at the very bottom of the ice, the pressure condition does not depend on the height of the ice above the boundary since this is constant, the depth of the water is the only thing that matters.
Acting down on the ice is the weight of the overlying water column is
p = ρgH
where ρ is the density of ice (917 kg/m^3), g is the gravitational constant of acceleration (9.81 m/s^2), and H is the thickness of the ice (for now constant at 700 m).
And the force acting up is the hydrostatic pressure which is equal to the weight of the overlying water column which can be understood with a little bit of fluid dynamics stuff. The way it works is we say that the liquid is always in hydrostatic equilibrium which means that there is no net acceleration happening in the fluid, which is a good assumption for water because it has low viscosity and therefore can flow much faster than the forces are being changed. The important condition that comes from hydrostatic equilibrium is that the pressure at any given depth is the same so the pressure is the same at the very bottom of the ice, where it is being acted upon by only ice and a little to the side at the same depth where it only is sensitive to the water above it (see points A and B in diagram below).
The height of the water is broken up into three components: (1) the average neutral height, (2) the tidal height change, and (3) the deflection from the neutral position. The average height is the position at which ice is in hydrostatic equilibrium at the average tidal height , determined with the assumption that the grounded ice is at the average height such that at average tidal height there is no deflection. this is calculated by setting the other two components to zero and solving the equation 0 = p_ice + p_water = -ρgH + ρ’gh for h.
This results in h = (ρ/ρ’)H meaning that h is the height of the ice times the ratio of the density of ice to the density of sea water. Having this constant makes it easy to represent the force on the bottom of the ice as a function of the tidal height and the deflection of the ice from equilibrium. This last part is a restoring force that increases the pressure acting upwards on the ice if it is below the equilibrium height and decreases that pressure if the ice is floating too high.
The pressure can then be expressed as
p = -ρgH + ρ’g(h + Δh – (H+y))
and substituting for h:
p = -ρgH + ρ’g((ρ/ρ’)H + Δh – (H+y)) = ρ’g( Δh – (H+y)) which is a buoyancy force that interestingly only depends upon the density of the sea water.
So, that was a lot of math, but what ultimately is important to take away from this is that the pressure is dependent upon the instantaneous deflection of the beam and varies along the surface, which I think is where I am having issues with PyLith.
While, at the moment, I still have some problems to work out but, as described above I think I have found their sources, the real issue is trying to fix them. I am getting a bit frustrated with the modelling, but I now understand it well enough that hopefully once the problems are resolved i will be able to progress quite quickly in increasing the complexity and hopefully accuracy of the model.
As usual, I am still working on getting the PyLith to properly run the simulations. I am really looking forward to the point at which the simulation runs and gives me values of the correct magnitude (right now the displacements are about 700 times larger than they should be) so that I can start actually looking at what happens when I vary the parameters.
PyLith is currently computing reasonable stresses for the first step. From there, I had to download a bunch of stuff for python so that I could run the script that lets me take the data from the first step to feed it to the time dependent part of the problem. This database is still giving me problems.
Hopefully next week I can get the 2D problem work.
I had a very productive week 3. I spent most of the week debugging and was finally able to get the program to run correctly after a few tweaks suggested by my mentor, Brad.
PyLith now computes the first time step to determine the initial stresses and displacements of a flexible beam with the properties of ice fixed in position at one end. I am very excited to work on the next step of the project which is to make pressure depth dependent and find the equilibrium position for the ice beam which should be close to the analytical solution found in Matlab.
I have made a lot progress in the last couple of days since my previous blog post.
I now have and can operate all of the software I need for my program including Trelis, which makes meshes to use for the finite element code PyLith, and Paraview which helps to visualize data. I have not yet written anything that works, but hopefully I can begin by replicating the model I made in Matlab using an analytical solution.
I now have Matlab running some fun animations of the elevation profiles of a simple elastic beam model for the ice. It draws the profile out from the grounding line (the zone where the ice last touches the ground before floating on the surface of the water) to the edge of the ice. The model calculates profiles at various different heights of the water, although it needs an inertial component in order to model the motion in time.
My next step is going to try to get the same solution out of the finite element code before trying to bring in some neglected terms and geometric complexities. Additionally, this may require some rethinking about the boundary conditions of the beam.
I have completed the first week of my SESUR project. I began working on Wednesday by downloading PyLith a finite element code which finds numerical solutions to equations with an emphasis on earthquake modeling.
I got familiar with the software by reading select portions of the user manual and running and reading through example problems provided with the download.
Similarly, I began working with ParaView which is a visualization software to help understand scientific data. It can run 2D and 3D renderings as well as animations for time dependent processes.
I still need to get software to generate mesh inputs for PyLith.
My graduate student mentor Brad Lipovsky and I discussed a model of an ice shelf that follows a simple beam flexure problem with numerical solutions. I began setting up the beam flexure problem for PyLith beginning with detailing the material properties. We also discussed the governing equations and boundary conditions which I will set up next week.