Solution of the Blatter-Pattyn model
The final form of the equations we'd like to solve is:
Again, note that for the x equation we've moved all the terms containing gradients in v to the right-hand side (RHS).
We've set it up this way in order to solve the equations using an operator splitting approach; for the x equation, we treat v as known (where we take the values of v from the previous iteration) and solve for u, and vice versa when we solve they y equation for v. The "splitting" refers to the fact that we are breaking the multi-dimensional divergence operation into multiple steps. Rather than solving one big matrix equation for u and v simultaneously we solve two smaller matrix equations in sequence with one of the unknowns treated as a known "source" term.
As with the 0-order model, we need to change from Cartesian to sigma coordinates. The first normal-stress term first term on the left-hand side becomes
where hatted values refer to the coordinate directions in sigma coordinates. Similarly, the first cross-stress term on the RHS is given by
This looks pretty scary eh? Yes, there is a lot of work to be done here. But, luckily, there is also a lot of symmetry here. Notice that if we wanted to design subroutines to discretize these equations, we could re-use a lot of them in multiple places by passing the appropriate arguments. For example, in many places, the only difference between the two equations is whether or not we are differentiating u or v (so apply the disrectization to either u or v) and/or whether or not differentiation is with respect to x or y (so pass the argument, either x, or y, where appropriate).