[发明专利]基于时空积分域处理时域边界元奇异积分的方法

摘要

本发明提供了一种基于时空积分域处理时域边界元奇异积分的方法，包括以下步骤首先，通过奇异分离法处理奇异子矩阵元素，采用线性单元进行离散时，将所有变量转变为对的函数，将形函数、，及线元用来表示。本发明的有益效果是对于奇异积分的奇异部分积分的处理不再采用非奇异积分处理时采用的先积时间后空间的积分顺序，变换积分顺序后重新进行积分，计算量大大减少，积分结果也变得简便很多，可提高时域边界元处理弹性动力学问题的计算精度和效率。

主权项

一种基于时空积分域处理时域边界元奇异积分的方法，其特征在于，包括以下步骤：首先，通过奇异分离法处理奇异子矩阵元素，采用线性单元进行离散时，将所有变量转变为对r的函数，将形函数N1、N2，及线元dΓ用dr来表示，则对角线元素为：2πρcs2h‾ik(m;e,1)=∫0Le[(Aik+Dik)esm+Bikdsm-cs2cd2(Bikddm+Dikedm)](1-rLe)dr---(1)]]>——tm瞬时e单元第1点位移对p点位移的影响系数；Le——单元的长度；ρ——材料密度；cs——剪切波波速，cd——压力波波速；μ——剪切模量；λ——拉梅常数；r——源点与场点之间的距离；Bik=-2μr3(δik∂r∂n+r,ink+r,kni-4∂r∂nr,ir,k);]]>ewm=|∫tm-1tmcwLw3Hwdτ‾;]]>dwm=∫tm-1tmcwLwNwHwdτ;]]>Lw=[cw2(t-τ)2-r2]-12;]]>将式(1)分成奇异部分和非奇异部分，表达式如下：2πρcs2h‾sik(m;e,1)=∫0Le[(Aik+Dik)esm+Kkdsm-cs2cd2(Bikddm+Dikedm)]dr---(2)]]>2πρcs2h‾nik(m;e,1)=-∫0Le[(Aik+Dik)esm+Bikdsm-cs2cd2(Bikddm+Dikedm)]rLedr---(3)]]>在计算公式(2)时，建立关于r‑τ的时空坐标系，取r1＝cw(t‑t1)，r2＝cw(t‑t2)，rτw＝cw(t‑τ)，tLw＝t‑r/cw；横轴τ表示脉冲作用于该单元结点的时间，纵轴r表示脉冲作用节点到计算点的距离，直线r＝L表示边界单元的长度，由于只在本单元计算，限定了积分上限，[t1，t2]表示一个时间单元；斜直线r＝cw(t‑τ)表示τ时刻从场点Q发出的脉冲在t时刻将要到达源点P，在源点P存在波前奇异性；斜直线r＝cw(t‑τ)上方表示脉冲未传播到的区域(r>cw(t‑τ)且r<L)，对源点P的响应无任何影响，不予考虑；斜直线r＝cw(t‑τ)下方表示脉冲已传过的区域(r<cw(t‑τ)且r<L)；根据r1、r2与L的相对位置关系，将积分域分成三种可能的积分域进行讨论；当r1≥L且r2≥L时，在r‑τ的坐标系内是矩形域，为第一种积分域；当r1>L且r2<L时，在r‑τ的坐标系内是混合域，为第二种积分域；当r1≤L且r2≤L时，在r‑τ的坐标系内是梯形域，为第三种积分域；由于三种积分域适合先r的积分，后τ的积分，因此与前面非奇异元素求解有所不同；(1)先对r积分对r积分的过程中可能会遇到空间奇异性和波前奇异性；对于第一种积分域，即矩形域，即r2≥L，所有空间积分可以计算如下：∫0LcwLwNwr3dr=dwu-dw0]]>dwu=-cwLw-1r2|Ldwl=-cwLw-1r2|0]]>下标中“u”表示积分上限值，“l”表示积分下限，下同；ew=∫0LcwrLw3dr=cwLw|0L]]>对于第二种积分域，即混合域，即r2≤L且r1≥L。τ∈[t1，tLw)，是矩形域，按矩形域的方法计算；τ∈[tLw，t2]，是梯形域，按梯形域的方法计算；对于第三种积分域，即梯形域，即r1≤L，所有空间积分可以计算如下：dwu=∫0rτwcwLwNwr3dr=-cwLw-1r2|rτw=0]]>ew=limr→cw(t-τ)(∫0rcwrLw3dr-cwLw)=-1t-τ]]>(2)再对τ积分积分过程可能会遇到波前奇异性，如果其中的波前奇异性在对r积分时出现了空间奇异性，那么这个奇异性就是双重奇异性；Aw1=cw(t-t1)-LAw2=cw(t-t1)+LAw5=cw(t-t1)]]>Aw3=cw(t-t2)-LAw4=cw(t-t2)+LAw6=cw(t-t2)]]>(1)矩形域，即r1≥L时，所有时间积分可以计算如下：Idwu1=∫t1t2dwuΨ1dt=∫t1t2-cwcw2(t-τ)2-L2L2t2-τdtdτ=-12γwL2cw2(t-τ)(t-t2)cw2(t-τ)2-L2+13γwL2(cw2(t-τ)2-L2)32-12γwcw(t-t2)ln(cw2(t-τ)2-L2-cw(t-τ))|t1t2=12γwL2-Aw62Aw3Aw4+23(Aw3Aw4)3-L2Aw6ln(Aw6-Aw3Aw4)Aw5Aw6Aw1Aw2-23(Aw1Aw2)3+L2Aw6ln(Aw5-Aw1Aw2)]]>Idwu2=∫t1t2dwuΨ2dτ=12γwL2cw2(t-τ)(t-t1)cw2(t-τ)2-L2-13γwL2(cw2(t-τ)2-L2)32+12γwcw(t-t1)ln(cw2(t-τ)2-L2-cw(t-τ))|t1t2=γw2L2Aw5Aw6Aw3Aw4-23(Aw3Aw4)3+L2Aw5ln(Aw6-Aw3Aw4)-Aw52Aw1Aw2+23(Aw1Aw2)3-L2Aw5ln(Aw5-Aw1Aw2)]]>Iew1=∫t1t2ewΨ1dτ=-cw(t-t2)cwΔtln(cw2(t-τ)2-L2-cw(t-τ))-cw2(t-τ)2-L2cwΔt-cwτ+cw(t-t2)ln(t-τ)cwΔt|τ=t1τ=t2=γw-Aw6ln(Aw3Aw4-Aw6)-Aw3Aw4+Aw6ln(Aw1Aw2-Aw5)+Aw1Aw2-[Aw6γwln(Aw6Aw5)+1]]]>Iew2=∫t1t2ewΨ2dτ=cw(t-t1)cwΔtln(cw2(t-τ)2-L2-cw(t-τ))+cw2(t-τ)2-L2cwΔt-cwτ+cw(t-t1)ln(t-τ)cwΔt|τ=t1τ=t2=Aw5γwln(Aw6-Aw3Aw4)+γwAw3Aw4-Aw5γwln(Aw5-Aw1Aw2)-γwAw1Aw2+[Aw5γwln(Aw6Aw5)+1]]]>(2)混合域积分，需要分两种情况分析；①当τ∈[t1，tLw)时是矩形域，只需将(1)中积分上限t2全部由tLw替换即可；Idwu1=∫t1tLwdwuΨ1dτ=∫t1tLw-cwcw2(t-τ)2-L2L2t2-τdtdτ=-12γwL2[γwAw6lnL-Aw5Aw6Aw1Aw2+23(Aw1Aw2)3-L2Aw6ln(Aw5-Aw1Aw2)]]]>Idw2=∫t1tLwdwΨ2dτ=-γw2L2[-L2Aw5lnL+Aw52Aw1Aw2-23(Aw1Aw2)3+L2Aw5ln(Aw5-Aw1Aw2)]]]>Iew1=∫t1tLwewΨ1dτ=γw[-Aw6lnL+Aw6ln(Aw1Aw2-Aw5)+Aw1Aw2]-[γw(Aw5-L)+Aw6γwln(LAw5)]]]>Iew2=∫t1tLwewΨ2dτ=[Aw5γwlnL-Aw5γwln(Aw1Aw2-Aw5)-γwAw1Aw2]+[γw(Aw5-L)+Aw5γwln(LAw5)]]]>②当τ∈[tLw，t2]时是梯形域，只需将(3)中积分下限值t1全部由tLw替换即可；Idwu1=∫tLwt2dwuΨ1dτ=0]]>Idwu2=∫tLwt2dwuΨ2dτ=0]]>Iew1=∫tLwt2ewΨ1dτ=-[γw(Aw5-L)+Aw6γwln(Aw6L)]]]>Iew2=∫t1t2ewΨ2dτ=Aw5γwln(Aw6L)+γw(Aw5-L)]]>当t2＝t时，Ii，Ij和Ie产生了时间上的奇异性，计算原积分的Hadamad主值：Iew1＝‑γw(Aw5‑L)Iew2=Aw5γwln(cwL)+γw(Aw5-L)]]>(3)梯形域，即r1≤L时，积分计算如下：Idwu1=∫t1t2dwuΨ1dτ=0]]>Idwu2=∫t1t2dwuΨ2dτ=0]]>Iew1=∫t1t2ewΨ1dτ=-cwτ+cw(t-t2)ln(t-τ)cwΔt|τ=t1τ=t2=-[Aw6γwln(Aw6Aw5)+1]]]>Iew2=∫t1t2ewΨ2dτ=cwτ+cw(t-t1)ln(t-τ)cwΔt|τ=t1τ=t2=Aw5γwln(Aw6Aw5)+1]]>当t2＝t时，Ii，Ij和Ie产生了时间上的奇异性，计算原积分的Hadamad主值：Iiw1＝Iew1＝‑1Iiw2＝Iew2＝‑Aw5γwln(t‑t1)+1上述三种情况中，t2≠t时：Idsd1=∫t1t2dsdΨ1dτ=csτ+cs(t-t2)ln(t-τ)2csΔt(1-cs2cd2)|t1t2=12(1-cs2cd2)[1+γsAs6ln(Aw6Aw5)]]]>Idsd2=∫t1t2dsdΨ2dτ=-csτ+cs(t-t1)ln(t-τ)2csΔt(1-cs2cd2)=-12(1-cs2cd2)[1+γsAs5ln(Aw6Aw5)]]]>t2＝t时：Idsd1=limt→0∫t1tdsdΨ1dτ=12(1-cs2cd2)]]>上式求取了Riemann积分，属弱奇异积分；Idsd2=limt→0[∫t1tdsdΨ2dτ+12(1-cs2cd2)γsAs5ln(t-t2)]=-12(1-cs2cd2)[1-γsAs5ln(t-t1)]]]>所有奇异单元的时‑空间积分系数都已求出，将得到的结果直接代入和计算，以提高时域边界元处理弹性动力学问题的计算精度和效率。