一种弹性动力学边界面计算方法

一种弹性动力学边界面计算方法
【专利摘要】本发明涉及一种弹性动力学边界面计算方法，针对传统时域法计算量大、分析时间较长时结果不稳定的问题，提出了一种拟初始条件法和卷积数值积分法相结合的方法来求解三维弹性动力学问题的时域边界积分方程。在新的方法中，首先利用拟初始条件法将整个分析时间划分为几个时间段，将每个时间段末的动态响应作为下一段的初始条件，不必再累加之前分析步的结果对当前分析步的影响，从而可以减少系数矩阵的存储量和计算时间，降低了计算量；然后在每个时间段内，再次划分为更小的时间步采用卷积数值积分法逐步求解。该方法既能够提高传统时域法的稳定性，又克服了卷积数值积分法计算量大的问题。
【专利说明】
一种弹性动力学边界面计算方法
技术领域
[0001] 本发明涉及计算机和力学数值分析技术，具体涉及一种弹性动力学边界面计算方 法。
【背景技术】
[0002] 随着计算机技术的迅速发展，数值计算方法使得我们能够对复杂的工程问题进行 仿真计算，从而得到其数值解。除有限元法以外，边界元法也是常用的数值计算方法之一， 特别是在处理无限域波动、裂纹扩展以及应力集中等问题时边界元法具有天然的优势。然 而，无论是边界元法还是有限元法，都需要单独生成一个离散化的网格近似模型才能进行 分析，存在着几何模型和分析模型不一致的问题。
[0003] 基于边界积分方程和计算机图形学，参考文献《基于边界面法的完整实体应力分 析理论与应用》(作者:张见明，期刊:计算机辅助工程，2010,19(3) :5-10)提出了一种边界 面法。边界面法具备边界元法的所有优点，同时又克服了几何模型和分析模型不一致的问 题。它直接利用几何造型系统中的边界表征数据结构对模型进行离散和分析，能够精确计 算分析过程中所有的几何数据，有效避免了几何误差的产生，是一种适应设计/分析一体化 发展趋势的新型数值分析方法。
[0004] 根据对时间处理方式的不同，弹性动力学问题的边界积分方程可以用以下三种方 法求解:时域法、频域法(Laplace变换法或Four ier变换法)和双重互易法。以下是弹性动力 学问题的时域边界积分方程：
[0005] 对于均匀、各向同性的线弹性材料而言，弹性动力学问题的位移运动方程(亦称为 控制方程)如下：
[0006]
(1)
[0007]式中i，j = l，2,3;Ui代表位移分量，Uj,ij(q，t)为t时亥Ijj方向位移的二阶导数，uMj (q，t)为t时刻i方向位移的二阶导数代表单位质量的体积力;P是材料密度A是位移对 时间的二阶导数，也就是加谏庠:λ和G是拉梅常数，表达式如下：
[0_
(2)
[0009] 分别表示材料的压缩性和剪切模量，E为弹性模量，V为泊松比。
[0010] 若雎件体可弄示为加图1所示的区域V，方程（1)满足边界条件：
[0011] (3)
[0012]
[0013] (4)
[0014] 其中，q表示弹性体上任意一点，可V)和MV)分别表示边界P和於上已知的位移 和面力，?
分别为初始位移和初始速度。
[0015] 控制方程(1)对应的时域边界积分方程为：
[0016]
[0017] 其中》:(P，qU-T)是位移基本解，表示τ瞬时作用在ρ点i方向的单位集中力脉冲所引 起的q点t瞬时j方向的位移分量为是面力基本解； 是速度基本解，为位移基本对时间的一次导数。当源点P位于光滑边界时(^(ρ) = δ^/2,位 于域内时Cij(p) = 3ij。
[0018] 相对而言，时域法是最为自然直接且应用范围最广的一种方法，它可以直接得到 时间相关的动态响应。特别是处理爆炸、冲击等瞬态问题、接触等非线性问题以及几何体形 状随时间变化的动边界问题时，时域法比其他方法更加有利。然而，当分析时间较长时，时 域法可能会出现结果不稳定的现象。Schanz曾将卷积数值积分法应用到时域边界积分方程 中，有效提高了计算结果的稳定性。该方法采用有限项求和近似计算边界积分方程中基本 解和节点值之间的卷积。它采用频域基本解，同时又可以直接得到真正的时域解，不需要将 边界条件进行变换，也不需要将结果进行反变换，这是一种介于时域法和频域法之间的一 种数值计算方法。它继承了传统时域法的优点，同时又能够提高时域法的稳定性。但是也存 在着一个问题，由于卷积数值积分法在计算系数矩阵时是对复数进行操作，又要进行傅里 叶变换，所以其计算量和存储量都比传统时域法大很多。

【发明内容】

[0019] 本发明提供了一种弹性动力学边界面计算方法，以解决现有技术中的弹性动力学 边界面计算方法计算稳定性差且计算量、存储量较大的缺陷。
[0020] 为解决上述问题，本发明的弹性动力学边界面计算方法包括如下步骤：
[0021 ] 1)将所求弹性动力学问题的整个分析时间[0，t ]划分为M个时间段，将对每个时间 段进行弹性动力学边界积分方程的求解过程称为一个子分析，采用拟初始条件法处理；
[0022] 2)将每一个子分析时间段[t^U]进一步划分为N个分析步;根据时域边界积分 方程，利用卷积数值积分法求出第一个子分析的每一个分析步的动态响应及N组权系数矩 阵；
[0023] 3)将第一个子分析的最后一个分析步计算出来的动态响应中的节点位移和速度 作为第二个子分析的虚拟初始位移和速度，调用第一个子分析的N组权系数矩阵，根据时域 边界积分方程，利用卷积数值积分法求出第二个子分析的每一个分析步的动态响应，以此 类推，最终求出M个子分析的每一个分析步的动态响应，实现弹性动力学边界面计算。
[0024] 在利用卷积数值积分法进行计算前，将当前子分析的节点位移和速度分别转换为 单位质量初始位移的等效虚拟力和单位质量初始速度的虚拟力，利用单位质量初始位移的 等效虚拟力和单位质量初始速度的虚拟力，将时域边界积分方程转换为初始位移加上基于 位移增量的边界积分方程。
[0025] 所述时域边界积分方程为：
[0026]
[0027] 其中，<(P，q;z-『)是位移基本解，表示τ瞬时作用在P点i方向的单位集中力脉冲所引起 的q点t瞬时j方向的位移分量为< O是面力基本解；〖〇似-〇) = M(P，叫_〇)?， 是速度基本解，为位移基本解对时间的一次导数。
[0028] 所怵拟系教的i+笪公忒加下，
[0029]
[0030]
[0031] 其中，4
为线性多步法中的二次向后差分公式， 为误差控制参数，计算公式为 分别为拉普拉斯变换后的位移基本解和面力基本解。
[0032] 所述基于位移增量的边界积分方程为：
[0033] ?=()
[0034] 其中
为权系数；APj(q，tm-i+aAt) 为边界S上任意一点的面力增量，Auj(q，tm-i+aA t)为边界S上任意一点的位移增量，bj(q, tm-1+a At)为单位质量的体积力，〇)为单位质量初始位移的等效虚拟力，广(q)为单位 质量初始速度的等效虚拟力，P是材料密度。
[0035]所述单位质量初始位移的等效虚拟力的公式为：
[0036]
[0037] 分别表示材料的压缩性和剪 切模量，E为弹性模量，V为泊松比，H( t-Ο)为Heavi s ide函数，P是材料密度，Uj, ij (q，0)为初 始时亥?方向位移的二阶导数，m,ij(q，0)为初始时刻i方向位移的二阶导数。
[0038] 所述单位质量初始速度的虚拟力的公式为：= vIl(tI)冲)，其中，δ(?)是Dirac delta函数，r〉'(q)为q点在i方向的初始速度。
[0039] 本发明的有益效果:针对传统时域法计算量大、分析时间较长时结果不稳定的问 题，本发明提出了一种拟初始条件法和卷积数值积分法相结合的方法来求解三维弹性动力 学问题的时域边界积分方程。在新的方法中，首先利用拟初始条件法将整个分析时间划分 为几个时间段，将每个时间段末的动态响应作为下一段的初始条件，不必再累加之前分析 步的结果对当前分析步的影响，从而可以减少系数矩阵的存储量和计算时间，降低了计算 量;然后在每个时间段内，再次划分为更小的时间步采用卷积数值积分法逐步求解，该方法 既能够提高传统时域法的稳定性，又克服了卷积数值积分法计算量大的问题。
[0040] 无论是拟初始条件法引入的虚拟初始条件，还是工程问题中实际存在的初始位移 和初始速度，都在边界积分方程中涉及到体积分的计算。由于弹性动力学问题的时域基本 解含有Dirac Delta函数，而且边界积分方程中的初始条件项不含有卷积，无法使用卷积数 值积分法，所以导致初始条件项的数值体积分变得非常复杂。针对该问题，本发明提出了一 种适用于三维线弹性动力学问题的虚拟力法，基于线性系统的叠加原理以及冲量定理将初 始条件转化为等效的体积力来处理，原控制方程被转换为一个初始条件完全为零的控制方 程，从而避免了复杂的几何求交操作以及球面积分的产生，大大简化了初始条件项的计算。
【附图说明】
[0041 ]图1为弹性体边界条件示意图；
[0042]图2为弹性动力学边界面计算方法流程图；
[0043] 图3为悬臂杆几何模型和边界条件示意图；
[0044] 图4(a)为外力作用下悬臂杆自由端长时间的位移响应结果对比整体曲线图；
[0045] 图4(b)为外力作用下悬臂杆自由端长时间的位移响应结果对比局部放大曲线图； [0046]图5(a)为外力作用下悬臂杆固定端长时间的面力响应结果对比整体曲线图；
[0047]图5(b)为外力作用下悬臂杆固定端长时间的面力响应结果对比局部放大曲线图；
[0048]图6为计算时间随分析步数的变化对比图。
【具体实施方式】
[0049]下面结合附图，对本发明的技术方案作进一步详细说明。
[0050]如图2所示，本实施例的弹性动力学边界面计算方法的步骤如下：
[0051 ] 1)将所求弹性动力学问题的整个分析时间[0，t ]划分为M个时间段，将对每个时间 段进行弹性动力学边界积分方程的求解过程称为一个子分析；
[0052] 2)将每一个子分析时间段[t^U]进一步划分为N个分析步;根据时域边界积分 方程，利用卷积数值积分法求出第一个子分析的每一个分析步的动态响应及N组权系数矩 阵；
[0053] 3)将第一个子分析的最后一个分析步计算出来的动态响应中的节点位移和速度 作为第二个子分析的虚拟初始位移和速度，调用第一个子分析的N组权系数矩阵，根据时域 边界积分方程，利用卷积数值积分法求出第二个子分析的每一个分析步的动态响应，以此 类推，最终求出M个子分析的每一个分析步的动态响应，实现弹性动力学边界面计算。
[0054] 下面详细阐述上述步骤：
[0055] 首先，利用拟初始条件法，将整个分析时间[0，t]划分为M个时间段，称每个时间段 的分析为子分析。在每个子分析内，将上一个子分析的计算结果作为当前子分析的初始条 件，那么当前子分析相应的时域边界积分方程为：
[0056]
[0057] 和速度。
[0058]然后将每一个子分析的时间段[t^U]进一步划分为N个更小的分析步，时间步 长设为△ t，用卷积数值积分法求出每一个子分析步的动态响应，时间1?散后的方程为：
[0059:
[0060：
[0061：
[0062:
[0063]其中，γ (Z) = 1.5-2z+0.5z2为线性多步法中的二次向后差分公式；当L = N+1以及 纪=石时，权系数ω n-a的误差量级为〇(')（1〇3 SidO111)，《(p，q.i)和Κ(Μ，?)分别是拉普拉 斯变换后的位移基本解和面力基本解。其中，权系数ω "1的值只与分析步之间的时间差n-a 的大小有关，与第几个子分析m无关，因此只需在第一个子分析中将N组系数矩阵计算出来 并存储，之后的所有子分析计算时只要调用即可，不需要重新计算。系数矩阵的存储量由原 来的M X N组降为N组，相应的计算时间也能大幅减少。
[0064] 方程(7)等号右边的后两项代表的是在上一个子分析的最后一步计算出来的节点 位移和速度，我们将其作为当前子分析的虚拟初始位移和初始速度。这两项体积分不包括 时间卷积，因此不能使用卷积数值积分法进行计算。为此，我们提出了一种适用于三维弹性 动力学问题的虚拟力法，将初始条件转换为等效的体积力处理，具体过程为：
[0065] 首先进行初始位移的等效。对于线弹性系统，可以设控制方程（1)的解等于初始位 移加上位移增量，如
[0066] Ui(q,t)=Ui(q,0)+Aui(q,t) (10)
[0067] 将公式（10)代入方程(1)可得：
[0068]
[0069] 定义单位质量初始位移的等效虚拟力为：
[0070]
(12)
[0071]之所以定义为单位质量的虚拟力，是为了与控制方程中的体积力单位保持一致。 这里H(t-〇)为Heaviside函数，只有当t>0的时候才有值。当t>0时，可以将上式简写为：
[0072]
Π3)
[0073] 若要求出域内任意一点初始位移虚拟力fWq)的值，可以通过以下两种方法：
[0074] (1)当m(q，0)的函数表达式已知时，可以直接对m(q，0)进行求导，然后求出相应 的fiu°(q)的值。
[0075] (2)当m(q，0)无法用具体的函数表达时，也就是说初始位移的分布是任意的、没 有规律的，那么可以通过求解下列Navier方程来得到fWq)的值。
[0076] 上述方程为弹性静力学控制方 程，与传统的控制方程求解不一样，在这个方程中，初始位移m(q，0)为已知量初始位移，而 体积力fWd)为需要求解的未知变量。
[0077] 接下来进行初始速度的等效。由动量定理"物体动量的增量等于它所受合外力的 冲量"可知，从t = (T到t = 0时刻物体动量的变化可以等效为单位质量的冲击力fV^qA)作 用在物体上，即
[0078]
(15}
[0079] 其中，衫(!,『) =·()。也就是说物体一开始处于静止状态，直到冲击力作用后才产生了 初始速度。由于作用时间很短，物体的位移还没来得及变化。将初始速度~他〇) = ^如)代入 方程(15)可以得到单位质量初始速度的虚拟力为
[0080]
(16)
[0081] 其中δ(?)是Dirac delta函数，代表./TdO除了t = 0时其他时刻值都为零。如果我 们将冲击力的作用平均到第一个时间步[0, At]内，就可以得到不含时间变量的初始速度 虚拟力^
[0082] (17)
[0083]最后，将初始位移和初始速度的等效虚拟力（12)和（16)代入原始控制方程(1)，可 以得到一个新的关于位移增量的控制方程：
[0089]原始控制方程（1)的解就等于新的控制方程（18)的解△ Ui(q，t)加上初始位移Ui
[0084]
[0085]
[0086]
[0087]
[0088] (q，〇) ο
[0090] 将上述虚拟力法应用于方程(7)，则方程（7)的解等于初始位移uKpAh)加上下 面关于位移增量的边界积分方程的解：
[0091]
[0092]
[0093]
[0094]
[0095] 本发明使用边界面法进行空间变量的离散和积分，边界面法是直接基于边界表征 实体造型数据结构实现其网格离散及单元积分的，实体的边界曲面都以参数形式表达。在 其实现中，无论是边界的数值积分还是场变量的插值都是在边界曲面的二维参数空间里进 行。该方法首先在曲面的二维参数空间里将每个参数曲面离散成若干个参数曲面单元，这 些单元同样都是定义在二维参数空间，而非三维物理空间，也就相当于CAD造型系统中的分 片曲面，从而曲面原始的几何信息都能够保留下来。积分过程使用到的几何数据，比如高斯 积分点的三维坐标、雅可比、外法线向量都是直接由参数曲面计算获得，而不是通过分段多 项式插值近似，从而避免了几何误差。
[0096] 将物体表面离散为BN个表面节点，DN个域内节点。对于任意一个源点qk，（ k = 1， 2,--,BN或者k = l .2. ···.),空间富散后方趕（21 )可弄示为加下钽陈方趕，
[0097]
[0098] 其中m=1，2,…，Μ，η = 1，2，···，Ν.系数矩阵只需要在第一个子分析中求出 即可，只要离散的分析步长不变，就不需要重新计算系数矩阵。求出当前子分析中每一个分 析步每个节点的位移增量'H1 +?△〇后，加上初始位移Uj (q，tm-i)的值就是每个节点真 正的位移响应。当一个子分析求解结束后，还需要求解出最后一个分析步所有域内节点的 位移和速度，为下一个子分析做准备。初始位移的计算直接求解域内边界积分方程即可，在 域内边界积分方程中Cij (p) = Sij(当i = j时，Sij = 1，否则Sij = 〇)，源点p取遍所有域内节点。 节点速度的计算需要用到最后一步以及前一步的节点位移，设在每一个分析步内位移呈线
f生变化，抓Λ晶臼一来的书占迪IiPff曾公才加下
[0099] (25)
[0100] 重复以上过程计算下一个子分析，直到M个子分析全部结束。加入拟初始条件法之 后，系数矩阵的存储量降为原来的1/M，虽然引入了体积分，但是计算时间仍然比之前要少， 划分的子分析数越多，内存需求和计算时间就越短。但是子分析数也不能过多，否则会影响 计算精度。
[0101]下面通过一个算例验证上述方法的正确性，对拟初始条件法的计算效率和计算精 度进行研究。所使用微机的处理器为3.4GHZ的Intel(R)C〇re(TM)i7-2600 CPU，内存空间为 8GB〇
[0102] 计算悬臂杆的受迫振动，悬臂杆几何模型及尺寸如图3所示，长L = 8m，宽和高均为 2m。左端面固定，右端面沿X轴承受一 Heavi side类型的均布拉力p = 1000N/m2。材料参数为： 弹性模量E = I. I e5N/m2，密度P = 2. Okg/m3，泊松比υ = 〇. 〇。本算例暂不考虑重力的影响，之 所以取零泊松比，是为了将数值结果和解析解进行比较，在程序中仍然将该问题当作三维 问题来计算。此工况下，悬臂杆上任意一点任意时刻的位移响应解析解可以用下面的公式 计算：
[0103]
[0104]
[0105]
[0106]
[0107] 在这个算例中，采用线性四边形单元和线性六面体单元离散模型，网格的单元尺 寸均为lm。加入拟初始条件法之前，分析步数如果超过400步，就会出现计算机内存不足的 情况。加入拟初始条件法之后，这里将分析时长延伸到1.424s，时间步长仍然为0.00089s， 总的分析步数为1600步，下面来验证本发明方法的性能：
[0108] 首先将整个分析过程划分为M个子分析，在每个子分析中，分析步个数为1600/M。 为了考察不同M的取值对结果的影响，分别取M=IO和M=20的数值结果与解析解进行对比。 当然，我们也考察了很多其它取值，计算结果相差并不特别明显，为了图片能够显示的清 楚，只取其中的两组值来观察。
[0109] 同时，我们也实现了Schanz的Reformulated CQM(改写的卷积数值积分法），在 Reformulated CQM中，只需开辟一组复数系数矩阵的存储空间，但是每一步需要重新计算 一次单元积分，未知量的系数矩阵也需要重新求逆。这个方法和频域法更加相似，相对来说 Reformulated CQM的计算量比传统的卷积数值积分法要小，但是却丧失了时域法的主要优 势。为了对比加入拟初始条件法前后卷积数值积分法的计算结果，我们将Reformulated CQM的结果也加入到对比图片中。
[0110]图4显示了 1.424s内悬臂杆自由端位移随时间的变化，图中曲线1代表一维杆纵向 振动的解析解，曲线2代表Reformulated CQM的数值解，曲线3、4代表卷积数值积分法和拟 初始条件法相结合的数值解（convolution quadrature method and pseudo-initial condition method，缩写为CQM&PICM)。可以看到，完全理想的工况下解析解一直保持同样 的振幅和频率，三种数值解都是稳定的，和解析解吻合较好，但是均有不同程度的衰减。为 了更清楚的观察数值解之间的差别，我们将前1/6和后1/6时间段的结果放大，图4(b)为图4 (a)的部分截断放大图，可以看到:在分析的一开始，三种数值解的精度都较高，随着分析时 间的增加，结果误差越来越大;相对而言，Reformulated CQM的数值解精度高于后两种数值 解，也就是说拟初始条件法的加入会对卷积数值积分法的结果精度有所影响，因为在虚拟 力法的应用中会引入额外的误差;子分析个数M=IO的数值解比M = 20的数值解精度高，可 见划分的子分析个数越少，对卷积数值积分法的结果精度影响就越小。同样的现象在固定 端面力的动态响应结果中也可以观察得到，如图5所示。
[0111]接下来对比三种数值解法的计算时间。本算例中统计的均为直接计算的时间，并 未加入快速算法。图6显示了三种数值解法的计算时间随分析步数的变化对比， Reformulated CQM随分析步数的增加呈线性增加趋势，当分析步数继续增加时，该方法的 计算时间将会非常可观;拟初始条件法中M=IO的计算时间为M = 20的两倍，因为前者需要 计算并存储的系数矩阵为后者的两倍。当总的分析步数较少时，本算例中拟初始条件法的 计算时间比Reformulated CQM的计算时间长，这是因为在这个例子中，初始条件和体积力 均为零，Reformulated CQM的边界积分方程中不含有体积分项，而拟初始条件法的加入需 要增加体积分的计算，增加了系数矩阵的计算时间。但是随着分析步数的增加，拟初始条件 法的计算时间增加非常缓慢，可见该方法能够显著减少卷积数值积分法的计算时间，提高 弹性动力学分析的计算效率。
[0112]表1对比列出了三种数值解法中系数矩阵所需的内存空间。因为Reformulated CQM只需要存储一组系数矩阵，所以其所需的内存空间最小，但是该方法不能算是真正的时 域法，因为需要对边界条件和计算结果进行傅里叶变换，只能用于边界条件随时间变化已 知的问题。
[0113]表1系数矩阵所需的存储空间对比
[0115]从以上算例分析可知，当分析步数较多时，拟初始条件法能够显著降低传统时域 法的计算时间和系数矩阵的存储量，与卷积数值积分法相结合能够分析较长时间的振动响 应。在拟初始条件法中划分的子分析数越多，计算效率越高，但是结果精度会有所降低。所 以，我们应该根据实际情况，如分析时长、节点个数、计算机能力等合理使用拟初始条件法， 划分适当的子分析数来平衡计算时间和计算精度之间的关系。
【主权项】
1. 一种弹性动力学边界面计算方法，其特征在于，该方法包括如下步骤： 1) 将所求弹性动力学问题的整个分析时间[〇，t]划分为M个时间段，将对每个时间段进 行弹性动力学边界积分方程的求解过程称为一个子分析，采用拟初始条件法处理； 2) 将每一个子分析时间段[t^U]进一步划分为N个分析步;根据时域边界积分方程， 利用卷积数值积分法求出第一个子分析的每一个分析步的动态响应及N组权系数矩阵； 3) 将第一个子分析的最后一个分析步计算出来的动态响应中的节点位移和速度作为 第二个子分析的虚拟初始位移和速度，调用第一个子分析的N组权系数矩阵，根据时域边界 积分方程，利用卷积数值积分法求出第二个子分析的每一个分析步的动态响应，以此类推， 最终求出M个子分析的每一个分析步的动态响应，实现弹性动力学边界面计算。2. 根据权利要求1所述弹性动力学边界面计算方法，其特征在于，在利用卷积数值积分 法进行计算前，将当前子分析的节点位移和速度分别转换为单位质量初始位移的等效虚拟 力和单位质量初始速度的虚拟力，利用单位质量初始位移的等效虚拟力和单位质量初始速 度的虚拟力，将时域边界积分方程转换为初始位移加上基于位移增量的边界积分方程。3. 根据权利要求1所述弹性动力学边界面计算方法，其特征在于，所述时域边界积分方 程为：其中，当源点P位于光滑边界时Ci j (p) = 5i」/2，位于域内时Ci j (p) = 5i」，当i = j时，& j = 1，否则 Sij = 〇;<(P，叫-0是位移基本解，表示谢时作用在p点i方向的单位集中力脉冲所引起的q点t 瞬时j方向的位移分量为< 叫-r)是面力基本解〇)/改，是 速度基本解，为位移基本解对时间的一次导数。4. 根据权利要求1所述弹性动力学边界面计算方法，其特征在于，所述权系数的计算公 式如下：Y (幻=1.5-22+0.522为线性多步法中的二次向后差分公式，货为误差控制参数，计算公式为： 为拉普拉斯变换后的位移基本解和面力基本解。5. 根据权利要求2所述弹性动力学边界面计算方法，其特征在于，所述基于位移增量的 边界积分方程为：其中，％ (化]],A/)、％ "(A;A/)及叫(％.p,A/)为权系数；A pj(q，tm-i+a S上任意一点的面力增量，A uj(q，tm-l+a A t)为边界S上任意一点的位移增量，bj(q，tm-l+a A t)为单位质量的体积力，flq)为单位质量初始位移的等效虚拟力，//Iq)为单位质量初 始速度的等效虚拟力，P是材料密度。6. 根据权利要求2所述弹性动力学边界面计算方法，其特征在于，所述单位质量初始位 移的等效虚拟力的公式为，分别表示材料的压缩性和剪切模 量，E为弹性模量，v为泊松比，H(t-O)为Heaviside函数，P是材料密度，Uj,ij(q，0)为初始时 亥Ijj方向位移的二阶导数，m,ij(q，0)为初始时亥Iji方向位移的二阶导数。7. 根据权利要求2所述弹性动力学边界面计算方法，其特征在于，所述单位质量初始速 度的虚拟力的公式为利0，其中，S(t)是Dirac delta函数，为q点在i方 向的初始速度。
【文档编号】G06F17/50GK106055743SQ201610340671
【公开日】2016年10月26日
【申请日】2016年5月20日
【发明人】李源, 孙林, 张见明, 王世勋, 张仕光, 王艳军
【申请人】河南师范大学