Numerical modeling of rock deformation 11 :: FEM in 2D for viscous materials

1. Numerical modeling of rock deformation11 :: FEM in 2D for viscous materials www.structuralgeology.ethz.ch/education/teaching_material/numerical_modeling Fallsemester 2013 Thursdays 10:15 – 12:00 NO D11 (lectures) & NO CO1 (computer lab) Marcel Frehner marcel.frehner@erdw.ethz.ch, NO E3 Assistant: Marine Collignon, NO E1

2. Goals of today • Repeat the finite element derivation • Remember the 4 steps in the FE-derivation:  Multiplication with test functions  Spatial integration  Integration by parts  Finite-element approximation (Galerkin method) ( Reorganization) • Do everything equally for 2D viscous rheology • Learn the v-p-formulation • Deal with incompressibility (4 steps have to be performed twice) Weak form

3. Recap from 2D continuum mechanics and rheology Equations Matrix notation Short form Incompressibility Conservation of linear momentum (i.e., force balance equation)Rheology Kinematic relation

4. Total system of equations ! NEW ! • Incompressibility • Short form • Matrixnotation • Written out(2 equations)

5. Force balance equation: Deriving the weak form • Test functions:  Multiplication with test functions:  Spatial integration

6. Force balance equation: Deriving the weak form  Integration by parts anddropping arising boundary terms This is the weak form of the force balance equation!

7. The finite element approximation 2D elasticity 2D viscosity • 9 bi-quadratic shape functionsfor velocity • 4 bi-linear shape functions for pressure (the same as for 2D elasticity) • Requirements for the shape functions: • Equal 1 at one nodal point • Equal 0 at all other nodal points • Sum of all shape functions has to be 1 in the whole finite element.

8. The Q2P1-element: bi-linear pressure shape functions

9. The Q2P1-element: bi-quadratic velocity shape functions Element for velocity

10. Force balance equations: The FE-approximation  Apply finite element approximation • Take nodalvalues out ofintegration

11. Force balance equation: Some reorganization • Writing everything as one equation in vector-matrix form • Reorganization

12. Force balance equation: The final equation Apply numerical integration using 9 integration points

13. Incompressibility equation: Steps  –  • Original equation: • Test functions (use pressure functions):  Multiplication with test functions:  Spatial integration: • NO integration by parts!!!  Apply FE-approximation:

14. Incompressibility equation: Some reorganization • Writing everything as one equation in vector-matrix form • Reorganization Apply numerical integration using 9 integration points

15. Final set of equations • Force balance: • Incompressibility: • Everythingtogether: • Voilà:

16. What’s new? • New physical parameters (viscosity instead of elastic parameters) • Numerical grid and indexing is different because of bi-quadratic shape functions (9 nodes per element).So, you need to change • EL_N • EL_DOF You need to newly introduce • EL_P, which provides a relationship between element numbers and pressure-node numbers • New set of shape functions (4 bi-linear for pressure stay the same,9 bi-quadratic for velocity are new) • New set of integration points (9 instead of 4) • Calculate G inside the loop over integration points in a similar way as K. For this you need the shape function derivatives in a vector format (instead of matrix B). • Put K, G, and F together to form the super-big K and super-big F. • Your solution is now a velocity. So, introduce a time increment Dt and update your coordinates as GCOORD = GCOORD + Dt*velocity. Then you can write a time-loop around your code and run it for several time steps.

17. Problem to be solved today and next week • Homogeneous medium • Boundary conditions • Left and bottom: Free slip (i.e., no boundary-perpendicular velocity, boundary-parallel velocity left free) • Right and top: Free surface (both velocities left free) • Calculate several time steps