半空间问题

半空间问题（共4篇）

半空间问题 篇1

0 引言

压电材料[1]是一种能够实现电能和机械能相互转化的智能材料, 已被广泛用来制造压电传感器和驱动器。优越的力电转化性能, 使其在智能结构和微电机系统中发挥着重要的作用[2]。接触问题是研究固体间相互接触时在接触区的作用状况, 包括位移、应变和应力, 以及由此引起的强度分析方面的问题。接触问题体在研究接触表面的破坏机理, 通过压痕实验测定材料力学性能等方面具有重要的意义。Johnson[3]则在其经典著作中给出了弹性体材料的接触问题分析方法。最近, Giannakopoulos等[4,5]用Hankel积分变换方法得到了压电材料半空间在导电体和绝缘体压头作用下的三维接触问题的解析解, 同时, 该文还针对三种商业压电材料进行了有限元分析, 并对理论解与Abaqus有限元解进行比较。该文主要侧重理论分析, 未给出有限元分析的详细过程, 以及应力分布, 压痕- 压力关系以及压力和电势的关系曲线。现在, 有限元软件的使用越来越广泛, 如王光建[6]基于ANSYS软件分析了双螺旋副斜齿轮强度, 为材料力学性能问题的解决提供了新的思路。有限元在有限元分析接触问题的方面, 孙林松等[7]列举了三种数值方法来分析接触问题;李一耕[8]通过ANSYS进行了推力球轴承赫兹接触问题的计算, 所得结果与赫兹理论解进行对比;王振波[9]通过有限元方法分析了MEMS涂层的接触问题, 为MEMS涂层的可靠性设计提供了理论基础。赵亚素[10]基于ANSYS软件分析了非均匀涂层在均布荷载作用下的接触问题

本文利用大型有限元分析软件ANSYS, 分析压电半空间的三维无摩擦接触问题。通过对三种目前可供商业使用的压 · 电材料进行分析和建模, 求解了在特定荷载值作用下, 三种压电材料在圆柱形和球形压头作用下的压入深度、表面接触应力和电势。

1 问题的描述

三维压电半空间接触问题的几何模型如图1所示, 本文主要研究平底圆柱 (图1a) 和球形 (图1b) 两种几何形状刚性压头作用于压电半空间的接触问题。柱坐标体系如图所示, 平面半径方向由r表示, 外法线方向由z表示。压电半空间的尺寸远大于压头尺寸。接触面为z=0 的平面, 由作用力P引起的压入深度为 δ0, 同时, 接触半径为a。压电材料受压的极化轴与z轴一致。

2 有限元模型建立及求解

本文针对三种可供商业使用的压电材料PZT-4、Ba Ti O3和PZT-5A进行分析。三种压电材料半空间在刚性压头作用下的接触问题, 三种压电材料性能参数参见表1。

2.1 刚性平底圆柱压头算例

刚性圆柱压头半径为0.2m, 压电半空间半径为2m, 厚度为2m。荷载P分别取0.2N, 2N, 5N, 10N, 15N, 20N, 25N, 30N, 35N, 40N, 50N共12 组模拟。NSYS建模过程中刚性压头采用SOLID185 单元, 压电半空间采用SOLID5 耦合场单元。如图2 所示, 因为是轴对称模型, 建模时将模型简化为1/4 模型, 提高运算效率。网格划分时采用六面体网格, 并在接触区域进行网格细化, 提升计算精度。接触对设置将压头下表面设置为目标面, 压电半空间上表面设置为接触面, 接触表现设置为Standard, 接触方向设置为On nodes-Normal from contact。PZT-5A压电半空间在P=2N时平底圆柱压头作用下的压入深度、表面接触应力与电势的云图分别见图3、图4 与图5。

2.2 刚性球形压头算例

刚性球形压头半径R=0.2m, 压电半空间半径为2m。施加的荷载P与刚性圆柱形压头相同。球形压头采用solid187 单元, 压电半空间采用耦合场solid98 单元, 建立1/4 模型。如图6 所示, 因为刚性球形压头与压电半空间接触区域较小, 需要在接触区域附近的网格划分足够密, 因此采用四面体网格。接触对设置将压头下表面设置为目标面, 压电半空间上表面设置为接触面, 接触表现设置为Standard, 接触方向设置为On nodes-Normal from contact。PZT-5A压电半空间在P=0.2N时平底圆柱压头作用下的压入深度、表面接触应力与电势的云图分别见图7、图8 与图9。

3 数值结果与分析

以下将对圆柱形压头和球形压头作用下的接触应力分布, 压力- 压痕曲线以及压力- 电势曲线进行分析。

首先, 为了验证结果的正确性, 本文计算了球形压头最大压入深度、接触半径。计算结过如表2 所示, 本文的结果与文献[13] 符合较好。

图10 和图11 给出了三种不同压电材料在圆柱形压头作用下的压痕和接触应力分布。从图10 中可以发现, 同样作用力作用下, 平底圆柱压头压入深度从大到小依次为PZT-5A材料、PZT-4 材料、BaTi O3材料。从图11 中可以发现, 相同作用力作用下, 三种材料的表面接触应力数值接近, 在接触区域内变化较为平缓, 接触边界处同时存在大幅度的急剧增高, 这一现象与解析解结论中接触边界存在的奇异性相吻合。

图12 和图13 给出了三种不同压电材料在平压头作用时, 压入深度 δ0与压力P以及最大电势 ϕ0与压力P的关系曲线。从图12 中可以发现, δ0与P呈线性关系。相同荷载作用下, PZT-5A材料产生的δ0最大, 而Ba Ti O3 产生的 δ0最小。从图12 中可以看到, ϕ0与P同样呈现出线性关系。在相同荷载作用下, PZT-5A材料产生较大的电势而Ba Ti O3产生较小的电势。

如图14 和图15 所示, 球形压头作用下的最大压入深度 δ0与最大接触应力P0被给出。从图中可以发现, 相同作用力作用下, 产生 δ0由大到小的压电材料依次为PZT-5A、PZT-4、Ba Ti O3。而产生的最大应力则刚好相反。

图16 和17 分别给出了当荷载P=0.2N时, 压电半空间表面的应力分布和电势分布。从图中可以发现, 最大接触应力P0出现在接触区域的中心。其中, Ba Ti O3材料产生的最大接触应力最大而接触区半径最小。PZT-5A材料产生的最大接触应力最小而接触区域半径最大。如图17 所示, 同样条件下, 最大电势 ϕ0同样出现在接触区域中心, 电势沿径向逐渐降低;电势由大到小依次为PZT-5A、PZT-4、Ba Ti O3, 整体趋势与理论结果相同。

图18 给出了最大电势 ϕ0与作用力P的关系曲线。从图中可以看出, 在相同荷载作用下, BaTi O3材料产生的电势最小, 而PZT-5A产生的结果最大。

4 结论

利用ANSYS软件模拟接触问题可以得到比较精确的结果。

相同荷载作用下, Ba Ti O3压电材料产生较大的应力分布和较小的最大压痕, 而PZT-5A压电材料产生较小的应力分布和较大的压痕。

在相同荷载作用下, PZT-5A材料产生较大的电势, PZT-4 材料产生电势次之, 而BaTiO3材料产生的电势较小。

半空间问题 篇2

边界元法因其具有边界离散简单、计算结果精度高且自动满足远场辐射条件等优点而被广泛应用于声学问题的求解中[1,2]。在求解半空间声学问题时,可依据镜像原理将半空间问题转化为全空间问题的方式来求解,但是该方法需要对原模型与镜像模型进行网格离散,如此就增大了半空间问题求解的计算量及存储量。考虑到边界元法本身求解声学问题时,所需计算量为O(n2)量级(n为自由度数),那么,随着求解问题规模的增大,全空间镜像方法的计算效率明显降低。

20世纪80年代,Rokhlin[3,4]提出了快速多极算法,并推导了求解二维和三维声学问题的快速多极算法理论。快速多极算法解决了传统边界元法计算量大的瓶颈问题,将计算量由O(n2)降低到O(n)量级,极大地提高了计算效率。快速多极算法也迅速成为声学领域的研究热点方法。此后,国内外学者在快速多极边界元方法改进与完善方面做了大量研究工作[5,6]。

本文在伽辽金多极边界元法的基础上,引入半空间格林公式[7],求解了二维半空间声学问题。在数值计算过程中,由半空间格林函数引入的镜像点与原模型对应点采用相同处理方式,避免了全空间方法中由于求解镜像模型带来的计算量与存储量过大的问题。在求解系统线性方程组时,使用近似求逆预处理技术处理系数矩阵,改善了矩阵性态,提高了系统方程的求解效率。二维半空间点声源辐射算例验证了伽辽金多极边界元法的计算精度,且与全空间镜像方法计算时间的对比结果表明,半空间伽辽金多极边界元法具有更高的计算效率。最后使用半空间方法求解了T形声屏障的声散射问题。

1 二维半空间声学伽辽金边界元法

在各向同性均匀介质中,稳态线性声学的边界积分方程可表示为

$\begin{array}{l}c(x)p(x)={\displaystyle {\int }_{\Gamma }\frac{\partial G(x,y)}{\partial \mathit{n}(y)}}p(y)\text{d}\Gamma (y)-\\ {\displaystyle {\int }_{\Gamma }G}(x,y)q(y)\text{d}\Gamma (y)\forall x\in \Gamma (1)\end{array}$

式中,c(x)为与点x处的边界几何特征相关的系数,边界光滑时c(x)取值为0.5;p(x)为x处的声压;n(y)为y点的外法向导数;Γ为求解问题的边界;G(x,y)为格林函数。

为便于书写,将式(1)中右端两个积分项分别记为(Kp)(x)和(Vq)(x)。

在二维声学问题中,自由空间格林函数可表示为

$G(x,y)=\frac{\text{i}}{4}{\rm H}{}_0^{(1)}(kr)$ (2)

k=ω/c

式中,r为源点x和场点y之间的欧氏距离;H(1)0()为第一类汉克尔函数;k为波数;ω为角频率;c为声音在均匀介质中的传播速度。

在求解半空间声学问题时,半空间格林函数的表达式为

$G(x,x^{\prime },y)=\frac{\text{i}}{4}[{\rm H}{}_0^{(1)}(k|\overline{xy}|)+$

$R{}_{\text{p}}{\rm H}{}_0^{(1)}(k|\overline{x^{\prime }y}|)]$ (3)

式中,x′为x的镜像点;Rp为对称面反射系数,Rp=1时,即表示对称面为刚性面。

半空间声学问题的边界条件可表示为

式中,ΓD为狄利克雷边界;ΓN为纽曼边界。

二维半空间声学问题示意图如图1所示。图1中E表示外部声场,E′表示内部声场,Sp表示对称面。x′表示x关于对称面Sp的镜像点。

对式(1)求偏导后,得到含有超奇异积分的边界积分方程为

$\begin{array}{l}c(x)q(x)={\displaystyle {\int }_{\Gamma }\frac{\partial G(x,y)}{\partial \mathit{n}(y)\mathit{n}(x)}}p(y)\text{d}\Gamma (y)-\\ {\displaystyle {\int }_{\Gamma }\frac{\partial G(x,y)}{\mathit{n}(x)}}q(y)\text{d}\Gamma (y)\forall x\in \Gamma (5)\end{array}$

式(5)中右边两个积分项依此记为-(Dp)(x)和(K′q)(x)。

在求解外场声学问题时,边界积分方程在特征频率处会出现解不唯一的问题,该问题属于纯数学问题,不具有任何物理意义。本文使用Burton-Miller方程[8]来保证在所有频率上解的唯一性,即使用式(1)和式(5)的线性组合求解外场声学问题,Burton-Miller方程可写为

$\begin{array}{l}\frac{\text{i}}{k}{\displaystyle {\int }_{\Gamma }\frac{\partial G(x,y)}{\partial \mathit{n}(y)\mathit{n}(x)}}p(y)\text{d}\Gamma (y)+{\displaystyle {\int }_{\Gamma }\frac{\partial G(x,y)}{\partial \mathit{n}(y)}}p(y)\text{d}\Gamma (y)+\\ c(x)p(x)=\frac{\text{i}}{k}[{\displaystyle {\int }_{\Gamma }\frac{\partial G(x,y)}{\partial \mathit{n}(x)}}q(y)\text{d}\Gamma (y)-c(x)q(x)]+\\ {\displaystyle {\int }_{\Gamma }G}(x,y)q(y)\text{d}\Gamma (y)(6)\end{array}$

式中,i/k为式(1)和式(5)线性组合中的耦合系数。

边界单元节点声压和法向速度分别使用插值函数φp(x)和φq(x)获得,即

式中,ps和qt分别为节点声压值及法向速度值。

在式(6)中引入线性测试函数,并使用式(1)和式(5)中积分的算符表示形式,可得到Burton-Miller方程的矩阵方程表达形式:

$(-\frac{1}{2}\mathit{{\rm I}}+\mathit{{\rm K}}+\frac{\text{i}}{k}\mathit{D})\mathit{p}=(\mathit{V}-\frac{\text{i}}{2k}\mathit{{\rm I}}-\frac{\text{i}}{k}{\mathit{{\rm K}}}^{\prime })\mathit{q}$(8)

在使用直接法(高斯消去法或LU分解法)求解系统线性方程组时,计算量达到O(n3)量级。即便是使用共轭梯度法(conjugate gradient squared method, CGS)或广义极小残差算法(generalized minimum RESidual,GMRES)等迭代法求解,计算量仍有O(n2)量级(n为系统矩阵的秩)。高计算量成为限制伽辽金边界元法求解大规模声学问题的主要原因。

2 半空间伽辽金多极算法

快速多极算法的主要思想是通过将边界积分方程中的基本解以级数形式展开,分离场点和源点之间的联系,实现加速矩阵与向量乘积运算的目的。本文将快速多极算法与伽辽金边界元法相结合,详细描述了伽辽金多极边界元法求解二维半空间声学问题的过程,包括核函数的多极展开、多极(局部) 展开系数的计算与转移以及伽辽金多极边界元法中矩阵向量乘积运算的实现过程。

2.1 核函数的多极展开

在使用半空间方法求解半无限域声学问题时,需在边界积分方程中引入半空间格林函数。基于汉克尔函数展开方法[9],当源点x、场点y以及镜像点x′满足$|\overline{xx{}_{\text{c}}}|<|\overline{yx{}_{\text{c}}}|$、$|\overline{yy{}_{\text{c}}}|<|\overline{xy{}_{\text{c}}}|$、$|\overline{x^{\prime }{x}^{\prime }{}_{\text{c}}}|<|\overline{y{x}^{\prime }{}_{\text{c}}}|$时,二维半空间格林函数的多极展开形式为

$G(x,x^{\prime },y)\approx \frac{\text{i}}{4}{\displaystyle \sum _{m=0}^{L-1}(}U{}_m(x,x{}_{\text{c}}){\rm T}{}_m(x{}_{\text{c}},y{}_{\text{c}})+$

Um(x′,x′c)Tm(x′c,yc))Vm(yc,y) (9)

$U{}_m(x,x{}_{\text{c}})=\frac{1}{\sqrt{L}}\text{e}{}^{\text{i}k\overline{xx{}_{\text{c}}}\cdot \widehat{\mathit{s}}}$ (10)

$\mathit{{\rm T}}{}_{\mathit{m}}(\mathit{x}{}_{\mathbf{c}},\mathit{y}{}_{\mathbf{c}})={\displaystyle \sum _{\mathit{l}=-\mathit{{\rm N}}{}_{\mathit{l}}}^{\mathit{{\rm N}}{}_{\mathit{l}}}\mathbf{i}}{}^{-\mathit{l}}\mathit{{\rm H}}{}_{\mathit{l}}^{(1)}(\mathit{k}|\overline{\mathit{x}{}_{\mathbf{c}}\mathit{y}{}_{\mathbf{c}}}|)\mathbf{e}{}^{\mathbf{i}\mathit{l}(\mathit{\theta }-\mathit{\beta }{}_{\mathit{m}})}$ (11)

$\mathit{V}{}_{\mathit{m}}(\mathit{y}{}_{\mathbf{c}},\mathit{y})=\frac{1}{\sqrt{\mathit{L}}}\mathbf{e}{}^{\mathbf{i}\mathit{k}\overline{\mathit{y}{}_{\mathbf{c}}\mathit{y}}\cdot \widehat{\mathit{s}}}$ (12)

$\widehat{\mathit{s}}=(\cos \beta {}_m,\sin \beta {}_m)$

βm=2πm/L m=0,1,…,L-1

式中,θ为向量$\overline{x{}_{\text{c}}y{}_{\text{c}}}$的极角;L和Nl分别为格林函数多极展开式的截断项数和计算汉克尔函数使用的积分点数。

镜像点与源点使用相同方式处理,同理可以得到含有镜像点x′的Um的表达式。

2.2 多极(局部)展开系数及转换关系

在伽辽金多极边界元法中,如图2所示,场点y (远离源点x)在展开点o3处的多极展开系数可以表示为

${\rm M}{}_m(o{}_3)={\displaystyle {\int }_{\Gamma {}_j}\text{e}}{}^{\text{i}k\overline{o{}_{3}y}\cdot \widehat{\mathit{s}}}q(y)\text{d}\Gamma (y)$ (13)

多极展开系数与源点x位置无关,当全部多极展开系数计算完成后,所有远离边界Γj形如式(1)中(Vq)(x)的积分都将通过多极展开系数快速计算出来,与镜像点x′相关的积分(Vq)(x′)采用相同方式处理。

在求解二维声学问题时,快速多极算法的实现过程包含对四叉树结构的上行遍历和下行遍历两个过程。在上行遍历过程中,需要计算各级栅格的多极展开系数,即从叶子栅格开始,多极展开系数从子栅格向父栅格传递。而在下行遍历过程中,从第二级栅格开始,局部展开系数从父栅格向子栅格传递。在下行遍历过程中,当所有局部展开系数计算完成时,对应源点积分通过局部展开系数计算。对所有源点循环一次,即完成一次系数矩阵与迭代矢量的乘法运算。

在对树结构进行上行遍历和下行遍历的过程中,多极展开系数间(M2M)与局部展开系数间(L2L)的传递公式如下:

$\mathit{{\rm M}}{}_{\mathit{m}}(\mathit{y}{}_{\mathbf{c}})=\mathbf{e}{}^{\mathbf{i}\mathit{k}\overline{\mathit{y}{}_{\mathbf{c}}\mathit{o}{}_3}\cdot \widehat{\mathit{s}}}{\rm M}{}_m(o{}_3)$ (14)

$\mathit{L}{}_{\mathit{m}}(\mathit{o}{}_1)=\mathbf{e}{}^{\mathbf{i}\mathit{k}\overline{\mathit{x}{}_{\mathbf{c}}\mathit{o}{}_1}\cdot \widehat{\mathit{s}}}L{}_m(x{}_{\text{c}})$ (15)

多极展开系数向局部展开系数转移(M2L)公式可表示为

Lm(xc)=TmMm(yc) (16)

在半空间声学问题中,快速多极边界元法中式(14)～式(16)三种传递方式详见图2。

在多极(局部)展开系数以及转移关系中,式(16)的求解计算量最大。为减小多极展开系数向局部展开系数转移的计算量,本文使用了自适应树结构[10],该自适应树结构可有效减少源点的相互作用列表中的单元数目,提高计算效率。此外,本文的截断项数L参照半经验公式选取[11],在保证算法计算精度的前提下,使用较小的截断项L,其中

L=kdl+plg(kdl+π) (17)

式中,dl为树结构第l层的栅格尺寸;p为与精度相关的常数。

由式(17)得到的展开系数截断项数会随栅格在树结构中所处层数的变化而变化。因此,在上下行遍历过程中,多极展开系数与局部展开系数在树结构不同层之间传递时,需使用有效的插值和滤波方法来保证传递的精度。本文使用基于快速傅里叶变换的插值方法[12]来保证不同层之间系数转移的计算精度。

2.3 矩阵向量乘积的实现

在使用伽辽金多极边界元法求解半空间声学问题时,以单元sj上节点τ为例,给出了式(8)中矩阵向量乘积(Dp)τ与(K′q)τ的求解过程。求解过程以与单元sj的距离关系分为两个部分:与sj相邻的栅格单元(记为Γ∈N)及与sj满足距离准则的栅格单元(记为Γ∈F),其中,N和F均为栅格集合。矩阵向量的乘积运算中相邻栅格以及满足距离准则的栅格单元对单元的作用分别通过直接求解以及伽辽金多极边界元法求解。实现过程如下:

$\begin{array}{l}({\mathit{{\rm K}}}^{\prime }\mathit{q}){}_{\tau }\approx {\displaystyle \sum _{\Gamma \in \mathit{{\rm N}}}(}{\mathit{{\rm K}}}^{\prime }\mathit{q}){}_{\tau }+{\displaystyle \sum _{\Gamma \in \mathit{F}}q}{}_{\Gamma }{\displaystyle \sum _{m=0}^{L-1}\frac{1}{4L}}U{}_m{\rm T}{}_m\cdot \\ k\widehat{\mathit{n}}(x){\displaystyle \sum _{i=1}^{G{}_n}\Delta }{}_{\Gamma }\omega {}_{\Gamma ,i}V{}_m(y{}_{\text{c}}^{\Gamma },y{}_{\Gamma ,i})\phi {}_{\Gamma ,\tau }^q(y{}_{\Gamma ,i})(18)\end{array}$

$\begin{array}{l}(\mathit{D}\mathit{p}){}_{\tau }\approx {\displaystyle \sum _{\Gamma \in {\rm N}}(}\mathit{D}\mathit{p}){}_{\tau }+{\displaystyle \sum _{\Gamma \in \mathit{F}}p}{}_{\Gamma }{\displaystyle \sum _{m=0}^{L-1}\frac{\text{i}}{4L}}U{}_m{\rm T}{}_m\cdot \\ k{}^2\widehat{\mathit{n}}(x)\cdot \widehat{\mathit{n}}(y){\displaystyle \sum _{i=1}^{G{}_n}\Delta }{}_{\Gamma }\omega {}_{\Gamma ,i}V{}_m(y{}_{\text{c}}^{\Gamma },y{}_{\Gamma ,i})\phi {}_{\Gamma ,\tau }^p(y{}_{\Gamma ,i})(19)\end{array}$

$\widehat{\mathit{n}}(x)=\mathit{n}(x)\cdot \widehat{\mathit{s}}\widehat{\mathit{n}}(y)=\mathit{n}(y)\cdot \widehat{\mathit{s}}$

式中,ΔΓ、ωΓ,i与yΓ,i分别为和单元Γ相关的雅可比行列式、加权系数以及高斯积分点;φpΓ,τ()为单元Γ上节点τ的插值函数;yΓc为单元Γ所在栅格的中心点。

基于式(9)中半空间格林函数的多极展开式,式(8)中其余矩阵向量乘积均可采用类似方式得到。对于Burton-Miller方程中含有的超奇异与奇异积分分别使用正则化方法以及奇异值消去法计算。

3 预处理方法

分离式(8)中的已知量和未知量,可得到形如Aλ=b的线性方程组表达式。本文采用广义极小残差算法(generalized minimal RESidual, GMRES)求解系统线性方程组,并使用近似求逆(approximate inverse, AI)预处理方法对系数矩阵进行预处理。近似求逆预处理方法主要思想即构造稀疏近似矩阵P-1,使得P-1≈A-1,并将P-1用作预处理矩阵。在构造矩阵P-1时,首先使用共节点相邻单元的方式确定矩阵P-1的稀疏模式,然后通过求解Frobenius范数的极小化问题来构造近似矩阵P-1。近似求逆预处理对迭代法求解效率的影响在方法数值验证时给出。

4 数值算例

4.1 数值验证

半空间点声源传播模型为:点声源S位于原点(0,0)处,声压幅值AS=1,声音接收点RS位于坐标(4,0)处, 单位为m。声速c=343m/s,ρ=1.225kg/m3,Sp表示刚性对称面。在验证数值方法精度时,引入单位圆边界Γ,其圆心位于原点(0,0)处,与点声源S的位置重合,如图3所示。

在图3所示半空间问题中,不考虑单位圆边界Γ,由点声源S引起的接收点RS处声压可表示为

$\text{ϕ}(R{}_S)=A{}_S\frac{\text{e}{}^{\text{i}\text{k}r}}{4\text{π}r}+A{}_{S^{\prime }}\frac{\text{e}{}^{\text{i}\text{k}{r}^{\prime }}}{4\text{π}{r}^{\prime }}$(20)

式中,声源的幅值AS=AS′=1N/m;r为声源S到接收点RS的距离;r′表示镜像声源S′到接收点RS的距离。

为验证本文半空间方法的计算精度,引入单位圆边界Γ,将半空间点声源的传播过程划分成两部分计算:首先,计算由点声源S辐射到单位圆边界节点上的声压值;然后,将得到的声压值作为边界条件,计算半空间圆形边界Γ的声辐射问题(第一部分使用点声源公式计算即可,第二部分使用半空间伽辽金多极边界元法计算);最后,将使用半空间伽辽金多极边界元法得到的接收点RS处声压与式(20)得到的解析解作对比。

本文将半空间伽辽金快速多极边界元计算结果记为pG,理论解记为pΓ,定义接收点RS处的声压相对误差e为

$e=\frac{|p{}^G-p{}^{\Gamma }|}{|p{}^{\Gamma }|}$(21)

图4给出了在1000Hz时,使用半空间伽辽金多极边界元法得到的接收点RS处声压相对误差变化曲线。当单位圆边界划分的自由度数增加到1084时,误差减小到0.004,随着自由度数的增加,伽辽金多极边界元法的计算误差持续减小,验证了伽辽金快速多极边界元法的计算精度。

图5给出了当单位圆边界Γ划分为9726个自由度时,使用和不使用预处理的迭代求解器收敛速度对比。未经过预处理的GMRES求解器迭代40步,经过近似求逆预处理后,只需18次迭代即可达到收敛精度。结果表明,近似求逆预处理可以有效改善矩阵的性态,大幅减少迭代次数,加速迭代求解器的收敛。

表1给出了基于伽辽金多极边界元法的半空间和全空间镜像方法求解效率对比。从表1中可以看出,在求解半空间问题时,由于镜像模型的存在,全空间镜像方法需求解原模型两倍的自由度数。随着半空间模型自由度增加,两种方法之间的计算时间差距不断增大,半空间方法的计算效率明显高于全空间镜像方法的计算效率。

4.2 T形声屏障声散射分析

使用半空间伽辽金多极边界元法求解二维声屏障散射问题。二维T形声屏障散射模型如图6所示。声屏障主要尺寸为:高度hB=2m,厚度t=0.1m,T形顶部尺寸wT=0.5m,声屏障底边中点为原点(0,0),点声源BS位于(2,0)m,声压幅值pBS=1Pa。声音在空气介质的传播速度c=343m/s,空气介质密度ρ=1.29kg/m3。

声屏障边界划分成2450个离散单元,以及3200个场点单元,图7给出了k=2时,半空间声屏障声散射问题的声压级等值线图。该算例说明伽辽金多极边界元法可有效地求解工程中的半空间声学问题。

5 结论

(1)使用经过近似求逆预处理的广义极小残差迭代求解器求解系统线性方程组时,可明显改善方程组系数矩阵的性态,减少求解器的迭代步数,提高计算效率。

(2)使用半空间伽辽金多极边界元求解半空间问题时,可略去全空间镜像方法中由于引入镜像模型而带来的计算量。随着求解模型规模增大,半空间镜像方法比全空间镜像方法更加高效。

参考文献

[1] Chen Z S, Hofstetter G, Mang H A. A Symmetric Galerkin Formulation of the Boundary Element Method for Acoustic Radiation and Scattering[J]. J. Comput. Acoust., 1997, 5(2): 219-241.

[2] Schenck H A. Improved Integral Formulation for Acoustic Radiation Problems[J]. J. Acoust. Soc. Am., 1968, 44(1):41-48.

[3]Rokhlin V.Rapid Solution of Integral Equations ofScattering Theory in Two Dimensional[J].J.Com-put.Phys.,1990,86:414-439.

[4] Rokhlin V. Diagonal Forms of Translation Operators for the Helmholtz Equation in Three Dimensions[J]. Appl. Comput. Harm. Anal., 1993,1(1): 82-93.

[5] Gumerov N A, Duraiswami R A. Broadband Fast Multipole Accelerated Boundary Element Method for the Three Dimensional Helmholtz Equation[J]. J. Acoust. Soc. Am., 2009,125(1):191-205.

[6]Wu H J,Liu Y J,Jiang W K.Analytical Integra-tion of the Moments in the Diagonal Form Fast Mul-tipole Boundary Element Method for 3-D AcousticWave Problems[J].Eng.Anal.Bound.Elem.,2012,36:248-254.

[7] Seybert A F, Soenarko B. Radiation and Scattering of Acoustic Waves From Bodies of Arbitrary Shape in a Three-dimension Half Space[J]. ASME Trans. J. Vib. Acoust. Stress Reliab. Des., 1988,110:112-117.

[8]Burton A J,Miller G F.The Application of the In-tegral Equation Methods to the Numerical Solutionof Some Exterior Boundary-value Problems[J].Proc.Roy.Soc.London,1971,A323:201-210.

[9] Labreuche C. A Convergence Theorem for the Fast Multipole Method for 2 Dimensional Scattering Problems[J]. Math. Comput., 1998, 67:553-591.

[10] Bapat M S, Liu Y J. A New Adaptive Algorithm for the Fast Multipole Boundary Elements Method[J]. CMES-Comp. Model Eng., 2010,58:161-183.

[11]Coifman R,Rokhlin V,Wandzura S.The Fast Multipole Method for the Wave Equation:a Pe-destrian Prescription[J].IEEE Antennas Propa-gation Mag.,1993,35(3):7-12.

半空间问题 篇3

1 材料与方法

1.1 实验动物与试剂

清洁级健康雄性ICR小鼠,体重18～22 g,购自吉林大学基础医学院动物实验中心,动物许可证号SCXK-(吉)-2007-0003,自由吃食和饮水,室温(20±2)℃。D-半乳糖购自上海蓝季科技发展有限公司(产品批号:070313);羧甲基纤维素钠购自美国Sigma公司;姜黄素购自国药集团化学试剂有限公司(产品批号:F20050802);AChe、ChAT考马斯亮蓝蛋白测试试剂盒均购自南京建成生物工程有限公司。

1.2 仪器与设备

水迷宫视频跟踪分析仪和BI 2000图像分析系统为成都泰盟科技有限公司产品,722分光光度计为上海奥普勒仪器有限公司产品,离心机为北京医用离心机厂产品,组织匀浆机为江苏金坛金城国胜实验仪器厂产品。

1.3 分组与给药

将50只小鼠按体重随机分为5组,即正常对照组、模型组、预处理组、同时处理组、后处理组,每组10只。对照组:上午颈部皮下注射等量生理盐水,下午给予1%羧甲基纤维素钠溶液灌胃1次;模型组:上午颈部皮下注射120 mg/kg D-半乳糖,下午给予1%羧甲基纤维素钠溶液灌胃1次;预处理组:上午姜黄素用1%羧甲基纤维素钠配成混悬液,按200 mg/(kg·d)灌胃给药,下午颈部皮下注射120 mg/kg D-半乳糖;同时处理组:上午颈部皮下注射120 mg/kg D-半乳糖,同时姜黄素用1%羧甲基纤维素钠配成混悬液,按200 mg/(kg·d)灌胃给药;后处理组:上午颈部皮下注射120 mg/kg D-半乳糖,下午姜黄素用1%羧甲基纤维素钠配成混悬液,按200 mg/(kg·d)灌胃给药;连续45 d。

1.4 Morris水迷宫行为测试[3]

水迷宫由圆形水池、安全岛和记录系统3部分组成。水池直径130 cm,高60 cm,水池内壁涂成乳白色。池中水深30 cm,水温保持在(24±2)℃,用奶粉使池水浑浊。池壁标记4个入水点,将水池分为4个均等的象限,池壁上固定位置贴有一幅画作为参照物。安全岛置于第1象限中点,直径10 cm,高28 cm,距水面2 cm。记录系统为成都泰盟公司提供的BI 2000系统的行为学检测模块。从第39天开始进行定位航行实验,连续测试5 d。将小鼠从安全岛之外的3个象限任一入水点面向池壁放入置有安全岛的水池中,经 BI 2000图像处理系统自动记录小鼠在120 s内找到安全岛的时间(逃避潜伏期)和游泳路径(搜索距离),并让小鼠在岛上停留10 s。超过120 s找不到安全岛的记为 120 s,并引导小鼠上安全岛。定位航行实验结束后,间隔1天,拆除安全岛,将小鼠从同一入水点放入池中,记录小鼠在120 s内寻找安全岛过程中在原安全岛所在象限内的搜索时间(空间探索时间)及探索距离。

1.5 脑匀浆AChE和ChAT活力的测定

行为学测试完成后,断头处死小鼠,在冰上迅速分离全脑,称重后用相当于脑组织质量9倍的生理盐水将全脑组织制成匀浆,于4 ℃以3 500 r/min(离心半径为16 cm)离心10 min,去上清液,用比色法测定脑组织中AChE、ChAT的活力,方法按照说明书进行。

1.6 统计分析

数据采用SPSS 10.0统计软件进行单因素方差分析。组间比较采用最小显著差异法检验(LSD法)。

2 结果

2.1 姜黄素对AD模型小鼠空间探索实验的影响

在Morris水迷宫实验中,模型组与对照组比较,小鼠空间探索时间明显延长,探索距离增大(P<0.05),学习 记忆能力明显减退;模型组与各处理组比较,探索时间延长,探索距离增加(P<0.05);预先处理组和同时处理组探索时间短,探索距离小,与后处理组比较,差异有统计学意义(P<0.05),结果见表1。

注:与模型组比较,aP<0.05;与后处理组比较,bP<0.05。

2.2 姜黄素对AD模型小鼠脑组织AChE、ChAT活力的影响

模型组小鼠脑组织中AChE活性升高,ChAT活力降低,与正常对照组比较,差异有统计学意义(P<0.05)。灌胃给予姜黄素后小鼠脑组织中AChE活力降低,ChAT活力升高,与模型组比较,差异有统计学意义(P<0.05);预先处理组和同时处理组脑组织中AChe活力比后处理组有显著降低(P<0.05),预先处理组ChAT活力与同时处理组和后处理组比较有显著增高,差异均有统计学意义(P<0.05)。

注:AChE—乙酰胆碱酯酶;ChAT—胆碱乙酰转移酶。与模型组比较,aP<0.05;与后处理组比较,bP<0.05。

3 讨论

目前,认为中枢胆碱能系统与学习记忆功能的关系极为密切,有资料表明,海马的胆碱能系统直接参与信息的储存和回忆。胆碱能神经元分泌的神经递质乙酰胆碱(ACh)是中枢系统中与认知功能最直接相关的递质,一旦缺乏即会出现记忆力下降,注意力不集中等症状[4]。ChAT和乙酰辅酶A结合后,活性中心的阴离子部位即与胆碱结合,将乙酰基转移到胆碱上,合成ACh,因而ChAT是ACh合成的限速酶,是胆碱能活力的重要标记,两者共同维持着中枢胆碱能系统重要的神经递质——ACh的动态平衡,在多种高级脑功能包括学习记忆过程中发挥重要作用,所以用AChE及ChAT的活力来反映胆碱能系统的功能较为可靠[5]。大量研究表明,AD患者ACh合成减少,ChAT的活力下降,且与痴呆的严重程度密切相关[6],海马和皮层锥体细胞的突触及突触烟碱受体、突触前M型受体亦减少。

本实验用Morris水迷宫观察了姜黄素对D-半乳糖所致小鼠学习记忆的影响,结果表明不同给药方式的姜黄素均能提高小鼠学习记忆成绩,小鼠脑组织中ChAT活力升高,AChE活力降低,但预处理组和同时处理组对痴呆小鼠学习记忆影响与后处理组差异有统计学意义(P<0.05),提示姜黄素改善中枢胆碱能系统功能是其增强模型小鼠学习记忆能力的基础。姜黄素对自由基系统的调节作用和对胆碱能神经系统的保护作用[7],可能是改善学习记忆障碍的主要机制。本研究表明,姜黄素对AD所致学习记忆障碍有明显的改善作用,其作用机制可能与抑制AChE作用密切相关。

参考文献

[1]Youssef KM,EI-Sherbeny MA.Synthesis and antitumor activityof somecurcumin analoga[J].Arch Pharm,2005,338:181-189.

[2]Ng TP,Chiam PC,Lee T,et al.Curry concumption and cognitivefunc-tion in the elderly[J].Am J Epidemiol,2006,164:898-906.

[3]陈奇.中药药理研究方法学[M].北京:人民卫生出版社,2004:33.

[4]Wu P,Tarasenko Y,Gu Y,et al.Region-specific genernation of cholin-ergic neurons from fetal human neural stem cells grafcedin adult rat[J].Nat Neurseci,2002,5:1271-1276.

[5]Oda Y.Choline ace tyltransferase:the structure,distribution and patho-logic changes in the central nervous system[J].Pathol Int,1999,49:921-929.

[6]Lopez OL,Hamilton RL,Becker JT,et al.Severity of cognitive impair-ment and the clinical diagnosis of AD with lewy[J].Bodies Neurol,2000,54:1780-1787.

半空间问题 篇4

许多工程实际问题都是抽象为搁置在弹性地基上的矩形板问题。代表性的工程如机场跑道面板、建筑物基础及其各种试验台。随着复合材料的广泛应用, 一些呈现正交各向异性的结构在实际工程中应用越来越广泛。因此, 研究弹性地基上的正交各向异性矩形板弯曲问题有重要的实际意义。

在分析地基板时, 常用的地基模型有Winkler模型、双参数模型和弹性半空间模型。而Winkler模型和双参数模型均需要对地基反力作出假设, 具有一定的局限性, 相对的弹性半空间模型则取消了地基反力的假设, 得到的解更加合理和精确。现有的研究, 要么采用Winkler弹性模型或双参数弹性模型;要么对地基反力做一些假设, 要么不能全部满足板的边界条件[2]等;要么用数值方法或半解析方法。现有很多大型商用软件可以求取弹性地基上的矩形板问题的数值解[3], 但是在许多情况下这个问题的解析解还是很有必要研究的。文献[4]研究了弹性半空间地基上四边自由矩形板的弯曲解析解, 但是均假定矩形板的材料是各向同性的, 而有关弹性半空间地基上正交各向异性矩形板的问题的解析解国内外研究甚少。

利用弹性半空间受任意竖向荷载作用下的静力积分变换解, 采用Fourier双三角正弦级数和单三角正弦级数, 分析了弹性半空间地基上四边自由正交各向异性矩形板的弯曲问题, 得出弹性半空间地基上四边自由正交各向异性矩形板受任意竖向荷载作用下的弯曲解析解。

1 控制微分方程及其一般解

弹性半空间地基上正交各向异性矩形板, 受垂直于板面横向力q (x, y) 作用, 地基作用于板的竖向反力为F (x, y) 时的控制微分方程, 若取x、y坐标轴同板的主方向平行时, 则有

式中, Dx、Dy分别为沿x、y方向的抗弯刚度, 它们同x、y方向的弹性模量Ex、Ey, 泊松比vx、vy和板的厚度h等的关系, 见文献[5]。

H=Dxvy+2Dxy为折算刚度, Dxy为抗扭刚度。

板的内力可以用挠度函数表示为

边界条件可表示为

根据控制微分方程 (1) , 可将其解表示成

其中w1, w2, w3为方程的齐次解;w0为特解。

这里 相应地将板面上的已知外荷载q (x, y) 和地基反力F (x, y) 也展为

将 (6) 、 (7) 代入方程 (1) , 比较系数, 得到

设矩形板四个角点的挠度分别为woo, wao, wob, wab;由式w3=Kxy+Lx+My+N得

讨论Ym (y) :

(1) H2

(2) H2=DxDy时, Ym (y) =A1mchεmy+B1mεmychεmy+

(3) H2>DxDy时, Ym (y) =A1mchθmy+B1mshθmy+

Xn (x) 与Ym (y) 类似, 这里不予列出。

为了简化计算, 设沿x=0, a;y=0, b各边的挠度均表示成带有补充项的正弦级数, 法向弯矩也表示成正弦级数, 即:

由文献[6]知:

(1) H2

(2) H2=DxDy时

(3) H2>DxDy时

以上三种情形, D12=Dxvy, w2分别与w1类似, 不予列出。

当弹性半空间地基上四边自由正交各向异性矩形板仅承受任意对称分布载荷q (x, y) 作用时, 由弯矩边界条件可知:En=Fn=Gm=Hm=0;由对称性知, woo=wao, wob=wab, 所以由式 (9) 知:K=L=M=0;同样由对称性可知:An=Bn, Cm=Dm。以下仅讨论H2=DxDy情形, H2DxDy情形类似, 不予列出。H2=DxDy时的w1、w2、w3简化成如下所示:

这样在式 (19) 、 (20) 、 (21) 中就只剩下待定系数Cm, An, N, Cmn。

2 求解方程的建立

方程的一般解为 (4) 式

同它们相应的内力素分别为Qx (1) 、Qx (2) 、Qx (3) 、Qx (0) 等等。利用剩下的剪力边界条件及角点条件, 令:在x=0边界上

在y=0边界上

在点 (0, 0) :

式 (25) 中, 进行了如下变换

其中

而1Θmn (2) 、1Θmn (3) 、1Θmn (4) 也有类似形式, 不予列出。

Qy与Qx类似, 这里不予列出。

若地基反力F (x, y) 已知, 由式 (8) 、 (22) ~ (24) 这四组方程可联立求解各待定系数Cm, An, N, Cmn, 各系数求得之后代回式 (4) 中得挠度, 进一步可求得弯矩等内力。若地基反力F (x, y) 未知, 先对F (x, y) 进行双重Fourier变换

再由 (7) 、 (35) 式得

由文献[4]知弹性半空间体受任意竖向荷载作用下的积分变换解为

式中, q12=ξ2+η2

将w|z=0看成是区域 (0≤x≤a, 0≤y≤b) 上的函数, 当然可将其展成双重正弦级数:

利用式 (37) 、 (39) , wzmn可以表示为

将板的挠度式也展成同样形式的双重正弦级数, 由于矩形薄板的挠度与弹性地基表面的竖向位移相等, 因此其对应项的系数也相等, 于是得以下变形协调方程:

式 (8) 、 (22) ~ (24) 以及 (41) 这5组方程, 即为弹性地基上矩形薄板弯曲的控制方程。联立求解各待定系数Cmn、Qmn、Cm, An, N。设m和n分别取最大M和N, 则共有2MN+M+N+1个方程, 用以解相同数目的未知数。各系数求得之后代回式 (4) , 得板的挠度, 进一步可求得弯矩内力值。

3数值算例

考虑一支承在弹性半空间上, 边长a=4m。厚度h=0.2m的弹性方薄板的弯曲。假设板与地基之间为光滑接触。地基泊松比为0.4, 弹性模量为343MPa。板的泊松比0.167, 弹性模量E=34300MPa。板上作用均布荷载q=0.98MPa, 采用本文的方法 (M=N=39) 、文献[4]的方法、四节点等参元的计算结果列于表1。

图1、图2分别是用本文方法所得到的矩形板的挠曲面图以及沿矩形板中心线和边界线的挠度, 图3、图4是弯矩Mx的三维图形及二维曲线。

由以上结果与文献[3]的结果对比可知, 本文的方法是切实可行的。

4结论

利用弹性半空间受任意竖向荷载作用下的静力积分变换解, 采用Fourier双三角正弦级数和单三角正弦级数, 分析了弹性半空间地基上四边自由正交各向异性矩形板的弯曲问题, 得出弹性半空间地基上四边自由正交各向异性矩形板受任意竖向荷载作用下的弯曲解析解, 包括板的挠度、内力。不仅弥补了数值法的一些不足, 同时取消了Winkler地基模型或双参数地基模型的假设, 从而得到板的内力更合理、更精确的分布规律。同样的, 也可利用挠度表达式的一般格式求解弹性半空间地基上任意边界条件下的正交各向异性矩形板受任意竖向荷载作用下的弯曲解析解。

参考文献

[1]Eatock Taylor E, Ohkusu M.Green function for hydroelastic analysis of vibrating free-free beams and plates[J].Applied Ocean Research, 2000, 22 (1) :12-19.

[2]曲庆璋, 梁兴复, 董学江.Winkler弹性地基上自由矩形板问题[J].强度与环境, 1998 (2) :31-37.

[3]王元汉, 邱先敏, 张佑启.弹性地基板的等参有限元法计算[J].岩土工程学报, 1998, 20 (4) :7-11.

[4]王春玲, 黄义.弹性半空间地基上四边自由矩形板的弯曲解析解[J].岩土工程学报, 2005, 27 (12) :1402-1407.

[5]曲庆璋, 章权, 季求知等.弹性板理论[M].北京:人民交通出版社, 2000.

[6]程银水, 朱钅奇.正交各向异性矩形薄板弯曲问题解析解的一般格式[J].北京建筑工程学院学报, 1991 (1) :9-18.

[7]李刚, 王贻荪, 尚守平.弹性半空间上矩形基础稳态振动积分变换解[J].湖南大学学报, 2000, 27 (4) :88-93.