CN110348158A - 一种基于分区异步长解的地震波动分析方法 - Google Patents

Abstract

本发明公开了一种基于分区异步长解的地震波动分析方法，包括S1、按照网格尺寸将计算区域分为I区域和II区域；S2、II区域按照计算时步取Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n时，对I、II区域进行求解，确定波动传播的初始位移、初始速度和初始加速度；S3、根据初始时刻t的加速度值，计算得到t+Δt时刻波动传播的位移与速度的预测值；S4、根据所述波动传播的位移和速度预测值，采用动力平衡方程求解t+Δt时刻的加速度，并确定该时刻波动传播的位移和速度；S5、重复步骤S2‑S4，直至全部时间结束。

Description

一种基于分区异步长解的地震波动分析方法

技术领域

本发明属于人工地震模拟的技术领域，具体涉及一种基于分区异步长解的地震波动分析方法。

背景技术

地震波动分析，目前常采用时域显式有限元方法^[1～3]，主要特点是不必形成总体刚度矩阵，利用中心差分和质量矩阵的对角化，使用统一时间步长，在“单元级”上求解，这样可以使用较少的内存，但需要频繁的内外存交换。由于实际问题的复杂性，有限元网格不可能均匀一致，而在整个求解区域使用的统一时间步长由最小单元尺寸确定，因此，这种解法计算效率较低。随着微机性能的提高，特别是微机内存配置的增大，这种不需形成总的刚度矩阵而频繁使用内外存交换的时域显式有限元方法，其缺点是显而易见的。

提高计算效率的有效途径之一是充分利用微机内存^[4]。另一提高计算速度的方法是分区算法，朱伯芳院士在求解不稳定温度场时采用了分区异步长算法^[5]，即在温度急剧变化的区域，采用小时间步长，在其余区域采用大时间步长。文献6、7在波动问题的计算中提出了一种显隐式求解方法，即在局部区域进行隐式积分，而在其它区域进行显式积分，这样在保证计算精度的前提下，可以使用较大的时间步长。

[参考文献]

[1]涂劲，有缝界面的混凝土高坝—地基系统非线性地震波动反应分析[D]，中国水利水电科学研究院博士学位论文，1999年7月.

[2]张晓志，谢礼立，采用显式有限元方法求解开放系统动力响应的主要环节和问题[J]，地震工程与工程振动，2005，25(2):10～15.

[3]李亮，杜修力，赵成刚，李立云，基于显式有限元方法的两相多孔介质地震反应研究[J]，应用力学学报，2007，24(4):550～553.

[4]张伯艳，高拱坝坝肩抗震稳定研究[D]，西安理工大学博士学位论文，2005年.

[5]朱伯芳，有限单元法原理与应用(第二版)[M]，北京：中国水利水电出版社，1998年.

[6]Hughes T J R,Liu W K.Implicit-explicit finite elements intransient analysis:implementation and numerical examples[J].Journal ofApplied Mechanics,1978,45:375～378.

发明内容

本发明的目的在于针对现有技术中的上述不足，提供一种基于分区异步长解的地震波动分析方法，以解决或改善上述的问题。

为达到上述目的，本发明采取的技术方案是：

一种基于分区异步长解的地震波动分析方法，其包括：

S1、按照网格尺寸将计算区域分为I区域和II区域；

S2、II区域按照计算时步取Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n时，对I、II区域进行求解，确定波动传播的初始位移、初始速度和初始加速度；其中，Δt为时间步长，Δt_s为网格尺寸较小区域的计算时步，Δt_n为网格尺度较大区域的计算步长，Δt_n＝nΔt_s；

S3、根据初始时刻t的加速度值，计算得到t+Δt时刻波动传播的位移与速度的预测值；

S4、根据波动传播的位移和速度预测值，采用动力平衡方程求解t+Δt时刻的加速度，并确定该时刻波动传播的位移和速度；

S5、重复步骤S2-S4，直至全部时间结束。

优选地，步骤S1的具体步骤包括：

按照网格尺寸对计算区域进行分区，生成所需要的数据文件，得到I、II两个区域，其中，I区域的网格尺度大于II区域网格尺寸；且I、II两个区域的计算时步由各自区域内的最小网格尺寸决定。

优选地，步骤S2的具体步骤包括：

采用PCA算法，计算波动传播的位移和速度的预测值：

其中，和为位移和速度的预测值，β和γ是Newmark法的控制参数，且β＝0、γ＝1/2，M为系统的质量矩阵，C为系统的阻尼矩阵，K为系统的刚度矩阵，F_t+Δt为t+Δt时刻的外力荷载向量；

当采用基于显式PCA分区异步长算法进行结构数值求解时，I区域的计算时步Δt_n是II区计算时间步长Δt_s的n倍，则：

其中，M^I,M^II,C^I,C^II,K^I,K^II分别为I、II两个区域的总体质量、阻尼和刚度矩阵，而F^I,F^II分别是I、II两个区域的外荷载向量；

II区按照计算时步为Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n对全部计算区域进行求解，在确定初始位移u₀和初始速度后，根据平衡方程得到波动传播的初始加速度；

其中，M^-1为系统质量矩阵的逆矩阵，F₀为初始时刻的外力荷载向量。

优选地，步骤S3的具体步骤包括：

对II区域建立求解t+Δt_s时刻的位移与速度的预测值：

当计算时步达到Δt_n＝nΔt_s，此时对全部计算区域进行求解，t+Δt_n时刻的位移预测值为速度的预测值为

优选地，步骤S4的具体步骤包括：

对II区求解，基于步骤S3中求解的位移和速度预测值，根据动力平衡方程求解t+Δt_s时刻的加速度为：

确定该时刻的位移和速度为：

当计算时步达到Δt_n＝nΔt_s，此时对全部计算区域进行求解，t+Δt_n时刻的加速度为：

确定该时刻的位移和速度为：

本发明提供的基于分区异步长解的地震波动分析方法，具有以下有益效果：

本发明利用总体刚度矩阵的稀疏性，采用非零存储方法，将总体刚度矩阵存储于一维数组，进一步提高了利用微机进行大型工程问题的动力计算效率，最大限度地使用了计算机内存，减少了内外存交换。同时，从波动方程出发，将非零存储技术与分区显式求解方法相结合，进行地震波动分析，该分析方法具有较高的计算效率，可以为大型工程地震反应分析提供了一个有效的方法和可供选择的工具。

附图说明

图1为基于分区异步长解的地震波动分析方法的流程图。

图2为包含两种材料一维波动传播问题的网格分区图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

根据本申请的一个实施例，参考图1，本方案的基于分区异步长解的地震波动分析方法，包括；

S1、按照网格尺寸将计算区域分为I区域和II区域；

S2、II区域按照计算时步取Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n时，其中，Δt_n＝nΔt_s，对I、II区域进行求解，确定波动传播的初始位移、初始速度和初始加速度；

S3、根据初始时刻t的加速度值，计算得到t+Δt时刻波动传播的位移与速度的预测值；

S4、根据所述波动传播的位移和速度预测值，采用动力平衡方程求解t+Δt时刻的加速度，并确定该时刻波动传播的位移和速度；

S5、重复步骤S2-S4，直至全部时间结束。

本发明从波动方程出发，将非零存储技术与分区显式求解方法相结合，进行地震波动分析方法，该方法具有较高的计算效率，可以为大型工程地震反应分析提供了一个有效的方法和可供选择的工具。

以下对上述各个步骤进行详细说明。

步骤S1，按照网格尺寸对计算区域进行分区，生成所需要的数据文件，包括I、II两个区域，其中，I区的网格尺度较大，II区网格尺寸较小。

参考图2，计算区域分为I、II两个区域，其中I区的网格尺度较大，II区网格尺寸较小，为了保证计算的精度以及稳定性要求，两个分区的计算时步由各自区域内的最小网格尺寸决定。

步骤S2，II区按照计算时步取为Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n(Δt_n＝nΔt_s)对全部计算区域进行求解，确定初始位移和初始速度，并根据平衡方程就算初始加速度，其步骤包括：

为实现波动问题的分区解算，采用了PCA(Predictor-corrector algorithms)算法：

其中，和为位移和速度的预测值，β和γ是Newmark法的控制参数，采用与显式有限元相结合的PCA算法，其中β＝0、γ＝1/2，M为系统的质量矩阵，C为系统的阻尼矩阵，K为系统的刚度矩阵，F_t+Δt为t+Δt时刻的外力荷载向量。

当采用基于显式PCA分区异步长算法进行结构数值求解时，I区的计算时步Δt_n是II区计算时间步长Δt_s的n倍，此时方程(1a)可以写为：

其中，M^I,M^II,C^I,C^II,K^I,K^II分别为I、II两个区域的总体质量、阻尼和刚度矩阵，而F^I,F^II分别是I、II两个区域的外荷载向量。

II区按照计算时步为Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n(Δt_n＝nΔt_s)对全部计算区域进行求解。在确定初始位移u₀和初始速度后，根据平衡方程计算得到初始加速度

步骤S3，根据初始时刻t的加速度值计算求得t+Δt时刻的位移与速度的预测值，包括：

对II区建立求解(小时步Δt_s)，t+Δt_s时刻的位移与速度的预测值为：

当计算时步达到Δt_n＝nΔt_s，此时对全部计算区域进行求解，t+Δt_n时刻的位移预测值为速度的预测值为

步骤S4，基于上述位移和速度预测值，根据动力平衡方程求解t+Δt时刻的加速度，进而确定该时刻的位移和速度，其具体步骤包括：

对II区求解，基于步骤S3中求解的位移和速度预测值，根据动力平衡方程求解t+Δt_s时刻的加速度为：

进而确定该时刻的位移和速度为：

当计算时步达到Δt_n＝nΔt_s，此时对全部计算区域进行求解，t+Δt_n时刻的加速度为：

进而确定该时刻的位移和速度为：

S5：重复步骤S2-S4，直至全部时间结束。

虽然结合附图对发明的具体实施方式进行了详细地描述，但不应理解为对本专利的保护范围的限定。在权利要求书所描述的范围内，本领域技术人员不经创造性劳动即可做出的各种修改和变形仍属本专利的保护范围。

Claims (5)

1.一种基于分区异步长解的地震波动分析方法，其特征在于：包括

S1、按照网格尺寸将计算区域分为I区域和II区域；

S2、II区域按照计算时步取Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n时，对I、II区域进行求解，确定波动传播的初始位移、初始速度和初始加速度；其中，Δt为时间步长，Δt_s为网格尺寸较小区域的计算时步，Δt_n为网格尺度较大区域的计算步长，Δt_n＝nΔt_s；

S3、根据初始时刻t的加速度值，计算得到t+Δt时刻波动传播的位移与速度的预测值；

S4、根据所述波动传播的位移和速度预测值，采用动力平衡方程求解t+Δt时刻的加速度，并确定该时刻波动传播的位移和速度；

S5、重复步骤S2-S4，直至全部时间结束。

2.根据权利要求1所述的基于分区异步长解的地震波动分析方法，其特征在于，所述步骤S1的具体步骤包括：

按照网格尺寸对计算区域进行分区，生成所需要的数据文件，得到I、II两个区域，其中，I区域的网格尺度大于II区域网格尺寸；且I、II两个区域的计算时步由各自区域内的最小网格尺寸决定。

3.根据权利要求1所述的基于分区异步长解的地震波动分析方法，其特征在于，所述步骤S2的具体步骤包括：

采用PCA算法，计算波动传播的位移和速度的预测值：

其中，和为位移和速度的预测值，β和γ是Newmark法的控制参数，且β＝0、γ＝1/2，M为系统的质量矩阵，C为系统的阻尼矩阵，K为系统的刚度矩阵，F_t+Δt为t+Δt时刻的外力荷载向量；

当采用基于显式PCA分区异步长算法进行结构数值求解时，I区域的计算时步Δt_n是II区计算时间步长Δt_s的n倍，则：

其中，M^I,M^II,C^I,C^II,K^I,K^II分别为I、II两个区域的总体质量、阻尼和刚度矩阵，而F^I,F^II分别是I、II两个区域的外荷载向量；

II区按照计算时步为Δt＝Δt_s进行求解，当计算时步达到Δt＝Δt_n对全部计算区域进行求解，在确定初始位移u₀和初始速度后，根据平衡方程得到波动传播的初始加速度；

其中，M^-1为系统质量矩阵的逆矩阵，F₀为初始时刻的外力荷载向量。

4.根据权利要求1所述的基于分区异步长解的地震波动分析方法，其特征在于，所述步骤S3的具体步骤包括：

对II区域建立求解t+Δt_s时刻的位移与速度的预测值：

当计算时步达到Δt_n＝nΔt_s，此时对全部计算区域进行求解，t+Δt_n时刻的位移预测值为速度的预测值为

5.根据权利要求1所述的基于分区异步长解的地震波动分析方法，其特征在于，所述步骤S4的具体步骤包括：

对II区求解，基于步骤S3中求解的位移和速度预测值，根据动力平衡方程求解t+Δt_s时刻的加速度为：

确定该时刻的位移和速度为：

当计算时步达到Δt_n＝nΔt_s，此时对全部计算区域进行求解，t+Δt_n时刻的加速度为：

确定该时刻的位移和速度为：