Slowness surface calculation for different media using the symbolic mathematics language MAPLE® EARTH SCIENCES RESEARCH JOURNAL Earth Sci. Res. J. Vol. 8, No. 1 (Dec. 2004) : 63 - 67 Manuscript received July 2004. 63 Paper accepted October 2004. SLOWNESS SURFACE CALCULATION FOR DIFFERENT MEDIA USING THE SYMBOLIC MATHEMATICS LANGUAGE MAPLE® Miguel Duarte1, Carlos Piedrahita1, Trino Salinas1, Hernando Altamar 2 and Karen Pachano 2 1 Instituto Colombiano de Petróleos , ECOPETROL. Piedecuesta – Colombia. E-mail: mudarte@ecopetrol.com.co; tsalinas@ecopetrol.com.co. 2 Universidad Industrial de Santander -UIS. Bucaramanga- Colombia ABSTRACT Starting from the equation in different media, we obtain the different type of waves that can exists in such media. The evaluation of the eigenvalues and eigenvectors let us construct the slowness surfaces. In general complex calculations case have to be made. In this work, routines were implemented in the symbolic language MAPLE® and the slowness surfaces were plotted. This work is relevant for the modelling of equivalent media that simulates natural fractured reservoirs, like those common in the Colombian foothills. It is important to understand the seismic response of this reservoirs for exploration of this areas. Key Words: Christoffel matrices, slowness vector, cubic symmetry, hexagonal symmetry, isotropic, symmetry, slowness surfaces, eigenvalues, eigenvectors RESUMEN Partiendo de la ecuación en diferentes medios, obtenemos los diferentes tipos de ondas que pueden existir en tales medios. La evaluación de los valores propios y los vectores propios permite construir las superficies de lentitud. En general, es necesario hacer complejos casos de cálculos. En este trabajo se implementaron algunas rutinas en el lenguaje simbólico MAPLE® y se señalaron las superficies de lentitud Este trabajo es pertinente para el modelado de medios equivalentes que simulan los depósitos fracturados naturales, como los que aparecen en las zonas de piedemonte colombianos. Es importante entender la reacción sísmica de estos depósitos para la exploración de estas áreas. Palabras clave: isotropía, matriz de christoffel, simetría cúbica, simetría hexagonal, superficies de lentitud, valores propios, vectores propios, vector de lentitud. © 2004 ESRJ -Unibiblos. Miguel Duarte et al. METHODOLOGY This work was developed entirely in the field of the linear theory of elasticity. For homogeneous media the solution of the general elastodynamic equation reduces to the solution of the Christoffel Matrix eigenvalues problem. The Christoffel equation is solved numerically using the symbolic mathematics programming language Maple®. The phase velocity is obtained as a function of the anisotropic parameters and the propagation direction. The phase velocity 1/v is used to obtain the slowness, which is plotted using the standard plotting Maple function in polar coordinates. The resulting graphic shows the relation between the slowness vector and the propagation angle. FORMULATION OF THE PROBLEM AND RESULTS Theory Research – Fundamentals The deformation laws for a linear elastic body are determined by the linear momentum equation (1), the constitutive equation or the Hooke’s Law (2) and the small strain – displacement linear relation (3) (Sen, 2002), (Krishnaswamy, 2004), (Pšenčík, 1994): 2 ij i j u u x t τ τ ρ ρ ∂ ∂ ∇⋅ = ↔ = ∂ ∂ && (1), : ij ijkl klC Cτ ε τ ε= ↔ = (2), ( )1 1 2 2 k l kl l k u u u u x x ε ε ⎛ ⎞∂ ∂ = ∇ + ∇ ↔ = +⎜ ⎟ ∂ ∂⎝ ⎠ (3). Where: τ is second order stress tensor, ρ is density u is displacement vector C is fourth order elasticity tensor or stiffness tensor ε is second order linearized finite strain tensor or tensor of small strain , ,i j kx x x – Cartesian coordinates. Strain and stress tensor are symmetric. This assumption leads us to a contracted form of the above equations, the Voigt notation. Equation (3) defines 6 different equations for only three components of the displacement vector ui, therefore some other conditions (4) need to be imposed over the small strain components in order to guarantee the existence of the displacement field (Mase, 1970). 2 22 2 0 0ij jmkm ik k m i j j i kx x x x x x x x ε εε ε ε ∂ ∂∂ ∂ ∇× ×∇ = ↔ + − − = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (4) Thus, combining (1), (2), (3) and (4) the most general form of the elastodynamic equation is obtained: 2 2 k ij j l u u C i x x t ρ ⎛ ⎞∂ ∂∂ =⎜ ⎟ ∂ ∂ ∂⎝ ⎠ (5) Note that the stiffness tensor reduces to a 6x6 symmetric matrix C i j k lC ij due to the symmetry of the stress and strain tensors. The condition (5) gives a simpler strain – displacement relation. The stiffness matrix can be obtained by combination of coordinate transformations corresponding to the material lattice symmetry properties. In the most general case the stiffness matrix has got 21 independent elastic constants. This is the monoclinic material. Hexagonal symmetry, for example, has got 5 different elastic constants. This system has got one of the axes of symmetry such that a rotation around this axis does not change the stiffness matrix. For homogeneous materials in which elastic constants and density are constant, (5) reduces to (6): 2 2 2 k i ij j l u u C x x t ρ ∂ ∂ = ∂ ∂ ∂ (6) First consider harmonic plane waves moving through homogeneous anisotropic medium. Then the solution must be sought in the next form: ( ) ( ) ( ) ( ), , i iik vt l xik vt l x i i iu x t Ue u x t U e − −− −= ↔ = (7) Where: U is the displacement amplitude vector, integration constant, i is an imaginary unit, 64 2 k π λ = is the wave number, λ is wavelength, v is phase velocity, t is time, l is unit vector along the propagation direction, x is position vector. Replacing the displacement vector (7) in the elastodynamic equation (6) and applying the Fourier transform, we get: 2 2 0ij jk kl i lL C L U Uk ρω − = (8) Where: ijL is Director cosines matrix of the propagation vector, 2 v kv π ω λ = = – angular frequency. The equation (8) can be rewritten as: ( )2 0ik ik kv Uδ ρΓ − = (9) Where: il ij jk klL C LΓ = ⋅ ⋅ – Christoffel matrix (10) The slowness 1 p v = is obtained by solving the next equation for the three values of v2. ( 2det 0ik ik vδ ρΓ − =) (11) Therefore, three types of plane harmonic waves along the propagation direction l are obtained, with three different phase velocities v(i). The corresponding displacement vectors are mutually orthogonal; they are also called polarization vectors. The vectors ( )i jU l and can be parallel, in that case the corresponding wave to that displacement direction is called pure longitudinal wave. Perpendicular directions of displacement correspond to pure shear waves. ( )i jU Body of the programming routine This routine plots the wave surfaces for all cases of anisotropy > restart: with(LinearAlgebra): unprotect(gamma,GAMMA): Stiffness matrices Stiffness matrix in the Voigt form > C:=Matrix(6,6,shape=symmetric): Now the user must choose the anisotropy case: Isotropic case (2 elastic constants) > C[1,1]:=lambda+2*mu: C[2,2]:=lambda+2*mu: C[3,3]:=lambda+2*mu: C[4,4]:=mu: C[5,5]:=mu: C[6,6]:=mu: C[1,2]:=lambda: C[1,3]:=lambda: C[2,3]:=lambda: evalm(C); Cubic anisotropy (3 elastic constants) > C[1,1]:=C11: C[2,2]:=C11: C[3,3]:=C11: C[4,4]:=C44: C[5,5]:=C44: C[6,6]:=C44: C[1,2]:=C12: C[1,3]:=C12: C[2,3]:=C12: evalm(C); Hexagonal anisotropy, VTI (5 elastic constants) > C[1,1]:=C11: C[2,2]:=C11: C[3,3]:=C33: C[4,4]:=C44: C[5,5]:=C44: C[6,6]:=C66: C[1,2]:=C11-2*C66: C[1,3]:=C13: C[2,3]:=C13: eval(C); Orthorhombic anisotropy (9 elastic constants) > C[1,1]:=C11: C[2,2]:=C22: C[3,3]:=C33: C[4,4]:=C44: C[5,5]:=C55: C[6,6]:=C66: C[1,2]:=C12: C[1,3]:=C13: C[2,3]:=C23: evalm(C); Monoclinic anisotropy (13 elastic constants) 65 Miguel Duarte et al. 66 > C[1,1]:=C11: C[1,2]:=C12: C[1,3]:=C13: C[1,4]:=C14: C[1,5]:=C15: C[1,6]:=C16: C[2,2]:=C22: C[2,3]:=C23: C[2,4]:=C24: C[2,5]:=C25: C[2,6]:=C26: C[3,3]:=C33: C[3,4]:=C34: C[3,5]:=C35: C[3,6]:=C36: C[4,4]:=C44: C[4,5]:=C45: C[4,6]:=C46: C[5,5]:=C55: C[5,6]:=C56: C[6,6]:=C66: evalm(C); Ray director cosines matrix > L:=Matrix(6,3): L[1,1]:=l[1]: L[2,2]:=l[2]: L[3,3]:=l[3]: L[4,2]:=l[3]: L[4,3]:=l[2]: L[5,1]:=l[3]: L[5,3]:=l[1]: L[6,1]:=l[2]: L[6,2]:=l[1]: evalm(L); Christoffel Matrix > GAMMA:=Transpose(L).C.L; Eigenvalues and Eigenvectors The elastic constants should be declared according to the chosen anisotropy case: consts:=C11=16.57*10^10,C12=6.39*10^10,C44 =7.956*10^10,rho=2332: (Cubic anisotropy) > consts:=C11=3.46E10,C13=1.07E10,C33 =2.84E10,C44=8.36E9,C66=1.26E10,rho=2500: (VTI elastic constants) >dir_l:=l[1]=sin(theta),l[2]=0,l[3]=cos(theta): >Ch:=CharacteristicMatrix(subs(dir_l,GAMMA), k); > d:=Determinant(Ch); > factor(d); > K:=[solve(d,k)]; > map(simplify,%,trig); > A:=subs(k=K[3],Ch); > A:=map(simplify,subs(k=K[3],Ch),trig); > solve({A[1,1]*u1+A[1,2]*u2=0},{u1}); > simplify(%,trig); > solve({A[2,1]*u1+A[2,2]*u2=0},{u1}); > simplify(%,trig); > Rank(A); > eg:=map(simplify,ColumnSpace(A)); > plot([subs(consts,dir_l,sqrt(rho/K[1])),subs(consts ,dir_l,sqrt(rho/K[2])),subs(consts,dir_l,sqrt(rho/K[ 3]))],theta=0..2*Pi,color=[red,blue,green],coords= polar); > plot([subs(consts,dir_l,sqrt(K[1]/rho)),subs(consts ,dir_l,sqrt(K[2]/rho)),subs(consts,dir_l,sqrt(K[3]/r ho))],theta=0..2*Pi,color=[red,blue,green],coords= polar); > phi_1:=theta+arctan(diff(sqrt(K[1]/rho),theta)/sqrt (K[1]/rho)); > v_1:=sqrt(sqrt(K[1]/rho)^2+diff(sqrt(K[1]/rho),the ta)^2); > phi_2:=theta+arctan(diff(sqrt(K[2]/rho),theta)/sqrt (K[2]/rho)): > v_2:=sqrt(sqrt(K[2]/rho)^2+diff(sqrt(K[2]/rho),the ta)^2): > phi_3:=theta+arctan(diff(sqrt(K[3]/rho),theta)/sqrt (K[3]/rho)): > v_3:=sqrt(sqrt(K[3]/rho)^2+diff(sqrt(K[3]/rho),the ta)^2): > plot({[subs(consts,dir_l,v_1),subs(consts,dir_l,phi _1),theta=0..2*Pi],[subs(consts,dir_l,v_2),subs(co nsts,dir_l,phi_2),theta=0..2*Pi],[subs(consts,dir_l, v_3),subs(consts,dir_l,phi_3),theta=0..2*Pi]},colo r=[red,blue,green],coords=polar); > plot({subs(consts,dir_l,sqrt(K[3]/rho)),[subs(const s,dir_l,v_3),subs(consts,dir_l,phi_3),theta=0..Pi/2] },theta=0..Pi/2,color=[red,blue],coords=polar); In the anterior routine the term VTI means Vertical and Transversally Isotropic media. Routine work examples This program plots the three wave surfaces for any anisotropic material: slowness surface, phase velocity surface and group velocity surface. The program was tested using the next stiffness matrix for the VTI case: 10 3 3.4 0.9 1.07 0 0 6 0.94 1.07 0 0 0 4 3.46 1.07 0 0 0 1.07 2.84 0 0 0 10 , 2500 0 0 0.836 0 0 0 0 0 0.836 0 0 0 0 0 0 1.26 kg C Pa m ρ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ = ⋅⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎝ ⎠ ij = Figure 1. Stiffness matrix for sandstone (Thomsen, 2000) CONCLUSIONS The theoretical study of anisotropic media is important for studying the effect of different type of media on wave propagation. Particularly the effect that fractures, typical of Colombian Foothills. We expect to continue our study of this concepts and the development of computational tools that simulate wave propagation in general Transversal Inhomogeneous media. Figure 2. Slowness surfaces for sandstone. Horizontal axis x3. Vertical axis x1. Slowness given in s/m. Figure 3. Slowness surfaces for sandstone. Horizontal axis x3. Vertical axis x1. Phase velocity given in m/s. ACKNOWLEDGMENTS The authors wish to thank the support received from the Colombian Institute of Petroleum (ICP), the research division of ECOPETROL, in particular the research unit UIN, which had let accessed laboratory data as well as information related to anisotropy. REFERENCES Sen, Mrinal. (2002) Lectures: “Seismic Wave Propagation in Anisotropic Media”. The University of Texas at Austin. Krishnaswamy, Sridhar. (2004) “Acousto-Optical Methods of Materials Characterization”. Northwestern University. Pšenčík, Ivan. (1994) “Introduction to Seismic Methods”. PPPG/UFBa, El Salvador. Mase, George E. (1970) “Theory and Problems of Continuum Mechanics”. Schaum’s Outline Series, McGraw-Hill Thomsen, Leon. (2000) “Weak elastic anisotropy”. Applied Seismic Anisotropy: Theory, Background, and Field Studies. Geophysics reprint series No. 20, SEG. 67