110 likes | 299 Views
Numerical modeling of rock deformation 11 :: FEM for heterogeneous media. www.structuralgeology.ethz.ch/education/teaching_material/numerical_modeling Fallsemester 2011 Thursdays 10:15 – 12:00 NO D11 & NO CO1 Marcel Frehner marcel.frehner@erdw.ethz.ch , NO E3 Assistant: Jonas Ruh, NO E69.
E N D
Numerical modeling of rock deformation11 :: FEM for heterogeneous media www.structuralgeology.ethz.ch/education/teaching_material/numerical_modeling Fallsemester 2011 Thursdays 10:15 – 12:00 NO D11 & NO CO1 Marcel Frehner marcel.frehner@erdw.ethz.ch, NO E3 Assistant: Jonas Ruh, NO E69
Goals of today • Go from homogeneous mediato heterogeneous media • Do a lot of exercises!
What we did so far... • 1D steady state • 2D steady state • 2D elastic (homogeneous media) • 2D viscous (homogeneous media) What we do today • 2D viscous (heterogeneous media) Almost nothing changes…
Homogeneous 2D elastic example • Elastic beam bending under gravity
Homogeneous 2D viscous example • Viscous block flowing on the floor
Heterogeneous 2D viscous example • Diapirism
Heterogeneous 2D viscous example • Folding
How to treat heterogeneous media with FEM...? The finite-element formulation was always done for one single element: • Multiplication with test functions • Spatial integration • Integration by parts(& dropping boundary terms) • FEM-approximation withshape functions(different for v and p) • Taking nodal values out of integration • Reorganization of equation Sum up local equations and assembleglobal system matrices Deriving weak form Local element
The finite-element equations for viscous flow Unknowns Local element Force vector
Heterogeneous medium • Now, we simply use different material parameters in the different elements! h2D2 r2 h1D1 r1
The Phase array % layers from bottom to top mu = [ 1 100 1 ]; % viscosity [Pa s] rho = [ 1 1 1 ]; % density [kg/m3] Ly_array = [ 15 1 15 ]; % thickness of layers nely_array = [ 10 10 10 ]; % #elements in y-dir ... for iely = 1:nely for ielx = 1:nelx iel = iely*nelx-nelx+ielx; EL_N(1,iel) = ... EL_PTS(1,iel) = ... Phase(iel) = sum(sign(iely-cumsum(nely_array))==1)+1; end end ... Kloc = Kloc + B‘ *D(:,:,Phase(iel)) *B *detjac*wts(k); Gloc = Gloc - B_G‘ *shape_p(k,:) *detjac*wts(k); Mploc = Mploc + shape_p(k,:)‘ *shape_p(k,:) *detjac*wts(k); Floc(2:2:dof_per_el) = Floc(2:2:dof_per_el) + shape(k,:)‘ ... *rho(Phase(iel))*gravity *detjac*wts(k); Phase = [1 1 1 2 2 2]