计算机执行确定特大型工程地震安全性预测的方法

摘要

一种计算机技术领域的计算机执行确定特大型工程地震安全性预测的方法，用有限元前处理软件建立有限元模型，并输入地震动时程所需参数；参数传递到超级计算机，用人造地震动时程并行模拟的混合任务分配方法得到地震激励动时程；地震激励动时程返回到有限元前处理软件进行约束加载，得到有限元分析软件的任务文件，并提交给有限元分析软件，得到数学模型，数学模型和时间步长传递至超级计算机，用瞬态响应分析和模态分析的并行处理方法或模态分析的并行处理方法并行计算；将数据返回有限元分析软件进行后续计算；有限元后处理软件读入计算结果，进行可视化处理。本发明采用软件方式实现，大大节省系统实施的人力、物力和财力。

说明书

技术领域

本发明涉及的是一种计算机技术领域的方法，具体是一种计算机执行确定特大型工程地震安全性预测的方法。

背景技术

地震安全性预测是指在对具体建设工程场址及其周围地区的地震地质条件、地球物理场环境、地震活动规律、现代地形变及应力场等方面深入研究的基础上，采用先进的地震危险性概率分析方法，按照工作所需要采用的风险水平，科学地给出相应的工程规划或设计所需要的一定概率水平下的地震动参数(加速度、设计反应谱、地震动时程等)和相应的资料。地震安全性预测分为两步：第一步，计算地面的地震动参数；第二步，计算工程结构的地震动参数。对特大型工程和结构进行地震安全性预测，可以采用三维有限元方法。由于规模巨大、结构复杂，特大型工程和结构的瞬态响应和模态计算以及可视化显示都需要海量计算、海量数据交换和存储来支撑，如果采用普通的计算机及串行的有限元软件，则有限元网格划分较粗，这样所得到的结果往往缺乏准确性，精度也就受到限制。另外，多点地震动模拟是一个复杂的计算过程，要使人工模拟的多点地震动输入更为精确，空间相关点的距离应尽量小。由于每个空间点都有几千个频率采样点，随着空间采样点的个数增多，计算量也随之增大很多。因此，必须采用超级计算机和并行有限元软件。

对于并行软件来说，不同的应用领域、不同的计算机硬件条件，并行软件的性能差异很大。由于软件技术封闭性强，其工程应用受到极大的限制。

经对现有技术的文献检索发现，沈建文等在《地震学刊》1997年第4期中发表了“XQH1.00地震安全性预测工作程序包简介”一文，该程序包是串行的，对大型工程结构的求解效率低下，地震安全性预测的精度也就不高。在进一步检索中还发现，丁海平等在《地震工程与工程振动》2004年第2期中发表了“三维地震动场数值模拟并行计算系统”一文，该系统仅适用于场地地震动场数值模拟，未涉及工程结构的地震动参数计算，而且是在若干台普通微机构成的“Beowulf并行系统”上实现的。

发明内容

本发明针对现有技术中存在的上述问题，提供一种计算机执行确定特大型工程地震安全性预测的方法，使其根据特大型工程地震安全性预测的特点和要求，借助超级计算机的高速计算和海量数据处理能力，提高地震安全性预测的准确性和精度。能为基于超级计算机的特的大型工程地震安全性预测并行软件的开发提供指导。

本发明是通过以下的技术方案实现的，本发明步骤如下：

(1)首先采用有限元前处理软件建立有限元模型，通过相应界面输入地震动时程所需要的参数：设计反应谱、近震或远震、抗震设防烈度、土壤类型、发生频率、阻尼常数、地震动持时、地震动步长，对于多点地震动模拟，所需参数还包括：自功率谱和互功率谱模型、空间点数目、相邻点距离、视速度；

(2)将步骤(1)得到的参数传递到超级计算机，根据工程实际情况，采用地震危险性概率分析方法，按照国家标准要求的风险水平，通过人造地震动时程并行模拟的混合任务分配方法，并行计算得到地震激励动时程；

(3)将步骤(2)中得到的地震激励动时程返回到有限元前处理软件进行约束加载，得到有限元分析软件的任务文件；

(4)将任务文件提交给有限元分析软件，由有限元分析软件进行有限元分析计算，得到相应的数学模型：刚度矩阵、阻尼矩阵、质量矩阵和载荷；接着将得到的数学模型和时间步长传递至超级计算机，采用瞬态响应分析和模态分析的并行处理方法进行地震响应的并行计算，得到位移、速度、加速度，或采用模态分析的并行处理方法并行计算特大型工程的固有频率、固有振型等动态特性；

(5)将步骤(4)中得到的数据返回有限元分析软件进行后续的计算，得到每一时刻的应力和应变；

(6)采用有限元后处理软件读入有限元分析软件计算的结果，进行可视化处理，通过可视化数值仿真，得到特大工程地震安全性预测的条件，为特大型工程的抗震设计提供指导。特大工程地震安全性预测的条件包括：工程结构任何位置在任意时刻的变形和受力状态，关键部位相对位移随时间和空间的变化规律，地震激励方向对相对位移的影响规律。以及对工程结构进行综合性评价条件，包括：地震对结构可能造成损坏的形式、地震对结构可能造成损坏的条件、地震对结构可能造成损坏的位置、地震对结构可能造成损坏的工况、结构抗震性能的评价。

所述的地震激励动时程，计算步骤如下：

(1)确定所采用的功率谱模型

选用Chough-Penzien谱，也称为双过滤白噪声谱，其模型公式为：

$>>S>>(>\omega >)>>=>>>>\omega >g>4>>+>4>>>\zeta >g>4>>>\omega >g>2>>>\omega >2>>>>>>(>>\omega >2>>->>\omega >g>2>>)>>2>>+>4>>\zeta >g>2>>>\omega >g>2>>>\omega >2>>>>·>>>\omega >4>>>>>(>>\omega >2>>->>\omega >f>2>>)>>2>>+>4>>\zeta >f>2>>>\omega >f>2>>>\omega >2>>>>>S>0>>->->->>(>1>)>>>s>$

式中ω_g为场地卓越频率。

(2)确定互功率谱

两点间的互功率谱与两点间的自功率谱、相干函数及波的视速度有关。小范围内地震动空间相干函数的方位角度变化是可以忽略的。基于现有的研究结果，从使用的角度考虑，可暂且认为小范围地震动相干函数不随方位变化。

$>>>S>cr>>>(>\omega >,>d>)>>=>>>S>i>>>(>\omega >)>>>S>j>>>(>\omega >)>>>\rho >>(>\omega >,>d>)>>exp>[>->i\omega >>d>>C>>(>\omega >)>>>>]>->->->>(>2>)>>>s>$

其中S_i(ω)和S_j(ω)分别为i，j两点的自功率谱密度函数；ρ(w，d)为两点的相干函数；C(ω)＝C_H(ω)/cosα，C_H(ω)为地震波的水平视波速，α为i→j方向与地震波传播方向的夹角。

(3)功率谱矩阵分解

功率谱矩阵为

$>>S>>(>>i\omega >k>>)>>=>\begin{array}{}>>>>S>1>>>(>>\omega >k>>)>>>>>S>12>>>(>i>>\omega >k>>)>>>>·>·>·>>>>S>>1>n>>>>(>i>>\omega >k>>)>>>>>>>S>21>>>(>>i\omega >k>>)>>>>>S>2>>>(>>\omega >k>>)>>>>·>·>·>>>>S>>2>n>>>>(>i>>\omega >k>>)>>>>>>·>>>·>>>·>·>·>>>·>>>>>·>>>·>>>·>·>·>>>·>>>>>·>>>·>>>·>·>·>>>·>>>>>>S>>n>1>>>>(>>\omega >k>>)>>>>>S>>n>2>>>>(>>\omega >k>>)>>>>·>·>·>>>>S>n>>>(>>\omega >k>>)>>>\end{array}>>->->->>(>3>)>>>s>$

进行Cholesky分解得到下三角矩阵

$>>L>>(>i>>\omega >k>>)>>=>\begin{array}{}>>>>l>11>>>(>>\omega >k>>)>>>>0>>>·>·>·>>>0>>>>>>l>21>>>(>i>>\omega >k>>)>>>>>l>22>>>(>>\omega >k>>)>>>>·>·>·>>>0>>>>>·>>>·>>>>>·>>>>>·>>>·>>>·>·>·>>>·>>>>>·>>>·>>>>>·>>>>>>l>>n>1>>>>(>>\omega >k>>)>>>>>l>>n>2>>>>(>>\omega >k>>)>>>>·>·>·>>>>l>nn>>>(>>\omega >k>>)>>>\end{array}>>->->->>(>4>)>>>s>$

求解相关幅值与相关相位谱

$>>>a>nm>>=>>4>\Delta \omega >>|>>l>nm>>>(>i>>\omega >k>>)>>|>->->->>(>5>)>>>s>$

$>>>\theta >nm>>=>>tg>>->1>>>>>>I>m>>[>>l>nm>>>(>i>>\omega >k>>)>>]>>>Re>[>>l>nm>>>(>i>>\omega >k>>)>>]>>>->->->>(>6>)>>>s>$

(4)确定相位差谱

脉动相位差谱

ε(f)＝Δ(f)-Δ(f)                                        (7)

其中Δ(f)为相位差谱的均值函数。Δ(f)与震级、震源距的关系为

log₁₀|Δ(f)|＝a₁(f)+a₂(f)·M+a₃(f)log₁₀(R+R₀)               (8)

Δ(f)＝-|Δ(f)|＜0                                        (9)

其中M为震级，R为震源距，R₀为常数(R₀＝15.0)，a1(f)，a2(f)，a3(f)为与频率有关的回归系数。

(5)由多点地震动合成公式求多点地震动时程

根据空间相关多点地震动合成公式，每个点的合成公式都表示成n项三角级数的和，即

……

上式可写成

……

由上式可生成各点的地震动傅立叶谱，然后用逆傅立叶变换技术得到各点的平稳地震动时程。

所述的地震激励动时程，采用并行处理混合任务分配方法，空间采样点与频率采样点相结合，具体地，相关功率谱求解按频率点分配的任务分配方法，相位差谱递推求解按空间点分配的任务分配方法。

所述的瞬态响应分析和模态分析的并行处理方法，具体采用以下的处理方法：(1)模块并行：当形成质量、刚度和阻尼矩阵以后，可以同时进行瞬态响应计算和模态分析计算。(2)工况并行、区间并行：在瞬态响应计算时，可以同时进行多工况的瞬态响应计算；在模态分析时，可以同时在多个频率区间[F(i-1)，F(i)]进行模态分析。

其中，有限元前、后处理软件、有限元分析软件可借助现有技术中的成熟CAD/CAE软件来完成。采用集成软件来实现上述各软件间的无缝集成，包括可视化操作、数据传递、软件远程启动、数据返回等。

有限元前处理软件的主要功能是：数据文件的输入、输出；实体建模与参数化建模；生成几何图形和图形显示及输出；构件的布尔运算；自动网格划分；载荷与材料参数值的直接输入；有限元模型信息自动生成；形成并提交有限元分析软件能识别的有限元模型数据文件等。因此，前处理软件要处理多种文件：读写描述有限元模型的数据库文件；生成分析软件可识别的输入数据文件；读写并记录各种操作的命令、编辑过程文件等。

有限元后处理软件的主要功能是：处理分析软件得到后的结果，包括识别分析软件生成的格式化及非格式化结果文件；将结果写成各种常用的图形文件以及以图形、动画、曲线、表格和文件形式显示分析软件分析后的结果等。

本发明利用有限元方法计算特大型工程地震响应，其中“求解位移方程组的地震响应”最耗时，对此类动力学问题，计算时间往往占整个过程总时间的90％以上。因此，通过并行化“求解位移方程组的地震响应”能有效地缩短整个分析过程的时间。另外，人造地震动时程模拟的计算量也很大，借助超级计算机和并行计算程序完成。

本发明集成了并行计算软件和现有的有限元前/后处理、分析软件。有限元前、后处理和一部分计算由现有有限元软件完成，大量的计算(系统方程组求解)由超级计算机完成。这样一种并行分析系统，不但拥有并行计算软件的高速计算功能，还拥有现有有限元软件强大的前、后处理功能，及其友好的客户化图形界面，方便用户的使用。本发明提供的方法，能大大节省系统实施的人力、物力和财力。

附图说明

图1为人造地震动时程并行模拟的混合任务分配方法示意图；

图2为地震安全性预测中瞬态响应分析的并行处理方法示意图；

图3为地震安全性预测中模态分析的并行处理方法示意图；

图4为本发明实施例实现原理图。

具体实施方式

如图1所示，人造地震动时程并行模拟的混合任务分配方法，其中i表示按空间采样点分配并行，k表示按频率采样点分配并行。对于空间相关多点动时程的人工模拟，整个任务由多个空间采样点组成，每个空间采样点对应一个时程，每个时程由若干频率采样点组成。在人造地震动时程并行计算软件中，可以按空间采样点分配任务，每个CPU计算几条完整的时程曲线；也可以按频率点分配任务，每个CPU计算每条时程曲线的一段频率范围的动时程。由于相干性和谱分析的需要，按空间采样点和按频率采样点进行任务分配时，CPU之间都需要不断进行通信以交互信息。为此，采用空间采样点与频率采样点交叉结合的任务分配方法进行并行计算，具体地，相关功率谱求解按频率点分配的任务分配方法，相位差谱递推求解按空间点分配的任务分配方法。

如图2和图3所示，地震安全性预测中瞬态响应分析和模态分析的并行处理方法。在地震响应并行计算软件中，运行效率是一关键问题。可采用以下处理方法：(1)模块并行：当形成质量、刚度和阻尼矩阵以后，可以同时进行瞬态响应计算和模态分析计算。(2)工况并行、区间并行：在瞬态响应计算时，可以同时进行多工况的瞬态响应计算；在模态分析时，可以同时在多个频率区间[F(i-1)，F(i)]进行模态分析。

结合本发明的内容提供以下具体实施例：

本实施例基于Unix平台和MPI采用Fortran语言开发了人造地震动时程并行计算软件、结构动态特性和地震响应并行计算软件。MPI(Message PassingInterface)是目前国际上最流行的并行编程环境之一，尤其是分布式存储的可缩放并行计算机和工作站网络的一种编程范例。人造地震动时程并行计算软件的主要功能是：根据工程场址及其周围地区的地震地质条件、地球物理场环境、地震活动规律、现代地形变及应力场等，采用先进的地震危险性概率分析方法，按照相应国家标准要求的风险水平，并行计算出地震激励动时程。结构动态特性和地震响应并行计算软件的主要功能是：并行计算特大型工程的固有频率、固有振型等动态特性；并行计算特大型工程上任意位置在任意时刻的位移、速度、加速度和应力等动态响应参数。

本实施例还采用了集成软件，其主要功能是：实现上述各软件间的无缝集成，包括可视化操作、数据传递、软件远程启动、数据返回等。

基于PATRAN、NASTRAN实现本实施例，其中NASTRAN、PATRAN分别是目前主流的有限元分析软件和前/后处理软件，由美国MSC公司开发。集成软件系统在“神威I”超级计算机(由江南计算技术研究所研制)上实现，该机具有6台前端机，与主机之间通过高速网络连接，所有对主机的操作均是通过前端机实现的。NASTRAN安装在前端机上，PATRAN安装在PC上，并行计算程序运行在“神威I”中。利用本实施例，对一大型隧道进行地震安全性预测，该隧道由若干管段拼接而成。

如图4所示，人造地震动时程并行计算软件、结构动态特性和地震响应并行计算软件运行于超级计算机上，有限元分析软件NASTRAN运行于“神威I”计算系统的前端机上，其他软件均运行于微机上。具体流程如下：

(1)利用PATRAN提供的PCL语言对PATRAN进行二次开发：把集成软件一集成到PATRAN软件系统中，并生成新的客户化图形界面。用户通过PATRAN建立大型隧道的有限元模型(采用六面体单元来划分网格，共划分单元数为706778，节点数为909878)，之后由图形界面调用集成软件一将地震动时程所需参数(设计反应谱选用China(GB50011-2001)，自功率谱选用Kanai-Tajimi模型，互功率谱选用屈铁军模型，IV类场地，震级7.0，距离震源为50m，场地卓越频率为15.71rad/s，视速度为500m/s，截止频率为25Hz，频率离散点为2048个，持续时间为40s，时间步长为0.024s，相邻空间点间距离为60m，阻尼为0.05)传递至超级计算机并远程启动人造地震动时程并行计算软件(根据图1，基于MPI开发)，完成后将人造地震动时程返回到PATRAN进行约束加载；

(2)利用PATRAN提供的PCL语言对PATRAN进行二次开发：①生成NASTRAN任务文件的同时对其进行修改(根据NASTRAN二次开发工具DMAP的规则)，在不影响求解序列其它部分(单元刚度矩阵求解、整体刚度矩阵装配、求解应力应变等)的前提下，调用并行计算软件实现位移方程组的求解(集成软件三，利用DMAP的ISHELL模块实现)，并代替NASTRAN的位移方程组求解模块，输出并行求解所需要的刚度、质量和阻尼矩阵，以及载荷向量；②把集成软件二集成到PATRAN软件系统中，并生成新的客户化图形界面。用户通过图形界面实现对NASTRAN任务文件的修改，并调用集成软件二将任务文件提交并远程启动NASTRAN；

(3)在NASTRAN进行有限元计算的过程中，当进行到“求解位移方程组的地震响应”时跳出NASTRAN，集成软件三将地震响应所需参数(时间步长0.024s、数据文件名)和数据文件(刚度、质量和阻尼矩阵文件，载荷向量文件)传递至超级计算机并远程启动结构动态特性或地震响应并行计算软件(根据图2～3，基于MPI开发)，完成后将相应的计算结果返回到NASTRAN继续进行后续计算；

(4)运行有限元后处理软件，读入有限元分析软件计算的结果，进行可视化处理。通过对该隧道地震响应进行多工况的数值模拟，对隧道地震安全性进行综合评价，并判断其满足7度抗震烈度设防要求。

集成软件一必须判断超级计算机上的人造地震动时程并行计算程序是否完成计算，以将其结果返回到微机。另外考虑到各软件运行于其上的系统平台不同，可采用批处理方式实现集成软件一。该软件运行于后台，用户不需要了解内部流程，即可利用并行计算软件进行人造地震动时程计算。

集成软件二必须判断有限元分析软件是否完成所有的分析计算，以将其结果返回到微机。另外考虑到各软件运行于其上的系统平台不同，可采用批处理方式实现集成软件二。该软件运行于后台，用户不需要了解内部流程，即可利用并行计算软件进行地震响应分析。

集成软件三嵌在有限元求解序列中，其功能是：调用并行计算软件实现位移方程组的求解，并代替NASTRAN的位移方程组求解模块。NASTRAN的后续计算需要并行计算软件的计算结果，因此在并行计算软件运行期间，NASTRAN必须等待。因此，集成软件三的工作原理同集成软件一和二类似，可采用批处理方式实现。

本实施例采用软件，将所有的实现集中到软件的运行中，能大大节省系统实施的人力、物力和财力。