Hi, I'm working on the dynamic of alpine glaciers and I'm trying to resolve the full equations of Stokes with finite differences. Is there anybody that knows an algorithm to solve the full equations of Stokes without approximations? For the moment I'm working in 2 dimensions (x,z) and I have so the equations:

-dp/dx + dtau_xz/dz = - dtau_xx/dy dp/dz + rho*g = dtau_xz/dy + dtau_zz/dz

where p is the pression, rho the density of the ice, g the acceleration of earth and tau_ij are the deviatoric stresses. I have the boundary condition v=0 at the bottom and tau*n=0 at the surface (tau is the matrice of the tau_ij's an n is a vector normal to the surface).

And finaly I have: eps_ij=1/2*(dv_i/dj+dv_j/di)=A*tau_star^m*tau_ij, where A is a const, m a paramenter (I tried 0 -> Newton, 2-> Glen) and tau_star the invariant of the stress tensor.

Can anybody help me? I tried an iterative algorithm to get the velocity fields proposed by K.Hutter, but for the moment it doesn't converge!

Is this possible with finite differences or have I to use other methods?

Thanks a lot.

For stokes flow you can try using boundary elements method.Also you need to be careful when applying finite difference schemes. I would suggest using an implicit scheme for non viscous space terms central difference scheme for viscous terms and explicit fourth order runga kutta for time marching .

-H

Thanks for your response. I do not yet trying to iterate in time, the aim for the moment is to determine the velocities for a given (fixed) geometry. What do you mean by "boundary elements method"? Ok for the different schemes, but with finite differences I can't solve my equations "straight forward" without makeing approximations, can you give me a hint? Thanks. Martina

For de transport-convectives terms : QUICK (warning in the numerical diffusion)