理解波动（二）：以有限计算无限——人工边界条件

发表于2015-12-25由范雨

在理解波动（一）：在有限结构上算出“波动”中，我们展示了如何简单地通过常用的有限元软件计算出结构的“波动现象”。可以发现，在由激振力的引起的扰动抵达边界前，清晰地呈现出“波传播”的特性。而在经历了边界的数次反射后，由于整个结构各处都有了变形，直观上更容易按照“自由振动”的特性来理解。在上一个算例中，我们初步探讨了“波动特性”是如何蕴含在有限结构中的。在本文中，我们将展示，“无限”是如何蕴含在“有限”中的。即如何通过修饰有限结构来计算无限结构的特性。由于没有边界，无限结构是研究波动的更理想模型。无限大结构虽然在实际中不严格存在，但是它仍有工程应用场合，如地震研究、隔声、无损探伤等。总的来说，在这些应用场合中：1）只关心到波未被边界反射之前的行为；2）由于阻尼或结构太大，波很难抵达边界。

人工边界条件（ABC）简介

熟悉有限元的工程师都明白，有限元分析的第一步是确定要分析的结构，及其边界条件。如果有有限大的结构，边界是自然存在的，可以是结构自身与其它结构无连接又不受力的部分（自由边界），与其他结构有连接的部分可以处理为力边界或者力-位移混合边界。但是“无限大”结构在其无限延展的那个维度上不具备“边界”，怎么办？我们可以采用一种类似于子结构的办法，如下图所示：
bg
我们在无限大的结构上划出一个有限大的部分，称为“近场”，包含所有感兴趣的部分。将剩下的部分（半无限大）划分为另一个（或多个，对于上图是左右两个）子结构，称为“远场”。这样，问题转换为，怎样对“近场”和“远场”建模并把它们整合成原始的无限大结构。

对于近场的建模相当于对一个具有自由边界的有限结构进行建模，并不难，可用的方法有很多，见上图，最常用的还是有限元方法。而对于远场的建模相当于对半无界结构的建模，目前并未有一个适用于所有情况下的方法，根据不同的场合，可以应用边界元、无限元、有限差分，以及这里要说的人工边界条件。

人工边界条件的英文是Artificial Boundary Condition，简称ABC；有的场合又称为Perfectly Matched Layer，简称PML。有针对稳态响应的频域ABC和针对瞬态响应的时域ABC。根据建立思路的不同又分为位移型ABC和应力型ABC。这里我们说最简单也是最有效的频域应力型ABC。

频域应力型ABC的推导思路

它的思路很简单，其最终目的是表达出近场与远场交界面处的力/位移关系，写作：
{\bf H}_\text{f}{\bf q}_\text{b} = -{\bf f}_\text{b}
近场通过有限元建模，其（频域）动力学方程为：
\begin{bmatrix}{\bf H}_\text{ii} & {\bf H}_\text{ib} \\{\bf H}_\text{bi} & {\bf H}_\text{bb} \end{bmatrix}\left(\begin{array}{c}{\bf q}_\text{i}\\ {\bf q}_\text{b}\end{array}\right)=\left(\begin{array}{c}{\bf f}_\text{i}\\ {\bf f}_\text{b}\end{array}\right)
将远场动力学方程带入近场中，消去边界处的内力{\bf f}_\text{b}，得：
\begin{bmatrix}{\bf H}_\text{ii} & {\bf H}_\text{ib} \\{\bf H}_\text{bi} & {\bf H}_\text{bb}+{\bf H}_\text{f} \end{bmatrix}\left(\begin{array}{c}{\bf q}_\text{i}\\ {\bf q}_\text{b}\end{array}\right)=\left(\begin{array}{c}{\bf f}_\text{i}\\ {\bf 0}\end{array}\right)

这样，就获得了原始无界结构的有限元模型，可以看到模型的自由度数和近场模型相当，远场的影响被考虑为一个边界传递函数。那么，问题的关键就是传递函数{\bf H}_\text{f}的获得。由于远场通常是均匀的，因此可以用解析模型推导出当波动到达边界时的位移和应力，从而获得传递函数。接下来分别用杆的一维杆方程和欧拉伯努力梁方程为例，展示人工边界条件的大致推导流程。

一维杆的ABC

对于一个原点受集中简谐力的杆，其动力学方程为：
\rho A\frac{\partial^2 u}{\partial t^2}-EA\frac{\partial^2 u}{\partial x^2}=P\delta_0e^{-j\omega t}
采用频域指标：
u(x,t)=U(x)e^{-j\omega t}
其频域动力学方程为：
\frac{\text{d}^2 u}{\text{d} x^2}+k^2U=-\frac{P}{EA}\delta_0
其中k = \omega \sqrt{\rho /E}是波数。对于这个问题求解的方式有很多，最简单粗暴的是作空间傅里叶变换，求解代数方程组后再作空间傅立叶反变换。或者直接采用达朗贝尔解。在正方向上的一点，其变形为：
U(x,x>0)= Ce^{jkx}
在严格的求解过程中，C的值是可以完全确定的，但是这里它并不重要，故没有写出详细形式。其内力为：
F(x,x>0)= -EA\frac{\text{d}U}{\text{d}x}=-jEAkU(x)
这样就获得了在任意点的频域传递函数：
H=jEAk=-j\omega A\sqrt{E\rho}

不难发现，在一维杆的情形下，无限远场的传函等效于一个阻尼，其阻尼系数A\sqrt{E\rho}为常值。这就为ABC与有限元的结合实现提供了极大的便利——即在近场的边界上设置阻尼单元就可以了！由于传递函数是常数，所以适用于各种频率，因此标准波动方程的频域ABC也是时域ABC。

上述推导和结论实际上适用于所有可以通过标准波动方程描述的现象，如弦的振动、扭转、声辐射等等。

欧拉伯努力梁的ABC

欧拉伯努力梁的频域动力学方程：
\frac{\text{d}^4U}{\text{d}x^4}-k^4U=P\delta_0
其中，波数：
k=\sqrt[4]{\omega^2\frac{\rho A}{EI}}
同样，可以用空间傅里叶变换的方法求解，其正方向上一点的位移为：
U(x)=Ce^{jkx}+De^{-kx}
转角为：
\theta(x)=\frac{\text{d}U}{\text{d}x}=jkCe^{jkx}-kDe^{-kx}
弯矩为：
M(x)=-EI\frac{\text{d}^2U}{\text{d}x^2}=EIk^2Ce^{jkx}-EIk^2De^{-kx}
内力为：
F(x)=-EI\frac{\text{d}^3U}{\text{d}x^3}=jEIk^3Ce^{jkx}+EIk^3De^{-kx}
注意到上述四个分量中都有一个随着空间快速衰减的分量（与e^{-kx}相关的那些项），可以将其忽略，之后得：
位移-内力传递函数：
H_\text{1}=\frac{F}{U}=-jEIk^3
转角-弯矩传递函数：
H_\text{2}=\frac{M}{\theta}=-jEIk

这样，就获得了欧拉伯努力梁的ABC。由于梁单元的每个节点是有两个自由度，所以获得了两个传递函数。在推导过程中忽略了快衰波的影响，因此要求ABC离激振力较远（至少2倍波长）。同样可以发现，无限长的梁对于近场的作用也相当于某种“阻尼”，相当于一个横向阻尼系数为EIk^3/\omega、扭转阻尼系数为EIk/\omega的阻尼器。但是阻尼系数与频率有关，因此仅适合用于频域或激振力只有一种频率成分的时域分析。

一个算例

采用理解波动（一）中同样的梁，在其两端设置人工边界条件,使其成为两端无界的均匀梁。然后再于1个固定的频率下激励，其结果为：

ABC

从图中可以清晰地观察到波动行为，而且可以看出，波在边界处并没有被反射，而似乎是“透过”了边界。仔细观察还能发现激振点附近的响应要稍稍大一些，那实际上是快衰波的作用。最后，通过上述推导也不难发现，远场对于近场来说实际上相当于一种“耗能”机制。这样，由近场流向交界面的能量才能“有去无回”。

讨论

在无限结构中截取有限个部分，并在其边界上施加人工边界条件，就精确地模拟无限结构的动力学行为。人工边界条件的建立需要事先知道无限远场在边界处的传递函数。这可以通过解析或者数值获得。本文举了杆和梁的例子，并用ANSYS为例说明了所推导的ABC的正确定。其实这些只是人工边界条件研究中最基础的一些观点。更大的挑战在于实现1）多维情况下的ABC和2）时域ABC。相关参考文献将会集中列在本系列文章的最后一篇中。

附：APDL代码


FINISH
/CLEAR
!作者：范雨
!电子邮件：fanyu.buaa@gmail.com
/PREP7
*SET,N_CIRCLE,40 !计算的周期数
*SET,NSTP_C,30 !每周期的载荷步数
*SET,FRE,10000 !载荷频率
*SET,CURR_T,0
*SET,DT,1.0/FRE/NSTP_C
*SET,PI,3.141592653

/PREP7
/eshape,1
!近场
BL = 1
BW = 2e-2
BH = 2e-2
NMESH = 100
ET,1,BEAM188
MP,EX,1,2.1E11
MP,PRXY,1,0.3
MP,DENS,1,7800
MAT,1
SECTYPE,1,BEAM,RECT
SECDATA,BW,BH,2,2
K,1,-BL,0,0
K,3,0,0,0
K,2,BL,0,0
L,2,3,NMESH
L,1,3,NMESH
LMESH,ALL

!ABC区域
N,1000,BL*1.1,0,0
N,1001,-BL*1.1,0,0

VI = BW*BH*BH*BH/12
VOMEGA = 2*PI*FRE
VK = SQRT(VOMEGA*VOMEGA*7800*BW*BH/2.1E11/VI)
VK = SQRT(VK)

!横向阻尼
ET,2,COMBIN14,0,2
R,2,0,2.1E11*VI*VK*VK*VK/VOMEGA
REAL,2
TYPE,2
E,1000,1
E,1001,102

!扭转阻尼
ET,3,COMBIN14,0,6
R,3,0,2.1E11*VI*VK/VOMEGA
REAL,3
TYPE,3
E,1000,1
E,1001,102

ALLSEL,ALL
D,1000,ALL,0
D,1001,ALL,0
ALLSEL,ALL
/SOLU
NSEL,S,LOC,X,0
F,ALL,FY,1
ALLSEL,ALL
ANTYPE,HARM
HARFRQ,FRE
NSUBST,1
KBC,1
HROPT,FULL
HROUT,OFF
LUMPM,OFF
EQSLV,,1E-8
PSTRES,0
SOLVE
FINISH
未知 的头像

About 范雨

Talk is cheap. Show me the code.
此条目发表在ANSYS, Dynamics, 有限元, 波动分类目录。将固定链接加入收藏夹。

1 Response to 理解波动（二）：以有限计算无限——人工边界条件

  1. Pingback： 理解波动（三）：从频域到时域 | Dynamics and Python

留下评论