边界积分法（精选4篇）
边界积分法 篇1
1 引 言
有限元法是近似求解数理边值问题的一种数值方法,他于20世纪40年代提出,随着计算机技术的迅速发展,已被广泛地应用在结构分析、机械制造、建筑设计、电磁分析等领域。由于他对场域的剖分具有很大的灵活性,因此对不规则、非均匀的物体具有广泛的适用性,而且形成的矩阵为稀疏矩阵,人们可以利用这一特点快速求解含未知量的稀疏矩阵方程[1,2]。但有限元在分析开域问题时遇到困难,因为区域的无限大,有限元剖分时将有无穷多个节点,所形成的矩阵方程将很难求解[3]。多年来对于二维开域问题的研究概括起来主要分为两大类[4,5,6,7,8]:一类是利用吸收边界条件或完全匹配层技术将无限大空间截断,使原散射问题限制在被研究物体周围一个有限的范围内,然后采用有限元进行分析;第二类是应用矩量法或边界积分方法来处理柱体外的无限大空间,柱体所在区域采用有限元法来分析。第一种方法缺点在于扩大了研究区域,产生了较多的节点数,特别是吸收边界的形状以及他与散射体间的距离对于解的精度有很大影响。第二种方法精度较高,产生的节点数较少,但生成的是部分稀疏,部分满秩的矩阵,如果能够快速求解该方程,则该方法有着较好的应用前景,特别对于电尺寸较大的散射体。
本文在已有文献工作的基础上,应用高阶有限元-边界积分法分析二维散射体的电磁散射特性。对散射体内、外区域分别应用高阶有限元和边界积分法进行分析。然后通过场的连续性进行耦合,形成待求矩阵方程组,求解该方程组,获得每个结点的场值,进一步求解雷达散射截面。作为算例,分别应用二阶和三阶有限元-边界积分法计算了几种二维方柱(导体和介质)的雷达散射截面,并将结果与矩量法的数值结果作了误差分析。数值结果表明,高阶有限元-边界积分法较之一阶线性有限元-边界积分法有着更高的计算精度、收敛速度和更高的计算效率。
2 理论分析
2.1 高阶插值基函数
传统的线性(一阶)有限元-边界积分法的缺点在于精度较低,对一定的单元或结点数,解的收敛性较慢。获得较高精度而不增加结点数的一个有效方法是应用高阶插值函数,下面将讨论以Lagrange插值法为基础的高阶插值函数的建立过程。
首先介绍高阶三角形单元插值函数的建立。Lagrange的k次多项式插值需要undefined项,需要undefined个结点。按照三角形各边的k等分点以及他们连线的交点作为结点,可以建立如图1所示的各阶三角形插值单元。设三角形单元e=ΔP1P2P3的面积为Δe,结点Pi的坐标为(xi,yi),对于e内任意一点P(x,y),定义他的面积坐标为:
undefined
由于三角形面积
undefined
,所以本文所要讨论的二阶和三阶三角形单元各个结点的面积坐标可以由以上定义的公式依次得到。具体见图1所示(为了讨论方便,图中所有面积坐标均取为整数)。
为了建立插值函数,可以用三个整数(I,J,K)来表示单元各个结点的面积坐标。那么插值函数可以写成[9]:
undefined
其中I+J+K=n,n为插值阶数,undefined。同样可以定义Pundefined(L2)和Pundefined(L3) 。对于式(1)可以证明:在结点i上,Nundefined=1;而在单元其他结点上,Nundefined=0。有了式(1)和以上讨论的结果,可以建立任意阶的三角形单元插值函数。对于二阶单元:
undefined
对于三阶单元有:
undefined
由此在每个高阶三角形面单元内,场可以表达为:undefined,其中对二阶单元,n=6;对三阶单元,n=10。对于线段的高阶插值函数可简单的按照Lagrange插值法建立,具体参见文献[10]。
二阶单元有3个结点:
undefined
三阶单元有4个结点:
undefined
由此每个线段内的场由高阶插值函数可以表达为:undefined。其中对二阶单元,m=3;对三阶单元,m=4。
2.2 有限元分析
对散射体所在区域用有限元进行分析。对Ez和Hz极化,散射体所在区域的场分别满足的Helmholtz方程可以统一表达为:
undefined
其中,对Ez极化,u=μr,υ=εr;对Hz极化,u=εr,υ=μr。Ω表示介质柱所在区域。在介质柱的表面上,假设场的法向导数为:
undefined
根据变分原理[11],与式(2),式(3)等价的变分表达式为:
undefined
其中:
undefined
式中Γ表示虚构边界。
为了离散式(4),散射体区域被高阶三角形单元划分为M个小单元(对于导体来说,需要在离导体目标0.1个波长到0.2个波长的外围区域建立新的虚构边界,离散区域是新的虚构边界到导体边界的区域)。相应地,虚构边界Γ被分成Ms个小线段。对每个三角形单元和线段进行高阶基函数插值后代入式(4),并利用变分法,再对所有矩阵单元进行组合,可得到下面的矩阵方程:
undefined
上式中矩阵[K]和[C]的具体计算公式参见文献[12],其中[K]是N×N的对称、稀疏、正定方阵,[C]是N×Ms矩阵,{Φ}表示N个剖分结点处电场或磁场的列向量,{Ψ}是表示在Ms个边界结点处Ψ值的列向量。
2.3 边界积分方程分析
对散射体外的无限大区域用边界积分方程分析。总场Φ满足Helmholtz方程,即:
undefined
其中,f(undefined)为产生入射场的源,Ω∞表示散射体外无限大区域。应用第二格林定理,并引入二维自由空间格林函数undefined,方程(6)变为下列形式:
undefined
在散射体虚构边界上,场Φ应满足下面的边界条件:
undefined
其中Γ-表示观察点从边界内接近Γ,Γ+表示观察点从边界外接近Γ。由于边界Γ被剖分为Ms个线段,对每个线段进行高阶插值,并在等式两边同时在Γ上进行积分,对矩阵单元进行组合可以得到下面的矩阵方程:
undefined
上式中矩阵[P′],[Q],{b}的具体计算公式参见文献[12],显然方程组(9)含有Ms个方程。方程组(5)和(9)可以构成一个完备的方程组,是一个稀疏度较低的矩阵方程。通过求解这个完备的方程组可以得到每个结点的Φ值和虚构面上每个结点的Ψ值。进一步雷达散射截面的计算公式可参考文献[12]。
3 计算实例
算例1 为了验证高阶有限元-边界积分法理论的正确性和优越性,首先计算边长为2个波长的导体方柱对EZ极化平面电磁波的散射(入射场为e-jk0x)。计算过程中一阶方法,剖分精度为λ/20,二阶方法,剖分精度为λ/10,三阶方法,剖分精度为λ/6,这样未知量个数基本相同。图2给出了一阶、二阶、三阶算法的RCS与矩量法数值结果的对比。图3给出了一阶、二阶、三阶算法的RCS与矩量法数值结果误差的比较。
表1给出计算导体方柱RCS在相同误差精度条件下,3种方法所需的未知量个数、占用的内存比较。通过比较可以清楚的看到,在相同误差精度条件下,二阶方法比一阶方法节省了约89%的内存,三阶方法更是节省了约96%的内存。本文中的全局误差由公式计算得出:
undefined
算例2 计算边长为1个波长的介质方柱对EZ极化平面电磁波的散射(介电常数为2-j,入射场为e-jk0x)。剖分精度与算例1相同。图4给出了一阶、二阶、三阶算法的RCS与矩量法数值结果的对比。图5给出了一阶、二阶、三阶算法的RCS与矩量法数值结果误差的比较。
表2给出在需要相同求解未知量个数条件下,3种方法的全局误差比较和CPU时间比较。通过比较可以看到,在相同求解未知量个数条件下,二阶方法节省了约57%的时间,且精度高。三阶方法则给出了更高的精度,节省了约75%的计算时间。
算例3 计算边长为5个波长的导体方柱对EZ极化平面电磁波的散射(入射场e-jk0x)。本算例采用三阶基函数方法,剖分精度为每个波长剖分6段,有3 721个未知量,耗时7.3 s。因为计算机内存的限制,一阶方法无法快速有效地计算电尺寸稍大的目标。本算例充分体现了高阶基函数的优越性。图6是5个波长的导体方柱的RCS,并与矩量法的结果做了比较。
所有数据均在P4 3.0 GHz,RAM为1 GB的PC机,Matlab环境下运行得出。
4 结 语
采用高阶插值基函数离散有限元方程和边界积分方程,构成高阶有限元-边界积分法,计算了典型二维散射体的雷达散射截面, 与矩量法的数值结果比较和一系列的误差分析表明了高阶有限元-边界积分法较之一阶方法有更高的计算精度、收敛速度和计算效率。由于高阶方法允许更大的剖分网格,因此在分析电尺寸较大的散射体时有高效快速的优势。
摘要：以二阶、三阶基函数为例,应用高阶有限元-边界积分法分析了二维散射体电磁散射特性。计算了几种二维方柱(导体和介质)的雷达散射截面,结果与矩量法一致,对三种数值结果进行了误差分析。结果表明:高阶有限元-边界积分法比一阶有限元-边界积分法有着更高的计算精度、收敛速度和计算效率。
关键词：高阶有限元,边界积分法,二维散射体,雷达散射截面
边界积分法 篇2
建立平面弹性薄板弯曲问题理论中具有直接变量的等价边界积分方程.传统的直接变量边界积分方程,它们都不是等价的,对此进行了深入的`讨论.
作 者：张耀明 张作泉 孙焕纯 吕和祥 作者单位：张耀明(南京航空航天大学,能源与动力系,江苏,南京,210016)
张作泉(北方交通大学,数理学院,北京,100044)
孙焕纯,吕和祥(大连理工大学,力学系,大连,116023)
边界积分法 篇3
比例边界有限元法(scaled boundary finite element method,SBFEM)[10,11]是由Wolf和Song于20世纪90年代提出的一种半解析算法.该方法首先将有限域边界沿径向外层克隆一个内外边界相似的边界单元,然后基于相似中心和比例边界坐标变换,将标准坐标下的无限域波动方程变换为比例边界坐标方程,最后采用与有限单元法相同的加权余量法或虚功原理建立半解析的比例边界有限元方程以及动力刚度方程.比例边界有限元法结合了有限元、边界元法的优点,边界按有限元进行离散,分析问题的维数降低一维,并自动满足无限远处的辐射条件,在解决无限域问题方面取得了良好效果.
本文基于比例边界有限元法在文献[12-13]基础上实现了一种局部的高阶透射边界,其中采用改进的连分式法求解无限域的动力刚度矩阵.然后,讨论了该高阶边界的稳定性问题及处理措施.最后,通过两个数值算例说明了该高阶边界的精确性.
1 比例边界有限元法位移控制方程
根据弹性动力学方程和比例边界坐标变换,可以推导出二维、三维的SBFEM位移控制方程[10,11]为
式中,s=2,3分别表示二维、三维情况;u(ξ)为比例边界坐标系中的位移向量;系数矩阵表示为
式中,B1,B2为应变-位移关系矩阵;D为弹性矩阵;J为边界上的雅可比矩阵;N为形函数;ρ为质量密度.
式(2)中的系数矩阵与ξ无关,它们只依赖于子域的几何形状和材料参数.E0,E1和E2类似有限元法中的刚度矩阵,M0类似质量矩阵.可以证明[10,11],E0和M0是正定矩阵,E2是半正定矩阵,E1是非对称矩阵.
2 高阶透射边界的建立
2.1 高阶透射边界
采用无限域动力刚度矩阵S∞(ω)表示的比例边界有限元方程[10,11]为
式(3)在整体上是严格的,是以动力刚度矩阵S∞(ω)为未知量关于ω的一阶常微分方程,可以采用龙格-库塔法[10]在频域内求出其半解析解.
无限域动力刚度矩阵S∞(ω)采用以下的连分式展开进行求解,即
式中,K∞和C∞为无限域的常数刚度矩阵和阻尼矩阵;Y0(i)和Y1(i)为常数矩阵;i=1,2,···,Mcf,Mcf为无限域的连分式展开阶数,当逐渐增大该阶数时,该算法可以在全频域范围内快速逼近准确解.
将式(4)代入式(3)中展开,关于(iω)0,iω和(iω)2项的系数矩阵为0时等式才成立,可分别求得C∞和K∞;然后结合式(5)建立递推关系,可依次求得Y1(i)和Y0(i),其计算过程详见文献[12-13].最初采用连分式法求解动力刚度矩阵时,选取Xu(i)=I,I为单位矩阵.对于自由度较少并且连分式阶数较低时,该算法能保持稳定;但对于自由度较多的系统当连分式阶数逐渐增大时,该算法可能会造成矩阵病态,求解结果甚至是错误的.Xu(i)的引入是为了增加系统的稳定性,因此称为改进的连分式法,采用LDLT分解[13]求解Xu(i).
根据边界ξ=1上的力-位移关系和引入辅助变量,得到一个表示外力载荷幅值与边界位移幅值关系的频域表达式[13],然后变换到时域建立了无限域时域分析的高阶透射边界,即
其中,系数矩阵Ku和Cu,待求向量z(t)以及右端载荷向量f(t)分别表示为
式中,Ku是无限域广义的刚度矩阵,为三块对角矩阵,主对角块由K∞和Y0(i)组成,次对角块由-Xu(i)及其转置矩阵填充.Cu是无限域广义的阻尼矩阵,为由C∞和Y1(i)组成的块对角矩阵.若结构自由度大小为nd,则Ku和Cu矩阵维数为[(Mcf+1)nd]×[(Mcf+1)nd].同时,由于系数矩阵K∞,C∞,Y0(i)和Y1(i)都是对称的,所以Ku和Cu是对称并且是稀疏的.实际编程计算时,根据稀疏矩阵性质,可以减小矩阵存储量,有效地提高计算效率.ub为高阶边界的位移向量,v(i)为辅助变量,R(t)为载荷矩阵.
2.2 稳定性措施
式(6)表示为一阶常微分方程组,其稳定性取决于其系数矩阵的广义特征值问题,即
其中,λ为矩阵B和A的广义特征值,λ=λr+iλi,λr和λi分别为特征值的实部和虚部;为对应的特征向量.
若λr0,式(8)中的响应z呈指数级衰减,表明系统是稳定的;反之,若λr>0,响应z则呈指数级增长,计算结果发散,表明系统存在虚假模态,是不稳定的.因此,必须采取有效措施来保证系统稳定.本文采用的是移谱法(spectral shifting technique)[14,15],分两种情况进行讨论:(1)特征值λj为正实数;(2)特征值λj为实部为正的复数.通过对刚度矩阵B作如下修正来实现,即
式中,zj为正实数特征值λj对应的特征向量;Zj为复数特征值λj对应特征向量zj的实部和虚部组成的矩阵,即Zj=[Re(zj)Im(zj)];εj为选择参数.
文献[14-15]中证明,只要选取
即可保证在不改变系统的模态阶数及正常模态的特征值的前提下,将含有正实部的特征值移动到复平面的负半轴.本文算法是在Matlab R2009b软件包中完成的.
3 算例分析
3.1 算例1:半无限楔形体平面内波动
考察如图1所示的二维楔形体,模型尺寸为x1=b,,x2=b,y2=0,其中b=1 m.介质的力学参数为:剪切模量为G=1 Pa,质量密度为ρ=1 kg/m3,泊松比为ν=0.25.假定在关键点1处施加集中力P(t),该载荷时程呈三角形变化,如式(11)所示
按平面应变问题进行分析.采用1个8结点的高阶单元[16]进行离散,共8个结点16个自由度,相似中心位于坐标原点O处,SBFE(scaled boundary finite element)网格如图2所示.采用比例边界有限元法分析时各边界自由.
为了说明高阶透射边界的精确性,针对该问题采用有限元扩展的网格进行了分析,有限元网格采用8结点四边形单元离散,共划分了2 000个单元6 241个结点,限于篇幅本文未给出网格图.有限元分析的边界条件为最外层的竖直边界固定,其余边界自由.时域计算总时间为20 s,积分步长为∆t=0.02 s.
将稳定性措施前后的结构响应进行了对比分析,以观测点3的水平位移为例,如图3所示.采取稳定性措施前其位移响应是发散的,而采用了移谱法后消除了系统的虚假模态,其位移响应是收敛的.可以看出,本文的稳定性方法是非常有效的.
稳定性措施后的观测点1和3的水平向位移时程如图4所示.可以看出,高阶边界(Mcf=15)的结果精度略差,若逐渐增大连分式展开的阶数,高阶透射边界(Mcf=20)的结果与扩展边界的结果吻合得很好.
3.2 算例2:半圆形河谷平面内波动
考虑如图5所示的半圆形河谷,其半径b=1 m.介质的力学参数为:剪切模量为G=1 Pa,质量密度为ρ=1 kg/m3,泊松比为ν=0.25.假定在A点作用竖直向下的力P(t),其随时间变化见式(11).
按平面应变问题进行分析.采用4个8结点的高阶单元进行离散,共29个结点58个自由度,相似中心位于坐标原点O处,SBFE网格如图6所示;比例边界有限元分析时各边界自由.
同理,针对该问题采用有限元扩展的网格进行了分析,考虑到对称性,有限元模型选取半圆形河谷的1/2.有限元网格采用8结点四边形单元离散,共划分了3 000个单元9 221个结点;有限元分析的边界条件为最外层的弧形边界固定,另一垂直边界考虑为对称约束.时域计算总时间为20 s,积分步长为∆t=0.02 s.
观测点A和B的竖向位移时程如图7所示.可以看出,高阶透射边界(Mcf=12)的结果与有限元扩展边界的结果基本一致,若继续增大连分式展开的阶数,高阶边界的结果将更精确.
4 结语
不定积分中分部积分法的新探究 篇4
关键词：不定积分,分部积分公式,微分
高等数学是所有高等学校的一门必修课, 而微积分是高等数学中的重要内容。不定积分研究的是与微分运算正好相反的问题:求一个可导函数, 使它的导函数等于已知的函数。不定积分是微分运算的逆运算, 是微分学和积分学的联系纽带。求不定积分的常规方法包括直接积分法、换元积分法和分部积分法, 其中的分部积分法是教学过程中的一个难点。
分部积分法是利用微分公式d (uv) =udv+vdu, 推导得到了分部积分公式[1].因此, 当不容易直接积出, 而较为容易求出时, 可以采用分部积分公式作为转换。一部分教材[1,2]指出分部积分法的关键是要正确地选择u和v, 选择时应兼顾如下两点: (1) dv容易求出; (2) 容易积出。但是, 这样的描述比较笼统。还有一些教材[3]列举了一些常用的选择方法, 如当被积函数为幂函数和指数函数相乘时, 选择幂函数为u;当被积函数为幂函数和对数函数相乘时, 选择对数函数为u;当被积函数为幂函数和反三角函数相乘时, 选择反三角函数为u;文[4]将基本初等函数分为“低、中、高”三类, 用“低等服从高等”的思路来解决u、v的选择问题;文[5]、[6]则总结出一套“口诀”, 这样都需要学生死记硬背, 不易于学生的理解和掌握。
本文将以例题的形式进一步分析分部积分公式, 不再局限于函数u、v的定义和选择, 而是将被积函数进行分类, 根据被积函数的特点来进行u、v的设定, 进而让学生能够灵活使用分部积分法进行不定积分的计算。
一、类型一:可降幂型
解:被积函数为两个函数的乘积, 利用直接积分法和换元积分法无法求出原函数, 可以考虑使用分部积分法。根据公式, 选择其中一个函数为u, 则剩下的部分就是dv.
若设u=x, 则, 应用分部积分公式 (1) 得
此类型的被积函数为幂函数与指数函数或三角函数的乘积。此时, 一般将幂函数设为u, 指数函数或三角函数则为dv.这是因为应用分部积分公式之后, 幂函数通过微分后次数降低一次, 使得转换后的不定积分较容易求积出。
二、类型二:直接型
解:被积函数为单个函数, 可以将被积表达式lnxdx直接看作udv, 即u=lnx、v=x, 于是
解:被积函数为指数函数与反三角函数的乘积。由于arctanx的原函数不能直接得到, 于是设u=arctanx、, 可得
此类型的被积函数为分对数函数或反三角函数 (即不能直接积分的函数) 与其他函数 (即能直接积分的函数) 的乘积。此时, 将分对数函数或反三角函数设为u, 其他函数则为dv.
三、类型三:循环型
解:被积函数为指数函数与三角函数的乘积, u、v可以任意设定。若设u=sinx、, 则
若设, 而sinxdx=-d (cosx) , 则
此类型的被积函数为指数函数与正弦 (或余弦) 函数的乘积。此时, 需要使用分部积分公式两次才能找到原函数。任意设定其中一部分函数为u, 其他函数则为dv.分部积分两次后会还原到原来的函数, 只是系数有一些相应的变化。因此, 等式两边就含有系数不同的同一积分。
一元函数的不定积分是微积分中的重要知识点, 对定积分的学习有非常重要的作用。教师通过对函数的类型、性质等的细致观察、理解和分析, 可以培养学生的发现问题、分析问题和解决问题的能力, 进而有助对于学生创新能力的培养。本文针对不定积分中的分部积分法, 通过对被积函数类型的分析将其分为三类:降幂型、直接型和循环型, 简化计算过程, 帮助学生合理、有效地使用分部积分公式。
参考文献
[1]祁爱琴, 邵珠艳, 胡西厚.医用高等数学[M].北京:科学出版社, 2013:77-97.
[2]王培承, 祁爱琴, 魏曼莎.医科高等数学[M].济南:山东人民出版社, 2010:76-98.
[3]同济大学应用数学系.高等数学[M].第六版.北京:高等教育出版社, 2007:184-222.
[4]范梅.不定积分的分部积分法探究[J].西安航空学院学报, 2015, 33 (1) :66-68.
[5]胡结梅, 郑华盛.不定积分计算方法注记[J].高等数学研究, 2014, 17 (6) :10-13.