一种基于时域边界元法分析半无限域双孔洞动力问题的数值方法与流程

本发明涉及动力学分析领域，特别涉及一种基于时域边界元法分析半无限域双孔洞动力问题的数值方法。

背景技术：

半无限域双孔洞是工程中常见的结构之一，例如，地铁小间距隧道和爆破工程中的炮孔。与单孔洞问题相比而言，影响双孔洞的动力响应的因素更加复杂，其中需要考虑孔洞间的相互作用。因此，往往采用数值方法分析该类问题。目前针对该方法主要采用的是有限元法研究。有限元法在处理半无限域问题时需要建立起人工边界或者粘弹性边界，非常复杂。

有限元法并不能很好地处理无限域或半无限域问题，而是非常粗糙地建立起人工边界或者粘弹性边界，因此会带来数值上的误差或者处理方式非常复杂，难以得到广泛的认可。

与有限元法相比而言，时域边界元法在解决半无限域动力问题上具有明显的优势。然而目前还没有一种针对半无限域双孔洞动力问题的时域边界元算法。其难点主要在于：半无限域的处理；双连通域边界积分的建立和求解。

技术实现要素：

本发明提供一种基于时域边界元法分析半无限域双孔洞动力问题的数值方法，旨在应用于分析列车振动荷载作用下地铁小间距隧道的动力响应及爆炸荷载作用下多个炮孔的响应等常见的实际工程问题。

本发明提供一种基于时域边界元法分析半无限域双孔洞动力问题的数值方法包括以下步骤：

s1.根据半无限域双孔洞问题建立位移边界积分方程、应力边界积分方程；

s2.将所解问题定义域的边界划分成若干个边界元，并引入无穷边界单元，求解出无穷边界单元的插值函数；

s3.根据步骤s2划分的边界元及其单元类型，对边界积分方程进行离散；

s4.将源点p遍历所有边界节点后，基于算法程序求解出边界积分方程的位移和动力。

作为本发明的进一步改进，所述步骤s1中，位移边界积分方程为：

其中，和分别表示单位脉冲荷载在τ时刻作用在源点p的i方向时，t时刻q点k方向产生的位移基本解和面力基本解，uk和pk分别表示k方向产生的位移和面力；

应力边界积分方程为：

其中，为位移影响系数核函数、为面力影响系数核函数。

作为本发明的进一步改进，所述步骤s1中的应力边界方程是根据胡克定律和位移边界积分方程的导数形式建立。

作为本发明的进一步改进，所述步骤s1中，边界积分方程求解的位移和面力可由基本解和迭加求得。

作为本发明的进一步改进，所述步骤s2中，无穷边界单元的插值函数满足以下条件：

(i)

(ii)∑ni＝1；

(iii)当场点q趋向无穷远时，uk(x)→0,pk(x)→0；

其中，i，j表示方向，取值为1或者2，即i,j＝1or2，δij表示狄拉克函数。

作为本发明的进一步改进，所述步骤s2中，无穷边界单元分为左边无穷单元和右边无穷单元；

左边无穷单元的位移和面力的插值函数分别表示为：

左边无穷单元的位移和面力的插值函数分别表示为：

其中，ξ表示高斯积分点，左右1、2、3号位置按照顺时针方向选取，l表示左边无穷单元，r代表右边无穷单元，上标u，p分别表示与位移和面力有关。

作为本发明的进一步改进，所述步骤s3中，对边界积分方程进行离散包括：

将边界离散为n个空间单元，分析并将时间步[0,t]离散为m个时间间隔，每个时间间隔为δt,经过离散边界积分返程中位移和面力表示为：

其中，

表示第m(m＝1,2,3m)时步边界单元e第a点面力对源点mδt时刻位移影响系数；和分别表示第(m)时步和(m+1)时步边界单元e第a点位移对源点mδt时刻位移影响系数；和分别表示边界单元e第a点面力和位移对源点mδt时刻面力影响系数。

作为本发明的进一步改进，所述步骤s3还包括：对离散的单元进行组装，在空间上的组装方式如下：

其中nei(i＝1,2,3)分别指相应边界γ1，γ2和γ3上边界单元的总数，表示第m(m＝1,2,3…m)时步边界单元e面力对源点mδt时刻位移影响系数；表示第m时步边界单元e位移对源点mδt时刻位移影响系数；和分别为边界单元e面力和位移对源点mδt时刻面力影响系数。式中上标“l”表示左边的无限边界元，“r”表示右边的无限边界元。p表示源点，m表示离散时间步长。表达式中cik为与源点位置有关的系数。

作为本发明的进一步改进，所述步骤s4具体包括：

将源点p遍历所有边界节点后，边界积分方程的矩阵形式表达如下：

h^mm{u}^m＝g^mm{p}^m+{a}^m

{σ}^m+s^mm{u}^m＝d^mm{p}^m+{b}^m，并基于matlab将边界积分方程的求解程序化；

其中，h、g分别为和形成的位移和面力总影响系数矩阵，s、d分别为和形成的位移和面力总影响系数矩阵，σ、u、p分别为节点应力，位移和面力向量。

本发明的有益效果是：本发明与有限元法相比而言，所提出的基于时域边界元法分析半无限域双孔洞动力问题的时域边界元法在解决半无限域动力问题上具有明显的优势。采用边界元法分析，可以降低求解问题的维度，只需在边界上划分网格，自动满足无穷远处的边界条件，无需施加人工边界。与其他数值方法相比，本发明提供的方法，显著地减少了计算工作量，提高了计算精度和效率。本发明的方法可广泛应用于地铁小净距隧道的动力分析和多孔洞爆破动力分析等实际工程问题中。

附图说明

图1是本发明中半无限域平行双隧道结构示意图；

图2是本发明中平行双隧道结构的边界离散图；

图3是本发明中左边无穷边界单元示意图；

图4是本发明中右边无穷边界单元示意图；

图5是本发明中隧道的几何性质和荷载分布示意图；

图6是本发明中动力荷载示意图；

图7是本发明图5中a点(2ro,0)处径向应力时程曲线图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

本发明提供一种基于时域边界元法分析半无限域双孔洞动力问题的数值方法，包括如下步骤：

步骤1：如图1所示，根据半无限域双孔洞问题建立相应的位移边界积分方程：

其中：

因此，位移边界积分方程可以简化为以下形式：

其中，和分别表示单位脉冲荷载在τ时刻作用在源点p的i方向时，t时刻q点k方向产生的位移基本解和面力基本解，uk和pk分别表示k方向产生的位移和面力。具体表达式如下：

式中，ρ密度，r是场点q和源点p之间的距离，cs是剪切波或者s波波速，cd为压缩波或者p波波速。基本解中包含参数具体表达式如下：

eik＝δik

hw＝h(mw)

mw＝cw(t-τ)-r

其中h为赫维赛德函数，δik为狄拉克函数，λ，μ为拉梅常数，v为泊松比，n为单元的方向导数，下标i,k,w,v取值为1或者2。

根据胡克定律和位移边界积分方程的导数形式，建立应力边界积分方程：

积分方程中位移影响系数核函数和面力影响系数核函数可通过下式求得：其中，λ和μ为拉梅常数，δij狄拉克函数，分别为面力基本解和位移基本解对相应方向的导数。

内点的应力积分方程可以简写为：

其中和分别代表剪切波和压缩波对与面力和位移有关积分贡献。具体表达式为：

式中张量eijk，fijk，jijk，aijk，bijk，dijk，和具体表达式为：

eijk＝0

其中h为赫维赛德函数，δik为狄拉克函数，λ，μ为拉梅常数，v为泊松比，n为单元的方向导数，下标i,k,w,v取值为1或者2。

步骤2：将所解问题定义域的边界划分成若干个边界元，其中包括引入的无穷边界单元以解决半无限域问题，如图2，图中标注1表示线性边界单元，标注2表示无穷边界单元。

根据积分方程中基本解和的表达式可知，当场点q趋向无穷远处，即r→∞时，基本解和分别为和的同阶无穷小量，如以下关系式：

而实际问题需要求解的位移和面力可由基本解迭加得到，因此，它们在场点q趋向无穷远的衰减速度与基本解同阶。无穷边界单元的插值函数除了满足以上的衰减性质，还需要满足收敛准则，即满足以下条件：

(i)

(ii)∑ni＝1；

(iii)当场点q趋向无穷远时，uk(x)→0，pk(x)→0；

其中，i，j表示方向，取值为1或者2，即i,j＝1or2，δij表示狄拉克函数。

无穷边界单元有两种形式，如图3所示，对于左边无穷单元，位移和面力的插值函数可以分别表示为公式(9)和(10)，如图4所示，右边无穷单元的插值函数分别表示为公式(11)和(12)。

其中、ξ表示高斯积分点，左右1、2、3号位置按照顺时针方向选取，l表示左边无穷单元，r代表右边无穷单元，上标u,p分别表示与位移和面力有关。

步骤3：根据步骤2划分的边界元及其单元类型，对边界积分方程进行离散。为了实现计算批量处理及程序化，需要对离散的单元进行组装。

对边界积分方程进行离散包括：

将边界离散为n个空间单元，分析并将时间步[0,t]离散为m个时间间隔，每个时间间隔为δt,经过离散边界积分返程中位移和面力表示为：

其中，

表示第m(m＝1,2,3…m)时步边界单元e第a点面力对源点mδt时刻位移影响系数；和分别表示第(m)时步和(m+1)时步边界单元e第a点位移对源点mδt时刻位移影响系数；和分别表示边界单元e第a点面力和位移对源点mδt时刻面力影响系数。

在空间上的组装方式如下：

其中nei(i＝1,2,3)分别指相应边界γ1，γ2和γ3上边界单元的总数。

表示第m(m＝1,2,3…m)时步边界单元e面力对源点mδt时刻位移影响系数；表示第m时步边界单元e位移对源点mδt时刻位移影响系数；和分别为边界单元e面力和位移对源点mδt时刻面力影响系数。式中上标“l”表示左边的无限边界元，“r”表示右边的无限边界元。p表示源点，m表示离散时间步长。表达式中cik为与源点位置有关的系数。

步骤4：将源点p遍历所有边界节点后，边界积分方程的矩阵形式可以表达如下：

h^mm{u}^m＝g^mm{p}^m+{a}^m(17)

{σ}^m+s^mm{u}^m＝d^mm{p}^m+{b}^m(18)

其中，

其中，h、g分别为和形成的位移和面力总影响系数矩阵，s、d分别为和形成的位移和面力总影响系数矩阵，σ、u、p分别为节点应力，位移和面力向量。

基于matlab将以上的边界积分方程的求解程序化。

实施例一：

如图5所示，突加荷载作用下半无限域平行隧道的动力响应分析。该问题的材料参数条件如下：材料参数为：弹性模量，e＝8x10⁹pa，泊松比ν＝0.25，密度ρ＝2350kg/m³。隧道的埋深与孔径比值(h/ro)为10，两个平行隧道的间距与孔径的比值(d/ro)为40。突加荷载如图6所示，p＝200mpa。

a点(2ro,0)作为动力响应的计算点，该处径向应力(σr)的时程曲线图如图6所示。

从图7中可以看出，在没有地面反射波和远处隧道压力波的影响，即0s-0.0098s时，该方法的计算结果与解析解的结果(特征性线法)相吻合，证明了该方法计算结果的准确。该处的解析解是通过特征性解法分析无限域单隧道问题求得。

此后，t＝0.0099s时a点受到地面反射波的影响，a点的径向应力值开始与无限域单隧道的解析解不同。而从0.0188s开始，远处隧道产生的压力波也到达a点，a点的径向应力会产生相应的变化。

本发明提供一种基于时域边界元法分析半无限域双孔洞动力问题的数值方法。与其他数值方法相比，本发明提供的方法，显著地减少了计算工作量，提高了计算精度和效率。本发明的方法可广泛应用于地铁小净距隧道的动力分析和多孔洞爆破动力分析等实际工程问题中。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。