Tìm hiểu phương pháp ma trận đối với sóng rayleigh trong môi trường phân lớp

ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN KHOA TOÁN - CƠ - TIN HỌC ————oOo———— Trần Ngọc Trung TÌM HIỂU PHƯƠNG PHÁP MA TRẬN ĐỐI VỚI SÓNG RAYLEIGH TRONG MÔI TRƯỜNG PHÂN LỚP KHÓA LUẬN TỐT NGHIỆP ĐẠI HỌC CHÍNH QUY Ngành: Toán - Cơ Cán bộ hướng dẫn: TS. Trần Thanh Tuấn Hà Nội - 2012 LỜI CẢM ƠN Lời đầu tiên của khóa luận này em xin gửi lời cảm ơn sâu sắc tới thầy giáo hướng dẫn TS. Trần Thanh Tuấn. Thầy đã giao đề tài và tận tình hướng dẫn em trong quá trình hoàn thành khóa luận này. Em cũng xin bày tỏ lòng biết ơn chân thành tới nhóm Seminar tại bộ môn Cơ học do PGS. TS Phạm Chí Vĩnh chủ trì, cùng toàn thể các thầy cô giáo trong khoa Toán - Cơ - Tin học, Đại học Khoa Học Tự Nhiên, Đại Học Quốc Gia Hà Nội đã dạy bảo em tận tình trong suốt quá trình học tập tại khoa. Nhân dịp này em cũng xin được gửi lời cảm ơn chân thành tới gia đình, bạn bè đã luôn bên em, cổ vũ, động viên, giúp đỡ em trong suốt quá trình học tập và thực hiện khóa luận tốt nghiệp. Hà Nội, ngày 20 tháng 05 năm 2012 Sinh viên Trần Ngọc Trung Mục lục Lời mở đầu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Chương 1. Phương trình tán sắc của sóng mặt trong môi trường đa lớp . 6 1.1. Dạng ma trận của bài toán cho sóng Rayleigh . . . . . . . . . . . . . . . . . . . 6 1.2. Một số tính chất tổng quát của nghiệm . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3. Tỉ số H/V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Chương 2. Dạng tiệm cận của phương trình tán sắc theo bước sóng . 17 2.1. Dạng tiệm cận của bước sóng dài . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2. Dạng tiệm cận cho bước sóng ngắn . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Chương 3. Tính toán số . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Kết luận . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Tài liệu tham khảo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Phụ lục . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2 Lời mở đầu Sóng Rayleigh [1] đã được Rayleigh phát hiện vào năm 1885. Từ đó đến nay đã có rất nhiều nghiên cứu về sóng này vì các ứng dụng rộng rãi của nó trong các lĩnh vực khác nhau. Một trong những lĩnh vực quan trọng sử dụng sóng Rayleigh đó là sự truyền sóng động đất, trong đó sóng mặt Rayleigh trong rất nhiều trường hợp là thành phần chính sóng động đất bên cạnh sóng Love và sóng khối, đặc biệt là trong các trường hợp khi tâm chấn (động đất) là khá xa so với địa điểm được khảo sát. Bề mặt trái đất ban đầu có thể được coi như là một bán không gian, nhưng mô hình chính xác của nó phải là môt hình phân lớp trong đó có một số lớp với các tính chất vật liệu khác nhau đặt trên một bán không gian. Hai tính chất cơ bản quan trọng nhất của sóng Rayleigh là vận tốc sóng và tỷ số H/V, là tỷ số của chuyển dịch ngang (horizontal) và chuyển dịch theo phương thẳng đứng (vertical). Sóng Rayleigh truyền trong bán không gian đàn hồi đẳng hướng là sóng không tán sắc và chỉ có duy nhất một mặt sóng. Trong mô hình đơn giản này hai tính chất cơ bản trên của sóng Rayleigh đã được nghiên cứu kỹ lưỡng (xem bài báo của Malischewsky (2004) [6]). Đối với mô hình phức tạp hơn, một lớp đặt trên bán không gian, trong trường hợp vật liệu là đàn hồi đẳng hướng, thì hai tính chất cơ bản trên cũng đã được nghiên cứu kỹ một cách giải tích trong nhiều công trình, ví dụ như trong bài báo của Trần Thanh Tuấn (2011) [11]. Nói chung phương pháp được dùng trong các mô hình đơn giản này là biểu diễn các đại lượng ứng suất và biến dạng của lớp và bán không gian phụ thuộc vào các tham số vật liệu và số sóng, sau đó sử dụng các điều kiện biên để nhận được một hệ phương trình thuần nhất. Phương trình tán sắc của sóng Rayleigh sau đó nhận được bằng cách cho định thức của hệ phương trình thuần nhất này bằng không để nhận được nghiệm không tầm thường. Phương pháp này có thể MỤC LỤC cho ta phương trình tán sắc của sóng Rayleigh dưới dạng hiển, thuận tiện cho việc nghiên cứu giải tích cũng như là các tính toán số. Tuy nhiên, phương pháp sẽ trở nên rất cồng kềnh khi được dùng để nghiên cứu mô hình phân lớp, khi số lớp là nhiều hơn hai. Một phương pháp thay thế để khảo sát mô hình phân lớp chính là phương pháp ma trận chuyển. Phương pháp này được đề xuất bởi Thomson (1950) [9]. Haskell (1953) [4] đã phát triển phương pháp này đối với môi trường đàn hồi đẳng hướng và Stuart Crampin (1970) [2] đã phát triển phương pháp này cho môi trường đàn hồi bất đẳng hướng. Phương pháp này sẽ cho ta phương trình tán sắc của sóng Rayleigh dưới dạng ẩn. Mặc dù khó có thể sử dụng phương trình tán sắc dạng ẩn này để nghiên cứu một cách giải tích các tính chất của sóng Rayleigh, nhưng nó được dùng một cách rộng rãi trong việc lập các chương trình tính toán số để khảo sát số các tính chất của sóng Rayleigh trong mô hình phân lớp này. Một trong các chương trình sử dụng phương pháp ma trận chuyển này là chương trình của Herrmann (1994) [5]. Chương trình này được viết bởi ngôn ngữ lập trình FORTRAN dưới dạng các gói lệnh và có thể sử dụng các ngôn ngữ khác, ví dụ như Matlab, để chạy chúng. Ưu điểm của chương trình là chạy ổn định và nhanh chóng. Tuy nhiên nhược điểm của nó là người sử dụng không thể can thiệp trực tiếp vào các code lệnh để thay đổi chương trình để phục vụ mục đích của mình. Một ví dụ là chương trình của Herrman tính toán vận tốc sóng và tỷ số H/V trên một miền tần số cho trước và nó chia miền tần số này thành một số khoảng rời rạc bằng nhau (nói chung là 29 hoặc 210 khoảng). Trong các tính toán thông thường thì số lượng chia rất lớn này nói chung là đủ để đáp ứng các yêu cầu của người dùng. Tuy nhiên với một số mục đích riêng biệt khác, như tính toán xung quanh điểm osculation point, là điểm hai mode của đường cong tán sắc gặp nhau, thì việc chia số khoảng lớn như trên vẫn chưa đủ. Hoặc như khi cần vẽ đường cong tán sắc và đường cong tỷ số H/V một cách rất mịn thì chúng ta cần phải can thiệp vào code của chương trình, và việc mà người dùng khó có thể làm khi sử dụng chương trình của Herrmann. Hoặc quan trọng nhất là giải số tìm các tần số của điểm không (điểm có H/V bằng không) và điểm singularity (là điểm có tỷ số H/V bằng vô cùng). Với lý do trên, việc thành lập code chương trình cho phương pháp ma trận chuyển là cần thiết. Vì vậy, mục tiêu của khóa luận tốt nghiệp này là đi tìm hiểu phương pháp ma trận chuyển và sử dụng ngôn ngữ lập trình Matlab để viết code cho phương pháp này. Nội dung chính của khóa luận là trình bày lại các kết quả 4 MỤC LỤC đã được đăng trong bài báo của Thomson (1950) [9] và của Haskell (1953) [4]. Thông qua các kết quả trong hai bài báo trên, công thức của tỷ số H/V đối với môi trường phân lớp cũng được dễ dàng tìm ra. Dựa trên các kết quả và công thức này, tác giả đã sử dụng ngôn ngữ lập trình Matlab để viết chương trình tính toán vận tốc sóng Rayleigh và tỷ số H/V đối với mô hình phân lớp có số lớp tùy ý. Cuối cùng, một số kết quả tính toán số sử dụng chương trình này cũng được trình bày. Khóa luận bao gồm ba chương : • Chương 1: Tìm phương trình tán sắc tổng quát và tỉ số H/V. • Chương 2: Tìm dạng tiệm cận của phương trình tán sắc theo bước sóng ngắn và bước sóng dài. • Chương 3: Áp dụng tính toán số. Do thời gian thực hiện khóa luận không nhiều, kiến thức còn hạn chế nên khi làm khóa luận không tránh khỏi những hạn chế và sai sót. Em mong nhận được sự góp ý và những ý kiến phản biện của quý thầy cô và bạn đọc. Xin chân thành cảm ơn! 5 Chương 1 Phương trình tán sắc của sóng mặt trong môi trường đa lớp 1.1. Dạng ma trận của bài toán cho sóng Rayleigh Ta xét sóng mặt có tần số góc p và vận tốc theo phương ngang c lan truyền theo phương ngang trong không gian với n lớp song song đồng nhất, đẳng hướng. +x +z (0) (1) (2) Phương truyền sóng 1 2 3 n-1 n (n-1) Hình 1 : Chiều của các trục tọa độ và cách đánh số của các lớp và các mặt phân cách. Tất cả các lớp giả sử đều là môi trường chất rắn. Trục x song song với các lớp và có chiều dương hướng theo chiều của phương truyền sóng. Trục z có chiều dương hướng vào trong môi trường. Các lớp khác nhau và các mặt phân cách được đánh số bắt đầu từ mặt tự do như hình 1. Ta chú ý dạng của sóng Rayleigh, tức là ở đây không có chuyển dịch theo phương y và biên độ giảm theo hàm mũ trong phương +z trong lớp bán không gian. 6 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP Xét lớp thứ m, ta đặt u,v,w = thành phần chuyển dịch theo chiều x,y và z u˙, v˙, w˙ = thành phần vận tốc theo chiều x,y và z ∆m = giãn nở khối ωm = sự quay αm = [(λm+2µm)/ρm]1/2 = vận tốc của sóng dọc βm = [µm/ρm]1/2 = vận tốc của sóng ngang λm = νE (1+ν)(1−2ν) E = môđun đàn hồi Young ν = hệ số Poisson λm,µm = hệ số đàn hồi Lame ρm = mật độ khối lượng dm = độ dày lớp thứ m p = tần số góc k = p/c = 2pi/bước sóng (theo phương ngang) rαm = { +[(c/αm)2−1]1/2,nếu c > αm −i[1− (c/αm)2]1/2,nếu c < αm rβm = { +[(c/βm)2−1]1/2,nếu c > βm −i[1− (c/βm)2]1/2,nếu c < βm γm = 2(βm/c)2 σ = ứng suất pháp τ = ứng suất tiếp. Giãn nở khối và sự quay được xác định qua phương trình sau ∆m = ∂u ∂x + ∂v ∂y + ∂w ∂ z , ωmx = 1 2 ( ∂w ∂y − ∂v ∂ z ) , ωmy = 1 2 ( ∂u ∂ z − ∂w ∂x ) , ωmz = 1 2 ( ∂v ∂x − ∂u ∂y ) . 7 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP Đạo hàm ∆m theo x,y và z ta được ba phương trình ∇2u = ∂∆m ∂x −2∂ωmz ∂y +2 ∂ωmy ∂ z , ∇2v = ∂∆m ∂y −2∂ωmx ∂ z +2 ∂ωmz ∂x , (1.1) ∇2w = ∂∆m ∂ z −2∂ωmy ∂x +2 ∂ωmx ∂y . Nghiệm tổng quát cho chuyển dịch theo ∆m và ωm trong đó chúng thỏa mãn phương trình sóng cho giãn nở khối và sự quay có dạng hàm điều hòa phụ thuộc thời gian là u = − ( αm p )2∂∆m ∂x +2 ( βm p )2∂ωmz ∂y −2 ( βm p )2∂ωmy ∂ z , v = − ( αm p )2∂∆m ∂y +2 ( βm p )2∂ωmx ∂ z −2 ( βm p )2∂ωmz ∂x , (1.2) w = − ( αm p )2∂∆m ∂ z +2 ( βm p )2∂ωmy ∂x −2 ( βm p )2∂ωmx ∂y . Ví dụ như, từ phương trình thứ nhất của hệ (1.2) trên ta có ∇2u = − ( αm p )2 ∂ ∂x ( ∇2∆m ) +2 ( βm p )2 ∂ ∂y ( ∇2ωmz )−2(βm p )2 ∂ ∂ z ( ∇2ωmy ) , mà khi so sánh với phương trình thứ nhất của (1.1) cho ta phương trình chuyển động của sóng ( ∇2 + ( p αm )2) ∆m = 0,( ∇2 + ( p βm )2) ωmz = 0,( ∇2 + ( p βm )2) ωmy = 0. Chọn hệ tọa độ đề Đề Các như hình 1, ωmx,ωmz sẽ bằng không. Xét lớp thứ m 8 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP thì giãn nở khối và sự quay mà thỏa mãn phương trình sóng là ∆m = (∂u/∂x)+(∂w/∂ z) = exp[i(pt− kx)][∆′m exp(−ikrαmz)+∆′′m exp(ikrαmz)], (1.3) ωm = (1/2)[(∂u/∂ z)− (∂w/∂x)] = exp[i(pt− kx)][ω ′m exp(−ikrβmz)+ω ′′m exp(ikrβmz)], (1.4) trong đó ∆′m,∆′′m,ω ′m và ω ′′m là các hằng số. Với quy ước về ký hiệu như trên, số hạng trong biểu diễn ∆′m biểu thị cho một sóng mặt mà hướng chuyển động của nó hợp với trục +z một góc cot−1 rαm với rαm là thực, và một sóng lan truyền theo phương +x với biên độ giảm dần theo hàmmũ theo chiều của +z khi rαm là số ảo. Tương tự, số hạng trong ∆′′m biểu thị cho một sóng mặt mà hướng chuyển động của nó hợp với trục −z một góc như trên khi rαm là số thực và một sóng lan truyền trong +x với biên độ giảm dần theo hàm mũ theo chiều của trục z khi rαm là số ảo. Áp dụng cách đánh dấu tương tự như trên, ta có được dạng ω ′m và ω ′′m với rβm thay thế cho rαm. Thành phần tương ứng của ứng suất và chuyển vị tương ứng với hàm giãn nở khối và sự quay được biểu diễn bởi (1.3) và (1.4) là, u = −(αm/p)2(∂∆m/∂x)−2(βm/p)2(∂ωm/∂ z), (1.5) w = −(αm/p)2(∂∆m/∂ z)+2(βm/p)2(∂ωm/∂x), (1.6) σ = λm∆m+2µm(∂u/∂x) = ρm[α2m∆m+2β 2m{(αm/p)2(∂ 2∆m/∂x2)+2(βm/p)(∂ 2ωm/∂x∂ z)}], (1.7) τ = µm(∂u/∂ z+∂w/∂x) = 2ρmβ 2m[− (αm/p)2(∂ 2∆m/∂x∂ z)+(βm/p)2{(∂ 2ωm/∂x2)− (∂ 2ωm/∂ z2)}]. (1.8) Các điều kiện biên thỏa mãn tại mặt phân cách giữa hai lớp là bốn đại lượng thỏa mãn điều kiện liên tục. Chuyển vị chắc chắn là liên tục nếu các thành phần vận tốc tương ứng u˙ và w˙ là liên tục, và vì c là như nhau trong tất cả các lớp nên chúng ta có thể xét các đại lượng không thứ nguyên u˙/c và w˙/c là liên tục. Thay các biểu diễn của (1.3) và (1.4) vào phương trình (1.5) đến (1.8), ta 9 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP có u =ik(αm/p)2 exp[i(pt− kx)][∆′m exp(−ikrαmz)+∆′′m exp(ikrαmz)] +2ikrβm(βm/p)2 exp[i(pt− kx)][ω ′m exp(−ikrβmz)−ω ′′m exp(ikrβmz)], (1.9) w =ikrαm(αm/p)2 exp[i(pt− kx)][∆′m exp(−ikrαmz)−∆′′m exp(ikrαmz)] +2ik(βm/p)2 exp[i(pt− kx)][ω ′m exp(−ikrβmz)+ω ′′m exp(ikrβmz)]. (1.10) Thành phần vận tốc tương ứng được biểu diễn như sau u˙ =− (α2m/c)exp[i(pt− kx)][∆′m exp(−ikrαmz)+∆′′m exp(ikrαmz)] −2rβm ( β 2m/c ) exp[i(pt− kx)][ω ′m exp(−ikrβmz)−ω ′′m exp(ikrβmz)], (1.11) w˙ =− rαm ( α2m/c ) exp[i(pt− kx)][∆′m exp(−ikrαmz)−∆′′m exp(ikrαmz)] −2(β 2m/c)exp[i(pt− kx)][ω ′m exp(−ikrβmz)+ω ′′m exp(ikrβmz)]. (1.12) Biểu diễn các hàm mũ của ikrz theo dạng lượng giác, ta tìm được u˙/c =− (αm/c)2[(∆′m+∆′′m)coskrαmz− i(∆′m−∆′′m)sinkrαmz] − γmrβm[(ω ′m−ω ′′m)coskrβmz− i(ω ′m+ω ′′m)sinkrβmz], (1.13) w˙/c =− (αm/c)2rαm[−i(∆′m+∆′′m)sinkrαmz+(∆′m−∆′′m)coskrαmz] + γm[−i(ω ′m−ω ′′m)sinkrβmz+(ω ′m+ω ′′m)coskrβmz], (1.14) σ =−ρmα2m(γm−1)[(∆′m+∆′′m)coskrαmz− i(∆′m−∆′′m)sinkrαmz] −ρmc2γ2mrβm[(ω ′m−ω ′′m)coskrβmz− i(ω ′m+ω ′′m)sinkrβmz], (1.15) τ = ρmα2mγmrβm[−i(∆′m+∆′′m)sinkrαmz+(∆′m−∆′′m)coskrαmz] −ρmc2γm(γm−1)[−i(ω ′m−ω ′′m)sinkrβmz+(ω ′m+ω ′′m)coskrβmz]. (1.16) Khi bất kỳ các giá trị rαm,rβm là ảo thì hàm lượng giác tương ứng sẽ được hiểu là hàm lượng giác hyperbolic. Đặt gốc tọa độ của z tại mặt phân cách thứ (m− 1) mối quan hệ tuyến tính giữa các giá trị của u˙/c, w˙/c,σ và τ tại mặt phân cách thứ (m−1) và các hằng số (∆′m +∆′′m),(∆′m −∆′′m) và (ω ′m +ω ′′m) có thể được biểu diễn bởi phép 10 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP biến đổi sau (u˙m−1/c, w˙m−1/c,σm−1,τm−1) = Em(∆′m+∆′′m,∆′m−∆′′m,ω ′m−ω ′′m,ω ′m+ω ′′m), (1.17) trong đó Em là ma trận Em =  −(αm/c)2 0 −γmrβm 0 0 −(αm/c)2rαm 0 γm −ρmα2m(γm−1) 0 −ρmc2γmrβm 0 0 ρmα2mγmrαm 0 −ρmc2γm(γm−1)  . (1.18) Đặt z = dm trong phương trình (1.13) đến (1.16) ta thu được giá trị của u˙/c, v.v... tại mặt phân cách m theo (∆′m+∆′′m) v.v... (u˙m/c, w˙m/c,σm,τm) = Dm(∆′m+∆′′m,∆′m−∆′′m,ω ′m−ω ′′m,ω ′m+ω ′′m), (1.19) trong đó Dm là ma trận Dm =  −(γm/c)2 cosPm i(αm/c)2 sinPm −γmrβm cosQm iγmrβm sinQm i(αm/c)2 rαm sinPm −(αm/c)2 rαm cosPm −iγm sinQm γm cosQm −ρmα2m (γm−1)cosPm iρmα2m (γm−1)sinPm −ρmc2γ2mrβm cosQm iρmc2γ2mrβm sinQm −iρmα2mγmrαm sinPm ρmα2mγmrαm cosPm iρmc2γm (γm−1)sinQm −ρmc2γm (γm−1)cosQm  , (1.20) với Pm = krαmdm và Qm = krβmdm. Các hằng số ∆′m+∆′′m, v.v... có thể được khử giữa hai phương trình (1.17) và (1.19), cho ta một mối quan hệ tuyến tính giữa các giá trị của u˙/c, w˙/c,σ và τ tại đỉnh và đáy của lớp thứ m mà có thể được biểu diễn một cách hình thức bởi phương trình, (u˙m/c, w˙m/c,σm,τm) = DmE−1m (u˙m−1/c, w˙m−1/c,σm−1,τm−1). (1.21) trong đó E−1m là ma trận nghịch đảo của Em và có dạng E−1m =  −2(βm/αm)2 0 ( ρmα2m )−1 0 0 c2 (γm−1)/α2mrαm 0 ( ρmα2mrαm )−1 (γm−1)/γmrβm 0 − ( ρmc2γmrβm )−1 0 0 1 0 ( ρmc2γm )−1  . (1.22) Từ các phương trình (1.20) và (1.22), các phần tử của ma trận tích am = DmE−1m 11 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP có thể được tính như sau (am)11 = γm cosPm− (γm−1)cosQm (am)12 = i [ (γm−1)r−1αm sinPm+ γmrβm sinQm ] (am)13 = − ( ρmc2 )−1 (cosPm− cosQm) (am)14 = i ( ρmc2 )−1 ( r−1αm sinPm+ rβm sinQm ) (am)21 = −i [ γmrαm sinPm+(γm−1)r−1βm sinQm ] (am)22 = −(γm−1)cosPm+ γm cosQm (am)23 = i ( ρmc2 )−1( rαm sinPm+ r−1βm sinQm ) (am)24 = (am)13 (am)31 = ρmc 2γm (γm−1)(cosPm− cosQm) (am)32 = iρmc 2 [ (γm−1)2 r−1αm sinPm+ γ2mrβm sinQm ] (am)33 = (am)22 (am)34 = (am)12 (am)41 = iρmc 2 [ γ2mrαm sinPm+(γm−1)2 r−1βm sinQm ] (am)42 = (am)31 (am)43 = (am)21 (am)44 = (am)11 . Bây giờ xét điều kiện biên với yêu cầu là các giá trị của u˙/c, w˙/c,σ và τ được tính tại đỉnh của lớp thứ m phải bằng các giá trị tương ứng được tính tại đáy của lớp thứ (m−1). Điều đó có nghĩa là ta có thể viết (u˙m/c, w˙m/c,σm,τm) = amam−1(u˙m−2/c, w˙m−2/c,σm−2,τm−2). (1.23) Trong bài báo của Thomson (1950) [9], đại lượng τ/2µ được thay cho τ như biến thứ tư của ma trận am. Đây chỉ là sự thay đổi trong ký hiệu với bất kỳ lớp nào có liên quan nhưng µ nói chung sẽ khác nhau trong các lớp, và đó là ứng suất trượt τ , không phải là biến dạng trượt τ/µ , là đại lượng liên tục qua các mặt phân cách. Quá trình lặp được chỉ ra bởi phương trình (1.23) do đó yêu cầu τ , hoặc một hằng số nhân với τ , không phải là τ/2µ , được coi như biến thứ tư. Dạng ma trận a của Thomson khi đó sẽ được hiệu chỉnh lại bằng việc nhân hàng thứ tư với 2µ (= 2G trong kí hiệu của Thomson) và cột thứ tư được nhân với (2µ)−1. 12 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP Bằng việc lặp lại phương trình (1.23) ta có, (u˙n−1/c, w˙n−1/c,σn−1,τn−1) = an−1an−2...a1(u˙0/c, w˙0/c,σ0,τ0), (1.24) và sử dụng phương trình nghịch đảo (1.17) cho lớp thứ n, ta có (∆′n+∆ ′′ m,∆ ′ n−∆′′n,ω ′n−ω ′′n ,ω ′n+ω ′′n ) = E−1n an−1an−2...a1(u˙0/c, w˙0/c,σ0,τ0). (1.25) Đến bước này mọi bước thực hiện đều là tổng quát, và phương trình (1.25) có thể được áp dụng đối với sóng mặt hoặc sóng truyền qua môi trường nhiều lớp. Trường hợp cụ thể mà chúng ta quan tâm trong đó không có ứng suất trên mặt tự do nên σ0 = τ0 = 0 và không có nguồn kích động nào ở vô cùng nên ∆′′n = ω ′′n = 0. Kí hiệu J là ma trận tích E−1n an−1an−2...a1, phương trình (1.25) trở thành, ( ∆′n,∆ ′ n,ω ′n,ω ′n ) = J (u˙/c, w˙/c,0,0) , hoặc, cụ thể là ∆′n = J11u˙0/c+ J12w˙0/c, ∆′n = J21u˙0/c+ J22w˙0/c, ω ′n = J31u˙0/c+ J32w˙0/c, ω ′n = J41u˙0/c+ J42w˙0/c. (1.26) Khử ∆′n và ω ′n ta có, u˙0 w˙0 = J22− J12 J11− J21 = J42− J32 J31− J41 . (1.27) Bởi vì những phần tử của ma trận J là các hàm của các tham số c và k, phương trình (1.27) cho ta mối liên hệ ẩn giữa c và k, đó chính là phương trình tán sắc của vận tốc. 1.2. Một số tính chất tổng quát của nghiệm Đặt A= an−1an−2...a1 và sử dụng phương trình (1.22) đối với E−1n , phương trình (1.27) có thể viết dưới dạng − (u˙0/w˙0) = K/L = M/N, (1.28) 13 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP trong đó K = γnrαnA12 +(γn−1)A22− rαnA32/ρnc2 +A42/ρnc2, L = γnrαnA11 +(γn−1)A21− rαnA31/ρnc2 +A41/ρnc2, M = −(γn−1)A12 + γnrβnA22 +A32/ρnc2 + rβnA42/ρnc2, N = −(γn−1)A11 + γnrβnA21 +A31/ρnc2 + rβnA41/ρnc2. (1.29) Trong trường hợp hai lớp, phương trình (1.28) và (1.29) có thể được biểu diễn dưới dạng tương đương với biểu diễn trước đó của Sazawa, Lee (1938) [7] và các tác giả khác. Trong các biểu diễn cho các phần tử của ma trận am, chúng ta thấy rằng các đại lượng sinPm,sinQm, rαm và rβm có thể là thực hoặc ảo phụ thuộc vào giá trị của c, chỉ xảy ra trong tích r±1αm sinPm và r ±1 βm sinQm. Vì sinPm là thực hoặc ảo khi rαm là thực hoặc ảo, và sinQm là tương tự đối với rβm, nên các tích trên luôn là thực đối với giá trị thực của c. Với các tính chất thực hoặc ảo của phần tử thuộc mà trận am trên thì am có dạng am =  R I R I I R I R R I R I I R I R  . trong đó R biểu thị cho các đại lượng thực và I biểu thị cho các đại lượng ảo. Tích của hai ma trận bất kỳ của dạng này cũng là một ma trận cùng dạng. Vì thế các phần tử của A trong phương trình (1.29) có dạng: A11,A22,A31,và A42 là thực; A12,A21,A32,và A41 là ảo. Với việc định nghĩa sóng "mặt" là sóng có biên độ giảm khi giá trị của z tăng lên, nghĩa là rαn và rβn trong trường hợp này là ảo. Nghĩa là ta chỉ xét trường hợp giá trị của c < βn. Khi đó đối với phương trình (1.29), tất cả các số hạng của K và N là thực và tất cả các số hạng của L và M là ảo. Do đó tỉ số u˙0/w˙0 sẽ luôn là số ảo, nghĩa là có sự lệch pha 900 giữa các thành phần theo trục ngang và trục dọc tại mặt tự do. Do đó chuyển động của các phần tử có dạng hình ellipse với các trục có phương ngang và phương thẳng đứng. Tuy nhiên độ lệch pha có thể có dấu bất kỳ và do đó chuyển động ellipse của phần tử không nhất thiết là ngược chiều kim đồng hồ đối với phương truyền sóng, là chuyển động của phần tử đối với sóng Rayleigh trong mô hình bán không gian. 14 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP Nếu môi trường chúng ta đang xét là tán sắc, thì sẽ cần thiết phải xét giá trị phức của c và k. Trong trường hợp này, tỉ số u˙0/w˙0 sẽ không nhất thiết là thuần ảo, nghĩa là độ lệch pha có thể khác ±900 thường xảy ra và các trục của ellipse có thể nghiêng một góc với trục. Do đó môi trường đàn hồi không hoàn hảo có thể dẫn đến độ nghiêng của trục, là trường hợp hay được quan sát trong trường hợp các vụ nổ được kích hoạt để tạo ra sóng mặt trên lớp không chặt. Rõ ràng là nếu hai lớp liền kề có tính chất vật lý giống nhau, chúng phải được coi là tương đương với một lớp có bề dày bằng tổng bề dày của hai lớp. Do đó, nếu ta gọi am (d) là ma trận am được tính cho lớp có độ dày d ta phải có am (d1)am (d2) = am (d1 +d2) (1.30) Mối liên hệ này có thể dễ dàng kiểm tra bởi việc nhân trực tiếp các ma trận. Và bởi vì k chỉ xuất hiện trong am dưới dạng kdm, phương trình (1.30) cho ta, am (k1)am (k2) = am (k1 + k2) , (1.31) trong đó k1 và k2 là giá trị bất kỳ của k và am được tính đối với c và dm cố định. 1.3. Tỉ số H/V Tỷ số H/V là tỷ số của thành phần chuyển vị theo phương ngang và phương thẳng đứng của chất điểm trên bề mặt của lớp trên cùng. Nghĩa là : χ = u(0) w(0) . Ta có, u và w được xác định như sau u =ik(αm/p)2 exp[i(pt− kx)][∆′m exp(−ikrαmz)+∆′′m exp(ikrαmz)] +2ikrβm(βm/p)2 exp[i(pt− kx)][ω ′m exp(−ikrβmz)−ω ′′m exp(ikrβmz)], w =ikrαm(αm/p)2 exp[i(pt− kx)][∆′m exp(−ikrαmz)−∆′′m exp(ikrαmz)] +2ik(βm/p)2 exp[i(pt− kx)][ω ′m exp(−ikrβmz)+ω ′′m exp(ikrβmz)]. Dễ dàng thấy được u và w chỉ phụ thuộc vào thời gian t qua nhân tử exp[i(pt− kx)], do đó u w = u˙ w˙ , với u˙, w˙ được biểu diễn theo công thức (1.11) và (1.12) 15 CHƯƠNG 1. PHƯƠNG TRÌNH TÁN SẮC CỦA SÓNG MẶT TRONG MÔI TRƯỜNG ĐA LỚP Ta lại có u˙0 w˙0 = −K L = −M N . trong đó K, L, M, N được biểu diễn như trong (1.29). Do đó công thức xác định tỉ số H/V là χ = u(0) w(0) = −K L = −M N . 16 Chương 2 Dạng tiệm cận của phương trình tán sắc theo bước sóng 2.1. Dạng tiệm cận của bước sóng dài Khi bước sóng rất lớn thì kdm → 0 và tất cả các ma trận am sẽ tiến về ma trận đơn vị. Do đó, Jm → E−1m và phương trình (1.27) trở thành u˙0/w˙0 = −(γn−1)/γnrαn = γnrβn/(γn−1) , (2.1) hoặc (γn−1)2 + γ2n rαnrβn = 0. (2.2) Đây chính là phương trình vận tốc sóng cho sóng Rayleigh của bán không gian. Nếu ta khai triển dạng của ma trận am theo lũy thừa của k và bỏ qua lũy thừa cao hơn bậc một thì ma trận A có dạng A→  1 ik n−1 ∑ 1 dm 0 ik n−1 ∑ 1 dm/ρmβ 2m ik n−1 ∑ 1 dm [ 1−2(βm/αm)2 ] 1 ik n−1 ∑ 1 dm/ρmα2m 0 0 ikc2 n−1 ∑ 1 dmρm 1 ik n−1 ∑ 1 dm ikc2 n−1 ∑ 1 dmρm [ 1−2γm +2γm (βm/αm)2 ] 0 ik n−1 ∑ 1 dm [ 1−2(βm/αm)2 ] 1  . (2.3) Do đó, đối với bước sóng rất dài thì k2 ( n−1 ∑ 1 dm )2 có thể được bỏ qua. Đối với xấp xỉ ở bậc này thì các đại lượng K,L,M và N là các hàm tuyến tính của k. Do đó, đối với một giá trị của c thì phương trình (1.28) là một phương trình bậc hai của k và có thể giải một cách tường minh. 17 CHƯƠNG 2. DẠNG TIỆM CẬN CỦA PHƯƠNG TRÌNH TÁN SẮC THEO BƯỚC SÓNG 2.2. Dạng tiệm cận cho bước sóng ngắn Sezawa và Kanai (1938) [7] đã chỉ ra rằng trong trường hợp hai lớp, dạng tiệm cận đối với tần số sóng lớn của phương trình vận tốc có thể được nhóm dưới dạng tích của các thừa số. Một trong những thừa số này có một nghiệm tương ứng với vận tốc sóng Rayleigh trên mặt tự do của lớp đầu tiên; Nghiệm còn lại là biểu diễn của Stoneley (1924) [8] cho vận tốc sóng truyền trên các mặt phân cách giữa hai lớp. Nghiệm thứ hai có thể có hoặc có thể không có nghiệm thực, phụ thuộc vào mối liên hệ giữa ρ,α,β trong hai lớp. Theo lý thuyết vật lý thì trong trường hợp đa lớp, phương trình tán sắc luôn có thể được nhóm dưới dạng tích các thừa số khi tần số là đủ lớn. Những thừa số này tương ứng với sóng Rayleigh trên mặt tự do và sóng Stoneley ở mỗi mặt phân cách. Để minh họa điều này, sẽ thuận lợi hơn khi ta viết ma trận J dạng J = bn−1bn−2...b1E−11 , (2.4) trong đó bm = E−1m+1Dm. (2.5) Nghĩa là thay vì nhóm các số hạng của ma trận J theo từng lớp, chúng ta sẽ nhóm chúng theo các mặt phân cách. Bây giờ giả sử rằng c < βn−1, sao cho Pn−1 và Qn−1 là các số ảo và các hàm sin,cos được hiểu là các hàm hyperbolic. Khi đó, đối với giá trị của kdn−1 lớn, sinPn−1 →−icosPn−1 và sinQn−1 →−icosQn−1. Trong giới hạn này, phần tử của bn−1 tiến tới các giá trị sau : (bn−1)11 = −(bn−1)12 = (αn−1/αn)2{γn− (γn−1−1)(ρn−1/ρn)cosPn−1, (bn−1)13 = −(bn−1)14 = (c/αn)2 γn−1rβ (n−1){γn− γn−1 (ρn−1/ρn)}cosQn−1, (bn−1)21 = −(bn−1)22 = (αn−1/αn)2 ( rα(n−1)/rαn ){(γn−1)− γn−1 (ρn−1/ρn)}cosPn−1, (bn−1)23 = −(bn−1)24 = −(c/αn)2 (γn−1/rαn){(γn−1)− (γn−1−1)(ρn−1/ρn)}cosQn−1, 18 CHƯƠNG 2. DẠNG TIỆM CẬN CỦA PHƯƠNG TRÌNH TÁN SẮC THEO BƯỚC SÓNG (bn−1)31 = −(bn−1)32 = −(αn−1/c)2 ( γnrβn )−1{(γn−1)− (γn−1−1)(ρn−1/ρn)}cosPn−1, (bn−1)33 = −(bn−1)34 = ( γn−1rβ (n−1)/γnrβn ){(γn−1)− γn−1 (ρn−1/ρn)}cosQn−1, (bn−1)41 = −(bn−1)42 = (αn−1/c)2 ( rα(n−1)/γn ){γn− γn−1 (ρn−1/ρn)}cosPn−1, (bn−1)43 = −(bn−1)44 = (γn−1/γn){γn− (γn−1−1)(ρn−1/ρn)}cosQn−1, Nếu chúng ta đặt Jn−1 = bn−2bn−3...b1E−11 , thì bởi vì (bn−1) j1 =−(bn−1) j2 và (bn−1) j3 = −(bn−1) j4 đối với tần số cao, ta có thể viết J22− J12 = [(bn−1)11− (bn−1)21] [(Jn−1)22− (Jn−1)12] +[(bn−1)21− (bn−1)23] [(Jn−1)42− (Jn−1)32] , J11− J21 = [(bn−1)11− (bn−1)21] [(Jn−1)11− (Jn−1)21] +[(bn−1)21− (bn−1)23] [(Jn−1)31− (Jn−1)41] , J42− J32 = [(bn−1)31− (bn−1)41] [(Jn−1)22− (Jn−1)12] +[(bn−1)33− (bn−1)43] [(Jn−1)42− (Jn−1)32] , J31− J41 = [(bn−1)31− (bn−1)41] [(Jn−1)11− (Jn−1)21] +[(bn−1)33− (bn−1)43] [(Jn−1)31− (Jn−1)41] . Bằng cách đặt  K′ = (Jn−1)22− (Jn−1)12 , L′ = (Jn−1)11− (Jn−1)21 , M′ = (Jn−1)42− (Jn−1)32 , N′ = (Jn−1)31− (Jn−1)41 , (2.6)  R = (bn−1)11− (bn−1)21 , S = (bn−1)13− (bn−1)23 , T = (bn−1)31− (bn−1)41 , U = (bn−1)33− (bn−1)43 . (2.7) thì mối liên hệ (1.27) giữa các phần tử của J có thể viết dưới dạng RK′+SM′ RL′+SN′ = TK′+UM′ TL′+UN′ . 19 CHƯƠNG 2. DẠNG TIỆM CẬN CỦA PHƯƠNG TRÌNH TÁN SẮC THEO BƯỚC SÓNG Nhân chéo và rút gọn phương trình trên dẫn về dạng tích như sau (RU −ST )(K′N′−L′M′)= 0. (2.8) Cho thừa số đầu tiên bằng không và sử dụng giá trị của các phần tử của bn−1 cho ở trên, thừa số chung cosPn−1 cosQn−1 có thể được bỏ đi và biểu thức nhận được có thể được biểu diễn như sau ρ2n [ (γn−1)2 + γ2n rαnrβn ][ 1+ rα(n−1)rβ (n−1) ]−2ρnρn−1 [(γn−1)+ γnrαnrβn] [ (γn−1−1)+ γn−1rα(n−1)rβ (n−1)]+ρnρn−1[rα(n−1)rβn+rαnrβ (n−1)] +ρ2n−1[(γn−1−1)2 + γ2n−1rα(n−1)rβ (n−1)][1+ rαnrβn] = 0, (2.9) là phương trình tương đương với phương trình Stoneley đối với mặt phân cách thứ (n−1). Thừa số thứ hai của phương trình (2.8), khi được cho bằng không, thì tương đương với biểu thức gốc ban đầu của phương trình tán sắc, ngoại trừ việc nó là phương trình cho (n−1) lớp thay vì là cho n lớp. Lặp lại quá trình dẫn đến kết quả cho mỗi phần tử tương đương với sóng Rayleigh trên một trong những mặt phân cách và thừa số cuối cùng dạng phương trình (2.2) cho lớp đầu tiên, do đó biểu diễn được sóng Rayleigh trên mặt tự do của lớp này. Cũng phải chú ý rằng việc phân tách thành nhân tử của phương trình tán sắc đối với tần số lớn theo quá trình trên chỉ có thể được thực hiện khi c là nhỏ hơn giá trị nhỏ nhất của βm. Đối với những giá trị của c lớn hơn giá trị nhỏ nhất của βm thì sẽ có ít nhất một lớp mà có Qm là đại lượng thực. Trong trường hợp này, các tỉ số K/L và M/N sẽ không tiến đến một giá trị tiệm cận khi k lớn, nhưng sẽ là một hàm dao động đối với k. Do đó, đối với một giá trị đã cho của c mà lớn hơn giá trị nhỏ nhất của βm và nhỏ hơn βn thì nói chung sẽ có vô số các giá trị của k sao cho KN−ML = 0. Những nghiệm này biểu diễn dãy các nghiệm của mode trực giao. Bởi vì chỉ có duy nhất một nghiệm tại tần số rất thấp, nên mỗi một mode trực giao này sẽ có một tần số giới hạn dưới. 20 Chương 3 Tính toán số Việc tính toán số hàm c(k) từ phương trình (1.28) sẽ được lập trình tính toán bởi ngôn ngữ lập trình Matlab. Nói chung thì chương trình tính toán sẽ tiện lợi hơn khi ta lấy độ dày của lớp đầu tiên là chiều dài đơn vị, ρ1 là mật độ khối lượng đơn vị, và β1 là vận tốc đơn vị. Khi đó kết quả tính toán của chúng ta sẽ là mối liên hệ giữa hai đại lượng không thứ nguyên c β1 và kd1. Với lưu ý trên, code của chương trình Matlab này được cho trong phần phụ lục của khóa luận. Để kiểm tra chương trình tính toán số này, chúng ta sẽ tính toán lại một số kết quả đã trình bày trong bài báo của Haskell (1953) [4]. Mô hình ta chọn để tính toán số là mô hình thứ I trong bài báo của Haskell, là một mô hình xấp xỉ của mô hình bề mặt trái đất có hai lớp đặt trên một bán không gian. Mô hình I này được Haskell sử dụng tính toán chỉ minh họa sự tán sắc của sóng Rayleigh cho miền vận tốc nhỏ được đưa ra bởi Guttenberg (1950) [3]. Các tham số cụ thể là Bảng I Lớp α(km/s) β (km/s) ρ(gm/cm3) d(km) 1 6.14 3.39 2.70 13.60 2 5.50 3.18 2.70 11.85 3 8.26 4.65 3.00 ∞ Hình 2 biểu diễn vận tốc sóng Rayleigh của mô hình trên cho khoảng hơn mười mode đầu tiên. Trong bài báo của Haskell, ông đã vẽ mode cơ bản (mode 0), và các kết quả của ông được quan sát bằng mắt rất giống kết quả trong hình vẽ. 21 CHƯƠNG 3. TÍNH TOÁN SỐ k(d1+d2) c/β 1 0 5 10 15 20 25 30 35 40 0.8 0.9 1 1.1 1.2 1.3 mode 0 mode 1 mode 2 Hình 2 : Một số mode đầu tiên của sóng Rayleigh. Trục tung là tỷ số của vận tốc sóng chia cho vận tốc sóng dài của lớp đầu tiên. Trục hoành là kd. Hình 3 vẽ một số mode đầu tiên của vận tốc sóng Rayleigh cho mô hình trên theo đại lượng vô hướng là tỷ số của độ dày của lớp đầu tiên và bước sóng của sóng ngang trong lớp đầu tiên này, tính theo chương trình của Herrmann và chương trình Matlab của khóa luận. Chúng ta có thể thấy rõ ràng rằng, do việc chia miền tính toán của tần số thành các đoạn rời rạc, nên chương trình Herrmann cho kết quả không mịn, mặc dù trong hình vẽ trên ta đã chia khoảng tần số thành 1024 khoảng con. Hình vẽ cũng minh họa được độ tin cậy của chương trình Matlab của chúng ta khi cho kết quả phù hợp với các kết quả tính toán của chương trình Herrmann. 22 CHƯƠNG 3. TÍNH TOÁN SỐ 0 0.5 1 1.5 2 0.8 0.9 1 1.1 1.2 1.3 d.f/β1 c/β 1 Hình 3 : Một số mode của đường cong tán sắc của sóng Rayleigh cho mô hình ở trên. Các đường nét liền được vẽ bằng chương trình Matlab của khóa luận. Các được nét đứt được vẽ bởi chương trình Herrmann. Trục tung là tỷ số của vận tốc sóng chia cho vận tốc sóng dài của lớp đầu tiên. Trục hoành là kd. 23 Kết luận Khóa luận đã tìm hiểu và trình bày một cách có hệ thống các kết quả trong hai bài báo của Thomson (1950) [9] và của Haskell (1953) [4]. Các ý tưởng và các phương trình cơ bản của phương pháp ma trận chuyển đã được trình bày rõ ràng trong khóa luận. Dựa vào các kết quả trong hai bài báo trên, tác giả cũng đã đưa ra công thức tính toán tỷ số H/V đối với môi trường phân lớp. Công thức này khá đơn giản và có thể suy ra một cách trực tiếp từ một số công thức trong hai bài báo. Kết quả quan trọng nhất của khóa luận là code chương trình tính toán vận tốc sóng và tỷ số H/V, là hai đại lượng cơ bản quan trọng nhất của sóng Rayleigh, đối với môi trường phân lớp có số lớp bất kỳ được viết bằng ngôn ngữ lập trình Matlab. Chương trình này chạy nhanh và có độ chính xác tương đương với chương trình của Herrman. Mặc dù nó chỉ có thể tính toán được hai đại lượng cơ bản trên, trong khi chương trình Herrman có thể tính toán được rất nhiều đại lượng khác của sóng Rayleigh như là chuyển dịch, ứng suất, hàm năng lượng, v.v... nhưng ưu điểm của nó là người sử dụng có thể thay đổi trực tiếp để phục vụ các tính toán đặc biệt của mình, ví dụ như tính toán xung quanh tần số điểm osculation point khi độ mịn của tần số yêu cầu là rất lớn. Hoặc là khi cần lập các bản đồ tổng thể của tính chất của sóng Rayleigh. Các bản đồ này khi dùng chương trình của Herrman có độ mịn kém và do đó cho một bản đồ trông không được liên tục, nhất là xung quanh các điểm đặc biệt. Khóa luận cũng trình bày một số các kết quả tính toán số sử dụng chương trình tính toán trên để so sánh với các chương trình tính toán khác. 24 Tài liệu tham khảo [1] Đào Huy Bích (2000), Lý thuyết đàn hồi, Nhà xuất bản Đại Học Quốc Gia Hà Nội, Hà Nội. [2] Crampin S., Taylor D. B. (1971), "The propagation of surface waves in anisotropic media", Geophys. J. R. astr. Soc., 25, pp.71-87. [3] Gutenberg B. (1950), "The structure of the crust in the continents", Sci- ence, 111:50. [4] Haskell N. A. (1953), "The dispersion of surface waves on multilayered media", Bulletin of the Seismological Society of America, 43(1), pp.17-34. [5] Herrmann R. B. (1994), Computer Programs in Seismology, St Louis Uni- versity. [6] Malischewsky P. G. (2004), "A note on Rayleigh-wave velocities as a func- tion of the material parameter", Geofisica Internacional, 45, pp.507-509. [7] Sezawa K., Kanai K. (1938), "Anomalous dispersion of elastic surface waves", Bull. Earthq. Res. Inst. Tokyo, 16: 683. [8] Stoneley R. S. (1924), "Elastic waves at the surface of separation of two solids", Proc. Roy. Soc., 106:416. [9] Thomson W. T. (1950), "Transmission of elastic waves through a stratified solid medium", Jour. Appl. Phys., 21, 89. [10] Tran Thanh Tuan (2009), The ellipticity (H/V-ratio) of Rayleigh surface waves, Dissertation in GeoPhysics, Friedrich-Schiller-University Jena, Germany. 25 TÀI LIỆU THAM KHẢO [11] Tran Thanh Tuan, Frank Scherbaum, Malischewsky P. G. (2011), "On the relationship of peark and troughs of the ellipticity (H/V) of Rayleigh waves anh the transmission response of single layer over half-space models", Geophys. J. Int., 184, pp.793-800. 26 27 PHỤ LỤC Code chương trình Matlab function run a=[6140 5500 8260];% van toc song doc b=[3390 3180 4650];% van toc song ngang rho=[2700 2700 3000];% he so rho d=[13600 11850]./13600; %chieu day cac lop figure(1) ff=@(kd,C) myfun(kd/sum(d),C*b(1),a,b,rho,d); h=ezcontour(ff,[0 20 0.1 1.4],60); set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2,'LevelList',[0 ]) figure(2) ff=@(kd1) velocity(kd1,a,b,rho,d)/b(1); fplot(ff,[0.0001 20]) figure(3) gg=@(kd1) HVratio(kd1,a,b,rho,d); fplot(gg,[0 3]) end function out=HVratio(kd1,a,b,rho,d) c=velocity(kd1,a,b,rho,d); n=length(a); A=eye(4,4); for j=1:(n-1) mtra=zeros(4,4); g=2*(b(j)/c)^2; if c>a(j) ra=sqrt(((c/a(j))^2-1)); else ra=-1i*sqrt((1-(c/a(j))^2)); end if c>b(j) rb=sqrt(((c/b(j))^2-1)); else rb=-1i*sqrt((1-(c/b(j))^2)); end 28 P=kd1*ra*d(j);Q=kd1*rb*d(j); mtra(1,1)=g*cos(P)-(g-1)*cos(Q); mtra(1,2)=1i*((g-1)/ra*sin(P)+g*rb*sin(Q)); mtra(1,3)=-1/(rho(j)*c^2)*(cos(P)-cos(Q)); mtra(1,4)=1i/(rho(j)*c^2)*(sin(P)/ra+sin(Q)*rb); mtra(2,1)=-1i*(g*ra*sin(P)+(g-1)/rb*sin(Q)); mtra(2,2)=-(g-1)*cos(P)+g*cos(Q); mtra(2,3)=1i/(rho(j)*c^2)*(ra*sin(P)+sin(Q)/rb); mtra(2,4)=mtra(1,3); mtra(3,1)=rho(j)*c^2*g*(g-1)*(cos(P)-cos(Q)); mtra(3,2)=1i*rho(j)*c^2*((g-1)^2/ra*sin(P)+g^2*rb*sin(Q)); mtra(3,3)=mtra(2,2); mtra(3,4)=mtra(1,2); mtra(4,1)=1i*rho(j)*c^2*(g^2*ra*sin(P)+(g-1)^2/rb*sin(Q)); mtra(4,2)=mtra(3,1); mtra(4,3)=mtra(2,1); mtra(4,4)=mtra(1,1); A=mtra*A; end g=2*(b(n)/c)^2; if c>a(n) ra=sqrt(((c/a(n))^2-1)); else ra=-1i*sqrt((1-(c/a(n))^2)); end if c>b(n) rb=sqrt(((c/b(n))^2-1)); else rb=-1i*sqrt((1-(c/b(n))^2)); end %------------------------------------------------------------------- K=g*ra*A(1,2)+(g-1)*A(2,2)- ra*A(3,2)/(rho(n)*c^2)+A(4,2)/(rho(n)*c^2); L=g*ra*A(1,1)+(g-1)*A(2,1)- ra*A(3,1)/(rho(n)*c^2)+A(4,1)/(rho(n)*c^2); ratio=-K/(1i*L); out=ratio; end function out=velocity(kd1,a,b,rho,d) ff=@(c) myfun(kd1,c,a,b,rho,d); n=50; kq=[]; Cmax=max(b);Cmin=0.6*min(b); dx=[Cmin:(Cmax-Cmin)/n:Cmax]; for i=1:length(dx)-1 if ff(dx(i))*ff(dx(i+1))<0 [x,fval,exitflag]=fzero(ff,[dx(i) dx(i+1)]); 29 if exitflag==1 out=x; return end end end while isempty(kq)==1 n=n+50; if n>300 out=NaN; return end dx=[0.8:(Cmax-0.8)/n:Cmax]; for i=1:length(dx)-1 if ff(dx(i))*ff(dx(i+1))<0 [x,fval,exitflag]=fzero(ff,[dx(i) dx(i+1)]); if exitflag==1 out=x; return end end end end end function out=myfun(kd1,c,a,b,rho,d) % d=[1 d2/d1 d3/d1 ......] chieu day cua cac lop n=length(a); A=eye(4,4); for j=1:(n-1) mtra=zeros(4,4); g=2*(b(j)/c)^2; if c>a(j) ra=sqrt(((c/a(j))^2-1)); else ra=-1i*sqrt((1-(c/a(j))^2)); end if c>b(j) rb=sqrt(((c/b(j))^2-1)); else rb=-1i*sqrt((1-(c/b(j))^2)); end P=kd1*ra*d(j);Q=kd1*rb*d(j); mtra(1,1)=g*cos(P)-(g-1)*cos(Q); mtra(1,2)=1i*((g-1)/ra*sin(P)+g*rb*sin(Q)); mtra(1,3)=-1/(rho(j)*c^2)*(cos(P)-cos(Q)); mtra(1,4)=1i/(rho(j)*c^2)*(sin(P)/ra+sin(Q)*rb); mtra(2,1)=-1i*(g*ra*sin(P)+(g-1)/rb*sin(Q)); mtra(2,2)=-(g-1)*cos(P)+g*cos(Q); mtra(2,3)=1i/(rho(j)*c^2)*(ra*sin(P)+sin(Q)/rb); mtra(2,4)=mtra(1,3); mtra(3,1)=rho(j)*c^2*g*(g-1)*(cos(P)-cos(Q)); mtra(3,2)=1i*rho(j)*c^2*((g-1)^2/ra*sin(P)+g^2*rb*sin(Q)); 30 mtra(3,3)=mtra(2,2); mtra(3,4)=mtra(1,2); mtra(4,1)=1i*rho(j)*c^2*(g^2*ra*sin(P)+(g-1)^2/rb*sin(Q)); mtra(4,2)=mtra(3,1); mtra(4,3)=mtra(2,1); mtra(4,4)=mtra(1,1); A=mtra*A; end g=2*(b(n)/c)^2; if c>a(n) ra=sqrt(((c/a(n))^2-1)); else ra=-1i*sqrt((1-(c/a(n))^2)); end if c>b(n) rb=sqrt(((c/b(n))^2-1)); else rb=-1i*sqrt((1-(c/b(n))^2)); end %------------------------------------------------------------------- K=g*ra*A(1,2)+(g-1)*A(2,2)- ra*A(3,2)/(rho(n)*c^2)+A(4,2)/(rho(n)*c^2); L=g*ra*A(1,1)+(g-1)*A(2,1)- ra*A(3,1)/(rho(n)*c^2)+A(4,1)/(rho(n)*c^2); M=-(g- 1)*A(1,2)+g*rb*A(2,2)+A(3,2)/(rho(n)*c^2)+rb*A(4,2)/(rho(n)*c^2); N=-(g- 1)*A(1,1)+g*rb*A(2,1)+A(3,1)/(rho(n)*c^2)+rb*A(4,1)/(rho(n)*c^2); bt=K*N-L*M; tg=double(bt); out=tg; end