雷克子波 因计算机资源有限,不便在太大的区域做PML 的计算,一般只在边界上取有限宽度的区域作为PML 计算区域。根据相关文献的研究,PML 边界区域最少长度应为半个波长[6]。本文综合考虑了效果与开销等因素,选取了边界上 50 层作为 PML 的计算区域。常规计算区域与 PML 边界区域的如图 2-4 所示。 xxxxxPvaPKtx zzzzzPvaPKtz xzxxvvPaPKtxz 衰减系数 2331lo g2xzPvRL 式中 L 为 PML 层的厚度, 为层内的点距 PML 与非 PML 的边界的距离,Pv 为纵波速度, R 那么在 PML 边界区域内,对于式(2-13)1xxxvPvtx 即为理论反射系数,一般取 0.001 较为合适, 为 方向的空间步长。 1xxxvPvtx,可看作为在常规的计算方程基础上,减去一项进行 PML 的阻尼修正项。因本文中只考虑各项同性介质中的地震波传播规律,故可做xz假设。在此利用一下三个假设: 2121,2121,21,21kjixkjixkjixvvv 21,2121,21,2121kjizkjizkjizvvv 1,1,21,21kjikjikjiPPP 因为以上三个近似精度均为时间方向上的近似,且时间精度均为二阶精度,因交错网格技术的时间精度为二阶,故以上近似不影响本式的计算精度。 故可得: NmkmjikmjimkjikjixkjixPPaxtvttv11,,21,2121,2121,15.015.011 (2-18a) 同理: NmkjmikjmimkjikjizkjizPPaztvttv1,1,21,21,2121,2115.015.011 (2-18b) NmkmjizkmjixmNmkjmizkjmizmkjikjikjivvaxtvvaztKPttP12121,2121,121,2121,2121,,1,5.015.011 (2-18c) 此为在PML边界区域内的弹性声波应力-速度方程组。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%交错网格---非均匀介质二维声波方程(一阶压力--速度)、2阶时间差分、2阶空间差分精度 %%加上边界 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close all;clea...