边界元方法(一)

前言

边界元法是在经典积分方程法和有限元方法基础上发展起来的一种偏微分方程数值解法。 Brebbia 于 1978 年用加权余量法统一边界积分方程的不同推导,并首次以“边界元法”这一名称统一称呼此前边界积分方程离散化的多种不同方法。

边界元素法的基本思路是:将计算区域内的微分方程变换为区域边界上的积分方程;然后将区域边界剖分为有限大小的边界元素,从而将边界积分方程离散化为线性代数方程组,并经求解线性方程组而得到边界节点物理量。

边界元素法与有限差分和有限元素法相比,具有以下显著优点:

  1. 计算精度高:边界元素法直接求解边界广义场源的分布。边界上任一点场量将通过线性叠加各离散广义场源的作用而得到,不必经微分计算间接得到。
  2. 计算速度快:由于仅对边界离散化,它降低了问题的维数,不必如有限差分及有限元法那样需引进大量区域变量,使变量数极大减少。
  3. 处理复杂边界能力强:由于只离散边界,大大降低了几何处理复杂性,从而增强了处理复杂边界的能力。

不过,由于边界元计算形成的线性方程组系数矩阵一般是稠密、不对称的,随着集成密度指数级增加,问题规模不断扩大,给方程求解带来越来越大的困难,特别是对三维问题。

Green 公式

设闭区域 D D D 由分段光滑的曲线 L L L 围成,函数 P ( x , y ) P(x, y) P(x,y) Q ( x , y ) Q(x, y) Q(x,y) D D D 上具有一阶连续偏导数,则有
∬ D ( ∂ Q ∂ x − ∂ P ∂ y ) d x d y = ∮ L P d x + Q d y \iint_D (\frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y}) dxdy = \oint_{L} Pdx + Qdy D(xQyP)dxdy=LPdx+Qdy
其中 L L L D D D 的取正向的边界曲线。

Green 第一公式

∬ S [ ( ∂ u ∂ x ) 2 + ( ∂ u ∂ y ) 2 ] d x d y = − ∬ s u Δ u d x d y + ∮ c u ∂ u ∂ n → d s \iint_{\mathrm{S}}\left[\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)^2+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)^2\right] \mathrm{dxdy}=-\iint_{\mathrm{s}} \mathrm{u} \Delta \mathrm{udxdy}+\oint_{\mathrm{c}} \mathrm{u} \frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{ds} S[(xu)2+(yu)2]dxdy=suΔudxdy+cun uds

证明:

不妨设 n → = ( cos ⁡ θ , sin ⁡ θ ) \overrightarrow{\mathrm{n}}=(\cos \theta, \sin \theta) n =(cosθ,sinθ);

由方向导数的定义有:
∂ u ∂ n → = ∂ u ∂ x cos ⁡ θ + ∂ u ∂ y sin ⁡ θ \frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}}=\frac{\partial \mathrm{u}}{\partial \mathrm{x}} \cos \theta+\frac{\partial \mathrm{u}}{\partial \mathrm{y}} \sin \theta n u=xucosθ+yusinθ
可知有
cos ⁡ θ = d y ( d x ) 2 + ( d y ) 2 sin ⁡ θ = − d x ( d x ) 2 + ( d y ) 2 d s = ( d x ) 2 + ( d y ) 2 \begin{aligned} \cos \theta &=\frac{d y}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}} \\ \sin \theta &=\frac{-\mathrm{dx}}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}} \\ \mathbf{d s} &=\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2} \end{aligned} cosθsinθds=(dx)2+(dy)2 dy=(dx)2+(dy)2 dx=(dx)2+(dy)2
故有
∮ c u ∂ u ∂ n → d s = ∮ c u ( ∂ u ∂ x d y ( d x ) 2 + ( d y ) 2 + ∂ u ∂ y d y ( d x ) 2 + ( d y ) 2 ) ( d x ) 2 + ( d y ) 2 = ∮ c u ∂ u ∂ x d y − u ∂ u ∂ y d x \begin{gathered} \oint_{\mathrm{c}} \mathrm{u} \frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{ds} \\ =\oint_{\mathrm{c}} \mathrm{u}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}} \frac{\mathrm{dy}}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}}+\frac{\partial \mathrm{u}}{\partial \mathrm{y}} \frac{\mathrm{dy}}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}}\right) \sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2} \\ =\oint_{\mathrm{c}} \mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{x}} \mathrm{dy}-\mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{y}} \mathrm{dx} \end{gathered} cun uds=cu(xu(dx)2+(dy)2 dy+yu(dx)2+(dy)2 dy)(dx)2+(dy)2 =cuxudyuyudx

由 Green 公式 :
∬ D ( ∂ Q ∂ x − ∂ P ∂ y ) d x d y = ∮ ∂ D P d x + Q d y \iint_{\mathrm{D}}\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{x}}-\frac{\partial \mathrm{P}}{\partial \mathrm{y}}\right) \mathrm{dxdy}=\oint_{\partial \mathrm{D}} \mathrm{Pdx}+\mathrm{Qdy} D(xQyP)dxdy=DPdx+Qdy

可以得到

∮ c u ∂ u ∂ x d y − u ∂ u ∂ y d x = ∬ S [ ∂ ∂ x ( u ∂ u ∂ x ) − ∂ ∂ y ( − u ∂ u ∂ y ) ] d x d y = ∬ S [ ∂ ∂ x ( u ∂ u ∂ x ) + ∂ ∂ y ( u ∂ u ∂ y ) ] d x d y = ∬ S [ ∂ ∂ x ( ∂ u ∂ x ) u + ( ∂ u ∂ x ) 2 + ∂ ∂ y ( ∂ u ∂ y ) u + ( ∂ u ∂ y ) 2 ] d x d y = ∬ S [ ( ∂ u ∂ x ) 2 + ( ∂ u ∂ y ) 2 ] d x d y + ∬ S [ ∂ ∂ x ( ∂ u ∂ x ) u + ∂ ∂ y ( ∂ u ∂ y ) u ] d x d y = ∬ S [ ( ∂ u ∂ x ) 2 + ( ∂ u ∂ y ) 2 ] d x d y + ∬ S u [ ∂ ∂ x ( ∂ u ∂ x ) + ∂ ∂ y ( ∂ u ∂ y ) ] d x d y = ∬ S [ ( ∂ u ∂ x ) 2 + ( ∂ u ∂ y ) 2 ] d x d y + ∬ S u Δ u d x d y \begin{aligned} &\oint_{\mathrm{c}} \mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{x}} \mathrm{dy}-\mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{y}} \mathrm{dx} \\ &=\iint_{\mathrm{S}}\left[\frac{\partial}{\partial \mathrm{x}}\left(\mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)-\frac{\partial}{\partial \mathrm{y}}\left(-\mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)\right] \mathrm{dxdy} \\ &=\iint_{\mathrm{S}}\left[\frac{\partial}{\partial \mathrm{x}}\left(\mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)+\frac{\partial}{\partial \mathrm{y}}\left(\mathrm{u} \frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)\right] \mathrm{dxdy} \\ &=\iint_{\mathrm{S}}\left[\frac{\partial}{\partial \mathrm{x}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right) \mathrm{u}+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)^2+\frac{\partial}{\partial \mathrm{y}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right) \mathrm{u}+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)^2\right] \mathrm{dxdy} \\ &=\iint_{\mathrm{S}}\left[\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)^2+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)^2\right] \mathrm{dxdy}+\iint_{\mathrm{S}}\left[\frac{\partial}{\partial \mathrm{x}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right) \mathrm{u}+\frac{\partial}{\partial \mathrm{y}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right) \mathrm{u}\right] \mathrm{dxdy} \\ &=\iint_{\mathrm{S}}\left[\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)^2+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)^2\right] \mathrm{dxdy}+\iint_{\mathrm{S}} \mathrm{u}\left[\frac{\partial}{\partial \mathrm{x}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)+\frac{\partial}{\partial \mathrm{y}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)\right] \mathrm{dxdy} \\ &=\iint_{\mathrm{S}}\left[\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)^2+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)^2\right] \mathrm{dxdy}+\iint_{\mathrm{S}} \mathrm{u} \Delta \mathrm{udxdy} \end{aligned} cuxudyuyudx=S[x(uxu)y(uyu)]dxdy=S[x(uxu)+y(uyu)]dxdy=S[x(xu)u+(xu)2+y(yu)u+(yu)2]dxdy=S[(xu)2+(yu)2]dxdy+S[x(xu)u+y(yu)u]dxdy=S[(xu)2+(yu)2]dxdy+Su[x(xu)+y(yu)]dxdy=S[(xu)2+(yu)2]dxdy+SuΔudxdy
即有
∮ c u ∂ u ∂ n → d s = ∬ S [ ( ∂ u ∂ x ) 2 + ( ∂ u ∂ y ) 2 ] d x d y + ∬ S u Δ u d x d y \oint_{\mathrm{c}} \mathrm{u} \frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{d} s=\iint_{\mathrm{S}}\left[\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)^2+\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)^2\right] \mathrm{dxdy}+\iint_{\mathrm{S}} u \Delta \mathrm{udxdy} cun uds=S[(xu)2+(yu)2]dxdy+SuΔudxdy
移项可得原式。

Green 第二公式

∬ S ∣ Δ u Δ v u v ∣ d x d y = ∮ C ∣ ∂ u ∂ n ⃗ ∂ v ∂ n ⃗ v v ∣ d s \iint_S\left|\begin{array}{cc} \Delta u & \Delta v \\ u & v \end{array}\right| d x d y=\oint_C\left|\begin{array}{ll} \frac{\partial u}{\partial \vec{n}} & \frac{\partial v}{\partial \vec{n}} \\ v & v \end{array}\right| ds S ΔuuΔvv dxdy=C n uvn vv ds

证明:

左边:
∬ S ∣ Δ u Δ v u v ∣ d x d y = ∬ S v Δ u − u Δ v d x d y \iint_S\left|\begin{array}{cc} \Delta u & \Delta v \\ u & v \end{array}\right| d x d y=\iint_S v \Delta u-u \Delta v d x d y S ΔuuΔvv dxdy=SvΔuuΔvdxdy
右边
∮ C ∣ ∂ u ∂ n → ∂ v ∂ n → ∣ d s \oint_{\mathrm{C}}\left|\begin{array}{cc}\frac{\partial \mathrm{u}}{\overrightarrow{\partial \mathrm{n}}} & \frac{\partial \mathrm{v}}{\partial \overrightarrow{\mathrm{n}}}\end{array}\right| \mathrm{ds} C n un v ds
= ∮ C ( ∂ u ∂ n → v − ∂ v ∂ n → u ) d s =\oint_C\left(\frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{v}-\frac{\partial \mathrm{v}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{u}\right) \mathrm{ds} =C(n uvn vu)ds
= ∮ C v ∂ u ∂ x d y ( d x ) 2 + ( d y ) 2 − v ∂ u ∂ y d x ( d x ) 2 + ( d y ) 2 =\oint_C \mathrm{v} \frac{\partial \mathrm{u}}{\partial \mathrm{x}} \frac{\mathrm{dy}}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}}-\mathrm{v} \frac{\partial \mathrm{u}}{\partial \mathrm{y}} \frac{\mathrm{dx}}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}} =Cvxu(dx)2+(dy)2 dyvyu(dx)2+(dy)2 dx
− u ∂ v ∂ x d y ( d x ) 2 + ( d y ) 2 -u \frac{\partial v}{\partial x} \frac{d y}{\sqrt{(d x)^2+(d y)^2}} uxv(dx)2+(dy)2 dy
+ u ∂ v ∂ y d x ( d x ) 2 + ( d y ) 2 ( d x ) 2 + ( d y ) 2 +u \frac{\partial v}{\partial y} \frac{d x}{\sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2}} \sqrt{(\mathbf{d x})^2+(\mathbf{d y})^2} +uyv(dx)2+(dy)2 dx(dx)2+(dy)2
= ∮ C v ∂ u ∂ x d y − v ∂ u ∂ y d x − u ∂ v ∂ x d y + u ∂ v ∂ y d x =\oint_C v \frac{\partial u}{\partial x} d y-v \frac{\partial u}{\partial y} d x-u \frac{\partial v}{\partial x} d y+u \frac{\partial v}{\partial y} d x =Cvxudyvyudxuxvdy+uyvdx
= ∮ C ( u ∂ v ∂ y − v ∂ u ∂ y ) d x + ( v ∂ u ∂ x − u ∂ v ∂ x ) d y =\oint_C\left(u \frac{\partial v}{\partial y}-v \frac{\partial u}{\partial y}\right) d x+\left(v \frac{\partial u}{\partial x}-u \frac{\partial v}{\partial x}\right) d y =C(uyvvyu)dx+(vxuuxv)dy
由 Green 公式:
∬ D ( ∂ Q ∂ x − ∂ P ∂ y ) d x d y = ∮ ∂ D P d x + Q d y ; \iint_{\mathrm{D}}\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{x}}-\frac{\partial \mathrm{P}}{\partial \mathrm{y}}\right) \mathrm{dxdy}=\oint_{\partial \mathrm{D}} \mathrm{Pdx}+\mathrm{Qdy} ; D(xQyP)dxdy=DPdx+Qdy;

P = ( u ∂ v ∂ y − v ∂ u ∂ y ) Q = ( v ∂ u ∂ x − u ∂ v ∂ x ) ∂ Q ∂ x = ∂ ( v ∂ u ∂ x − u ∂ v ∂ x ) ∂ x = ∂ v ∂ x ∂ u ∂ x + v ∂ 2 u ∂ x 2 − ∂ v ∂ x ∂ u ∂ x − u ∂ 2 v ∂ x 2 = v ∂ 2 u ∂ x 2 − u ∂ 2 v ∂ x 2 \begin{aligned} &\mathrm{P}=\left(\mathrm{u} \frac{\partial \mathrm{v}}{\partial \mathrm{y}}-\mathrm{v} \frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right) \\ &\mathrm{Q}=\left(\mathrm{v} \frac{\partial \mathrm{u}}{\partial \mathrm{x}}-\mathrm{u} \frac{\partial \mathrm{v}}{\partial \mathrm{x}}\right) \\ &\frac{\partial \mathrm{Q}}{\partial \mathrm{x}}=\frac{\partial\left(\mathrm{v} \frac{\partial \mathrm{u}}{\partial \mathrm{x}}-\mathrm{u} \frac{\partial \mathrm{v}}{\partial \mathrm{x}}\right)}{\partial \mathrm{x}} \\ &=\frac{\partial v}{\partial x} \frac{\partial u}{\partial x}+v \frac{\partial^2 u}{\partial x^2}- \frac{\partial v}{\partial x} \frac{\partial u}{\partial x}-u \frac{\partial^2 v}{\partial x^2} \\ &=v \frac{\partial^2 u}{\partial x^2}-u \frac{\partial^2 v}{\partial x^2} \end{aligned} P=(uyvvyu)Q=(vxuuxv)xQ=x(vxuuxv)=xvxu+vx22uxvxuux22v=vx22uux22v

同理

∂ P ∂ y = u ∂ 2 v ∂ y 2 − v ∂ 2 u ∂ y 2 \frac{\partial \mathrm{P}}{\partial \mathrm{y}}=\mathrm{u} \frac{\partial^{2} \mathrm{v}}{\partial \mathrm{y}^{2}}-\mathrm{v} \frac{\partial^{2} \mathrm{u}}{\partial \mathrm{y}^{2}} yP=uy22vvy22u

故有:

∬ D ( ∂ Q ∂ x − ∂ P ∂ y ) d x d y = ∬ D ( v ∂ 2 u ∂ x 2 − u ∂ 2 v ∂ x 2 − u ∂ 2 v ∂ y 2 + v ∂ 2 u ∂ y 2 ) d x d y = ∬ D v Δ u − u Δ v d x d y = ∬ S ∣ Δ u Δ v u v ∣ d x d y \begin{array}{l} \iint_{\mathrm{D}}\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{x}}-\frac{\partial \mathrm{P}}{\partial \mathrm{y}}\right) \mathrm{dxdy} \\ =\iint_{\mathrm{D}}\left(\mathrm{v} \frac{\partial^{2} \mathrm{u}}{\partial \mathrm{x}^{2}}-\mathrm{u} \frac{\partial^{2} \mathrm{v}}{\partial \mathrm{x}^{2}}-\mathrm{u} \frac{\partial^{2} \mathrm{v}}{\partial \mathrm{y}^{2}}+\mathrm{v} \frac{\partial^{2} u}{\partial \mathrm{y}^{2}}\right) \mathrm{dxdy} \\ =\iint_{\mathrm{D}} \mathrm{v} \Delta \mathrm{u}-\mathrm{u} \Delta \mathrm{v} d x d y=\iint_{\mathrm{S}}\left|\begin{array}{cc} \Delta \mathrm{u} & \Delta \mathrm{v} \\ \mathrm{u} & \mathrm{v} \end{array}\right| \mathrm{dx} d \mathrm{y} \end{array} D(xQyP)dxdy=D(vx22uux22vuy22v+vy22u)dxdy=DvΔuuΔvdxdy=S ΔuuΔvv dxdy

Green 第三公式

若 u 为有界闭区域 S 中的调和函数, 则有: u ( x , y ) = 1 2 π ∮ C ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s u(x, y)=\frac{1}{2 \pi} \oint_{C}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s u(x,y)=2π1C(un lnrlnrn u)ds 其中 C C C S S S 边际, ∂ u ∂ n ⃗ \frac{\partial u}{\partial \vec{n}} n u u u u 沿着 C C C 的外法线方向的方向导数; r = ( ξ − x ) 2 + ( η − y ) 2 r=\sqrt{(\xi-x)^{2}+(\eta-y)^{2}} r=(ξx)2+(ηy)2 ( x , y ) (x, y) (x,y) 到边界 C C C 上动点 ( ξ , η ) (\xi, \eta) (ξ,η) 的距离;
证明:
由 Green 第二公式得到

∮ C ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s = ∬ D v Δ u − u Δ v d x d y \oint_{C}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s=\iint_{D} v \Delta u-u \Delta v d x d y C(un lnrlnrn u)ds=DvΔuuΔvdxdy

由于 u u u 为有界闭区域 S \mathrm{S} S 中的调和函数, Δ u = 0 \Delta \mathrm{u}=0 Δu=0

Δ v = Δ ln ⁡ r = Δ ln ⁡ ( ξ − x ) 2 + ( η − y ) 2 = 0 \Delta v=\Delta \ln r=\Delta \ln \sqrt{(\xi-x)^{2}+(\eta-y)^{2}}=0 Δv=Δlnr=Δln(ξx)2+(ηy)2 =0

可知 ln ⁡ r \ln \mathrm{r} lnr 也是调和函数,
故在没有奇点的情况下, S \mathrm{S} S 内的任何区域
∮ C ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s = ∬ D v Δ u − u Δ v d x d y = 0 \oint_{C}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s=\iint_{D} v \Delta u-u \Delta v d x d y=0 C(un lnrlnrn u)ds=DvΔuuΔvdxdy=0
故有设以 ( x , y ) (x, y) (x,y) 为中心, t t t 为半径的一个领域 D D D ,
∮ C ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s = ∮ ∂ D ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s \oint_{C}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s=\oint_{\partial D}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s C(un lnrlnrn u)ds=D(un lnrlnrn u)ds
有在 ∂ D \partial \mathrm{D} D 上,
∮ ∂ D ln ⁡ r ∂ u ∂ n → d s = ln ⁡ t ∮ ∂ D ∂ u ∂ n → d s = ln ⁡ t ∬ D Δ u d s = 0 ∮ ∂ D u ∂ ln ⁡ r ∂ n → d s = ∮ ∂ D u ∂ ln ⁡ r ∂ r d s = ∮ ∂ D u 1 t d s = 1 t ∮ ∂ D u d s = 2 π u ( ξ 1 , η 1 ) \begin{array}{c} \oint_{\partial \mathrm{D}} \ln r \frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{ds}=\ln \mathrm{t} \oint_{\partial \mathrm{D}} \frac{\partial \mathrm{u}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{ds}=\ln \mathrm{t} \iint_{\mathrm{D}} \Delta \mathrm{uds}=0 \\ \\ \oint_{\partial \mathrm{D}} \mathrm{u} \frac{\partial \ln \mathrm{r}}{\partial \overrightarrow{\mathrm{n}}} \mathrm{ds}=\oint_{\partial \mathrm{D}} u \frac{\partial \ln \mathrm{r}}{\partial \mathrm{r}} \mathrm{ds}=\oint_{\partial \mathrm{D}} \mathrm{u} \frac{1}{\mathrm{t}} \mathrm{d} s=\frac{1}{\mathrm{t}} \oint_{\partial \mathrm{D}} \mathrm{uds}=2 \pi \mathrm{u}\left(\xi_{1}, \eta_{1}\right) \end{array} Dlnrn uds=lntDn uds=lntDΔuds=0Dun lnrds=Durlnrds=Dut1ds=t1Duds=2πu(ξ1,η1)

故由 u u u S S S 上的连续性得到
lim ⁡ t → 0 ∮ C ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s = lim ⁡ t → 0 2 π u ( ξ 1 , η 1 ) = 2 π u ( x , y ) \lim _{t \rightarrow 0} \oint_{C}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s=\lim _{t \rightarrow 0} 2 \pi u\left(\xi_{1}, \eta_{1}\right)=2 \pi u(x, y) t0limC(un lnrlnrn u)ds=t0lim2πu(ξ1,η1)=2πu(x,y)
故得证
u ( x , y ) = 1 2 π ∮ C ( u ∂ ln ⁡ r ∂ n ⃗ − ln ⁡ r ∂ u ∂ n ⃗ ) d s u(x, y)=\frac{1}{2 \pi} \oint_{C}\left(u \frac{\partial \ln r}{\partial \vec{n}}-\ln r \frac{\partial u}{\partial \vec{n}}\right) d s u(x,y)=2π1C(un lnrlnrn u)ds

二重积分的分步积分公式

∫ Ω ∇ ⋅ ( c ∇ u ) v d x d y = ∫ ∂ Ω ( c ∇ u ⋅ n ⃗ ) v d s − ∫ Ω c ∇ u ⋅ ∇ v d x d y \int_{\Omega} \nabla \cdot (c \nabla u)v dxdy = \int_{\partial \Omega} (c \nabla u \cdot \vec{n})v ds -\int_{\Omega} c \nabla u \cdot \nabla v dxdy Ω(cu)vdxdy=Ω(cun )vdsΩcuvdxdy

∬ Ω f ∂ g ∂ x d x d y = ∮ ∂ Ω ( f g ) d y − ∬ Ω g ∂ f ∂ x d x d y \iint_{\Omega} f \frac{\partial g}{\partial x} dxdy = \oint_{\partial \Omega} (fg) dy - \iint_{\Omega} g \frac{\partial f}{\partial x} dxdy Ωfxgdxdy=Ω(fg)dyΩgxfdxdy

∬ Ω f ∂ g ∂ y d x d y = − ∮ ∂ Ω ( f g ) d x − ∬ Ω g ∂ f ∂ y d x d y \iint_{\Omega} f \frac{\partial g}{\partial y} dxdy = - \oint_{\partial \Omega} (fg) dx - \iint_{\Omega} g \frac{\partial f}{\partial y} dxdy Ωfygdxdy=Ω(fg)dxΩgyfdxdy

证明:(由Green公式) ∬ Ω ∂ F ∂ x d x d y = ∮ ∂ Ω F d y \iint_{\Omega} \frac{\partial F}{\partial x} dxdy = \oint_{\partial \Omega} Fdy ΩxFdxdy=ΩFdy , 把 F F F 换为 f g fg fg,则

∬ Ω ∂ ( f g ) ∂ x d x d y = ∮ ∂ Ω ( f g ) d y \iint_{\Omega} \frac{\partial (fg)}{\partial x} dxdy = \oint_{\partial \Omega} (fg)dy Ωx(fg)dxdy=Ω(fg)dy

∬ Ω ( f ∂ g ∂ x + g ∂ f ∂ x ) d x d y = ∮ ∂ Ω ( f g ) d y \iint_{\Omega} (f \frac{\partial g}{\partial x} + g \frac{\partial f}{\partial x}) dxdy = \oint_{\partial \Omega} (fg)dy Ω(fxg+gxf)dxdy=Ω(fg)dy

∬ Ω f ∂ g ∂ x d x d y = ∮ ∂ Ω ( f g ) d y − ∬ Ω g ∂ f ∂ x d x d y \iint_{\Omega} f \frac{\partial g}{\partial x} dxdy = \oint_{\partial \Omega} (fg)dy - \iint_{\Omega} g \frac{\partial f}{\partial x} dxdy Ωfxgdxdy=Ω(fg)dyΩgxfdxdy

注: ∫ u v ′ d x = u v − ∫ v u ′ d x \int u v'dx = uv - \int vu'dx uvdx=uvvudx

积分方程

以二维的 Poisson 方程为例将其转换成积分方程。

∇ 2 u ( x , y ) = f ( x , y ) ,       i n    Ω            ( 1 ) \nabla^2 u(x, y) = f(x, y) , \ \ \ \ \ in \ \ \Omega \ \ \ \ \ \ \ \ \ \ (1) 2u(x,y)=f(x,y),     in  Ω          (1)
u ( x , y ) = u 0 ,       o n    Γ 1                     ( 2 ) u(x, y) = u_0 , \ \ \ \ \ on \ \ \Gamma_1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2) u(x,y)=u0,     on  Γ1                   (2)
q ( x , y ) = q 0 ,       o n    Γ 2                     ( 3 ) q(x, y) = q_0 , \ \ \ \ \ on \ \ \Gamma_2 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3) q(x,y)=q0,     on  Γ2                   (3)

其中 ∂ Ω = Γ 1 ∪ Γ 2 \partial\Omega = \Gamma_1 \cup \Gamma_2 Ω=Γ1Γ2, q = ∂ u n ⃗ q = \frac{\partial u}{\vec{n}} q=n u n ⃗ \vec{n} n 为边界上的外法向量。

对方程 (1) 两端各乘一个具有二阶导数的任一连续函数 v ( x , y ) v(x, y) v(x,y) 并在 Ω \Omega Ω 上积分得到如下方程:
∫ Ω ( ∇ 2 u − f ) v d Ω = 0                 ( 4 ) \int_{\Omega}(\nabla^2 u - f ) v d\Omega = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4) Ω(2uf)vdΩ=0               (4)

对等式 (4) 运用分步积分公式可得:

− ∫ Ω ∇ u ⋅ ∇ v d x d y + ∫ ∂ Ω ( ∇ u ⋅ n ⃗ ) v d s − ∫ Ω ( f v ) d x d y = 0                 ( 5 ) -\int_{\Omega} \nabla u \cdot \nabla v dxdy + \int_{\partial \Omega} (\nabla u \cdot \vec{n})v ds - \int_{\Omega} (fv) dxdy = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (5) Ωuvdxdy+Ω(un )vdsΩ(fv)dxdy=0               (5)

对等式 (5) 再次运用分步积分公式可得:

∫ Ω ( ∇ 2 v u − f v ) d x d y + ∫ ∂ Ω ( ∇ u ⋅ n ⃗ ) v d s − ∫ ∂ Ω ( ∇ v ⋅ n ⃗ ) u d s = 0                 ( 6 ) \int_{\Omega} (\nabla^2 v u - fv) dxdy + \int_{\partial \Omega} (\nabla u \cdot \vec{n})v ds - \int_{\partial \Omega} (\nabla v \cdot \vec{n})u ds = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6) Ω(2vufv)dxdy+Ω(un )vdsΩ(vn )uds=0               (6)

由于等式 (4) 和 等式 (6) 是等价的, 故可写成如下形式( 将等式 (4) 带入等式 (6) ):

∫ Ω { ( ∇ 2 u ) v − ( ∇ 2 v ) u } d Ω = ∫ ∂ Ω ( ∇ u ⋅ n ⃗ ) v d s − ∫ ∂ Ω ( ∇ v ⋅ n ⃗ ) u d s                 ( 7 ) \int_{\Omega} \{(\nabla^2 u ) v -(\nabla^2 v ) u \}d\Omega = \int_{\partial \Omega} (\nabla u \cdot \vec{n})v ds - \int_{\partial \Omega} (\nabla v \cdot \vec{n})u ds \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (7) Ω{(2u)v(2v)u}dΩ=Ω(un )vdsΩ(vn )uds               (7)

等式 (7) 也是常见的 Green 第二公式

Green 第二公式也常作为建立BEM(边界元)边界积分方程的出发点。

将等式 (2) 和 (3) 带入等式 (6),则等式 (6) 可写为:

∫ Ω { ( ∇ 2 v ) u − f v } d x d y + ∫ Γ 1 q v d s + ∫ Γ 2 q 0 v d s − ∫ Γ 1 u 0 ∂ v ∂ n ⃗ d s − ∫ Γ 2 u ∂ v ∂ n ⃗ d s = 0                 ( 8 ) \begin{array}{c} \int_{\Omega}\left\{\left(\nabla^{2} v\right) u-f v\right\} dxdy+\int_{\Gamma_{1}} q vds +\int_{\Gamma_{2}} q_0 v ds - \\ \\ \int_{\Gamma_{1}} u_0 \frac{\partial v}{\partial \vec{n}} ds-\int_{\Gamma_{2}} u \frac{\partial v}{\partial \vec{n}} ds=0 \end{array} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (8) Ω{(2v)ufv}dxdy+Γ1qvds+Γ2q0vdsΓ1u0n vdsΓ2un vds=0               (8)

为了恢复 L a p a c e Lapace Lapace 算子项 ∇ 2 u \nabla^2u 2u,对等式 (8) 进行逆向分部积分有:
∫ Ω ( − ∂ v ∂ x ∂ u ∂ x − ∂ v ∂ y ∂ u ∂ y − f v ) d x d y + ∫ Γ ∂ v ∂ n ⃗ u d s + ∫ Γ 1 q v d s + ∫ Γ 2 q 0 v d s − ∫ Γ 1 u 0 ∂ v ∂ n ⃗ d s − ∫ Γ 2 u ∂ v ∂ n ⃗ d s = 0                 ( 9 ) \begin{array}{c} \int_{\Omega}\left(-\frac{\partial v}{\partial x} \frac{\partial u}{\partial x}-\frac{\partial v}{\partial y} \frac{\partial u}{\partial y}-f v\right) dxdy+\int_{\Gamma} \frac{\partial v}{\partial \vec{n}} u ds+\int_{\Gamma_{1}} q v ds+ \\ \\ \int_{\Gamma_{2}} q_0 v ds-\int_{\Gamma_{1}} u_0 \frac{\partial v}{\partial \vec{n}} ds - \int_{\Gamma_{2}} u \frac{\partial v}{\partial \vec{n}} ds = 0 \end{array} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (9) Ω(xvxuyvyufv)dxdy+Γn vuds+Γ1qvds+Γ2q0vdsΓ1u0n vdsΓ2un vds=0               (9)


∫ Ω ( − ∂ v ∂ x ∂ u ∂ x − ∂ v ∂ y ∂ u ∂ y − f v ) d x d y + ∫ Γ 1 ∂ v ∂ n ⃗ u d s + ∫ Γ 1 q v d s + ∫ Γ 2 q 0 v d s − ∫ Γ 1 u 0 ∂ v ∂ n ⃗ d s = 0                 ( 10 ) \begin{array}{c} \int_{\Omega}\left(-\frac{\partial v}{\partial x} \frac{\partial u}{\partial x}-\frac{\partial v}{\partial y} \frac{\partial u}{\partial y}-f v\right) dxdy + \int_{\Gamma_{1}} \frac{\partial v}{\partial \vec{n}} u ds+\int_{\Gamma_{1}} q v ds+ \\ \\ \int_{\Gamma_{2}} q_0v ds - \int_{\Gamma_{1}} u_0 \frac{\partial v}{\partial \vec{n}} ds=0 \end{array} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (10) Ω(xvxuyvyufv)dxdy+Γ1n vuds+Γ1qvds+Γ2q0vdsΓ1u0n vds=0               (10)

对等式 (10) 再次分部积分有:
∫ Ω { ( ∇ 2 u ) v − f v } d x d y − ∫ Γ v q d s + ∫ Γ 1 ∂ v ∂ n ⃗ u d s + ∫ Γ 1 q v d s + ∫ Γ 2 q 0 v d s − ∫ Γ 1 u 0 ∂ v ∂ n ⃗ d s = 0                 ( 11 ) \begin{array}{c} \int_{\Omega}\left\{\left(\nabla^{2} u\right) v - f v\right\} dxdy-\int_{\Gamma} v q ds+\int_{\Gamma_{1}} \frac{\partial v}{\partial \vec{n}} u ds+\int_{\Gamma_{1}} q v ds+ \\ \\ \int_{\Gamma_{2}} q_0 v ds - \int_{\Gamma_{1}} u_0 \frac{\partial v}{\partial \vec{n}} ds=0 \end{array} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (11) Ω{(2u)vfv}dxdyΓvqds+Γ1n vuds+Γ1qvds+Γ2q0vdsΓ1u0n vds=0               (11)

等式 (11) 中 Γ 1 \Gamma_1 Γ1 的边界积分项被抵消后得:

∫ Ω { ( ∇ 2 u ) v − f v } d x d y − ∫ Γ 2 v q d s + ∫ Γ 1 ∂ v ∂ n ⃗ u d s + ∫ Γ 2 q 0 v d s − ∫ Γ 1 u 0 ∂ v ∂ n ⃗ d s = 0                 ( 12 ) \begin{array}{c} \int_{\Omega}\left\{\left(\nabla^{2} u\right) v - f v\right\} dxdy-\int_{\Gamma_{2}} v q ds+\int_{\Gamma_{1}} \frac{\partial v}{\partial \vec{n}} u ds+ \\ \\ \int_{\Gamma_{2}} q_0 v ds - \int_{\Gamma_{1}} u_0 \frac{\partial v}{\partial \vec{n}} ds=0 \end{array} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (12) Ω{(2u)vfv}dxdyΓ2vqds+Γ1n vuds+Γ2q0vdsΓ1u0n vds=0               (12)

经整理得:

∫ Ω { ( ∇ 2 u ) v − f v } d x d y − ∫ Γ 2 ( q − q 0 ) v d s + ∫ Γ 1 ( u − u 0 ) ∂ v ∂ n ⃗ d s = 0                 ( 13 ) \begin{array}{c} \int_{\Omega}\left\{\left(\nabla^{2} u\right) v - f v\right\} dxdy-\int_{\Gamma_{2}}(q-q_0) v ds + \\ \\ \int_{\Gamma_{1}}(u-u_0) \frac{\partial v}{\partial \vec{n}} ds=0 \end{array} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (13) Ω{(2u)vfv}dxdyΓ2(qq0)vds+Γ1(uu0)n vds=0               (13)

(13) 式表明了基本函数 u u u 在定义域内要满足微分方程, 同时还要满足基本及自然边界条件。

从Green定理出发可以强调边界元法的严密数学基础, 可见上面的等式 (13) 。但是边界积分方程和边界元法是两回事, 边界元法仍然是近似解法,Green定理也可以从分部积分得到, 如等式 (7)。

直接边界积分方程

以上节介绍的 G r e e n Green Green 第二公式 (7) 为出发点,

∫ Ω { ( ∇ 2 u ) v − ( ∇ 2 v ) u } d Ω = ∫ ∂ Ω ( ∇ u ⋅ n ⃗ ) v d s − ∫ ∂ Ω ( ∇ v ⋅ n ⃗ ) u d s                 ( 14 ) \int_{\Omega} \{(\nabla^2 u ) v -(\nabla^2 v ) u \}d\Omega = \int_{\partial \Omega} (\nabla u \cdot \vec{n})v ds - \int_{\partial \Omega} (\nabla v \cdot \vec{n})u ds \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (14) Ω{(2u)v(2v)u}dΩ=Ω(un )vdsΩ(vn )uds               (14)

假设 Ω \Omega Ω 为一个封闭的三维均匀介质区域, Γ = ∂ Ω \Gamma = \partial \Omega Γ=Ω 为其边界, u u u 代表该区域的电势分布, ∂ u n ⃗ \frac{\partial u}{\vec{n}} n u 为边界面上法向电场强度(边界法向定义为向外),v 为另一个一、 二阶偏导连续的标量函数。

取三维拉普拉斯方程基本解 u ∗ ( r ) = 1 4 π r u^*(r) = \frac{1}{4\pi r} u(r)=4πr1 代替公式 (14) 中的 v,又根据均匀介质静电场中 u u u 满足拉普拉斯方程(即 ∇ 2 u = 0 \nabla^2 u = 0 2u=0 代入等式 14),可得如下边界积分方程:

∫ Γ ( u ∂ u ∗ n ⃗ − u ∗ ∂ u n ⃗ ) d s = ∫ Ω u ( ∇ 2 u ∗ ) d Ω = − c s u s                 ( 15 ) \int_{\Gamma} (u \frac{\partial u^*}{\vec{n}} - u^* \frac{\partial u}{\vec{n}})ds = \int_{\Omega} u (\nabla^2 u^* ) d\Omega = -c_s u_s \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (15) Γ(un uun u)ds=Ωu(2u)dΩ=csus               (15)

式中,带 ∗ * 变量为源点 s s s 处基本解对应的场分布函数, u s u_s us 是源点 s s s 的电势。若源点在区域内,系数 c s = 1 c_s = 1 cs=1; 若源点在区域边界上,系数 c s = 0.5 c_s= 0.5 cs=0.5 。令 q = ∂ u n ⃗ q = \frac{\partial u}{\vec{n}} q=n u ,即面上的电场强度, q ∗ = ∂ u ∗ n ⃗ = − 1 4 π r 2 ⋅ ∂ r n ⃗ = − ( r ⃗ n ⃗ ) 4 π r 3 q^* = \frac{\partial u^*}{\vec{n}} = - \frac{1}{4\pi r^2} \cdot \frac{\partial r}{\vec{n}} = - \frac{(\vec{r} \vec{n})}{4\pi r^3} q=n u=4πr21n r=4πr3(r n ) 式 (15) 变为:

∫ Γ u q ∗ − u ∗ q d s = − c s u s                 ( 16 ) \int_{\Gamma} uq^* - u^*qds = -c_s u_s \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (16) Γuquqds=csus               (16)

式 (16) 被称为直接边界积分方程。(16) 式中所含的未知量 u u u 的物理意义是电势, 未知量 q q q 的物理意义是电场强度。从中可以看出,若已知边界上的物理量 u , q u,q u,q, 就可以积分算得域内及边界上任一评价点的电势(只需将基本解的源点定在该 评价点上)。反过来,若边界上场量未知,也可通过设置一系列源点构造一组方程 (16) 并求解。这种基于直接边界积分方程的离散化求解方法称为直接边界元法。

直接边界积分方程离散化

任选一介质区域 i i i,直接边界积分方程 (16) 可写为

c s u s ( i ) + ∫ Γ i q s ∗ u ( i ) d s = ∫ Γ i u s ∗ q ( i ) d s ,          ( i = 1 , … , M )         ( 17 ) c_{s} u_{s}^{(i)}+\int_{\Gamma_{i}} q_{s}^{*} u^{(i)} ds=\int_{\Gamma_{i}} u_{s}^{*} q^{(i)} ds, \ \ \ \ \ \ \ \ (i=1, \ldots, M) \ \ \ \ \ \ \ (17) csus(i)+Γiqsu(i)ds=Γiusq(i)ds,        (i=1,,M)       (17)

其中 u s ( i ) u_s^{(i)} us(i) 是介质 i i i 中源点 s s s 的电势, c s c_s cs 是与源点附近边界几何形状有关的常数。 u s ∗ = 1 / 4 π r u_s^{*} = 1/4 \pi r us=1/4πr L a p a c e Lapace Lapace 方程基本解, 其沿单位外法向 n ⃗ \vec{n} n 的方向导数 q s ∗ = ∂ u s ∗ / ∂ n ⃗ = − ( r ⃗ , n ⃗ ) / 4 π r 3 q_s^{*} = \partial u_s^{*} / \partial \vec{n} = -(\vec{r}, \vec{n}) / 4 \pi r^3 qs=us/n =(r ,n )/4πr3, r ⃗ \vec{r} r 为源点到被积分场点的向量( r r r 为相应的欧式距离)。 Γ i \Gamma_i Γi 是包围介质区域 i i i 的边界, 设共有 M M M 个介质区域。变量的上标 ( i ) (i) (i) 为其所属介质区域 i i i 的标识。

为了求解方程 (17),须离散化区域边界,即将介质 i ( i = 1 , . . . , M ) i (i=1, ..., M) i(i=1,...,M) 的边界 Γ i \Gamma_i Γi 离散为 N i N_i Ni 个边界元(介质交界面上的元同时属于相邻两种介质)。下用常数元推导,每个元的中心点为变量结点,同时取其为源点,可得一系列离散化边界积分方程 :

c k u k ( i ) + ∑ j = 1 N i ( ∫ Γ j ( i ) q k , ( i ) ∗ d s ) u j ( i ) = ∑ j = 1 N i ( ∫ Γ j ( i ) u k , ( i ) ∗ d s ) q j ( i ) , ( k = 1 , … , N i ,    i = 1 , … , M )           ( 18 ) \begin{array}{c} c_{k} u_{k}^{(i)}+\sum_{j=1}^{N_{i}}\left(\int_{\Gamma_{j}^{(i)}} q_{k,(i)}^{*} ds\right) u_{j}^{(i)}=\sum_{j=1}^{N_{i}}\left(\int_{\Gamma_{j}^{(i)}} u_{k,(i)}^{*} ds\right) q_{j}^{(i)}, \\ \\ \left(k=1, \ldots, N_{i},\ \ i=1, \ldots, M\right) \end{array} \ \ \ \ \ \ \ \ \ (18) ckuk(i)+j=1Ni(Γj(i)qk,(i)ds)uj(i)=j=1Ni(Γj(i)uk,(i)ds)qj(i),(k=1,,Ni,  i=1,,M)         (18)

其中, Γ j ( i ) \Gamma_j^{(i)} Γj(i) Γ i \Gamma_i Γi 上第 j j j 个边界元。 u j ( i ) u_j^{(i)} uj(i) q j ( i ) q_j^{(i)} qj(i) 分别表示 Γ j ( i ) \Gamma_j^{(i)} Γj(i) 上的电势与法向电场强度, c k = 0.5 c_k = 0.5 ck=0.5。这种逐一取变量结点为源点构造离散方程组的方法常被称为“点配置”(collocation)法。

在边界元上进行积分后, 可得如下线性方程组:

1 2 u k ( i ) + ∑ j = 1 N i H k j ( i ) u j ( i ) = ∑ j = 1 N i G k j ( i ) q j ( i ) ,    ( k = 1 , … , N i ,    i = 1 , … , M )           ( 19 ) \frac{1}{2}u_k^{(i)} + \sum_{j=1}^{N_i} H_{kj}^{(i)} u_j^{(i)} = \sum_{j=1}^{N_i} G_{kj}^{(i)} q_j^{(i)}, \ \ \left(k=1, \ldots, N_{i},\ \ i=1, \ldots, M\right) \ \ \ \ \ \ \ \ \ (19) 21uk(i)+j=1NiHkj(i)uj(i)=j=1NiGkj(i)qj(i),  (k=1,,Ni,  i=1,,M)         (19)

其中 H k j ( i ) H_{kj}^{(i)} Hkj(i) G k j ( i ) G_{kj}^{(i)} Gkj(i) 分别是在边界元 Γ j ( i ) \Gamma_j^{(i)} Γj(i) 作积分得到:

H k j ( i ) = ∫ Γ j ( i ) q k , ( i ) ∗ d s ,       G k j ( i ) = ∫ Γ j ( i ) u k , ( i ) ∗ d s           ( 20 ) H_{kj}^{(i)} = \int_{\Gamma_j^{(i)}} q_{k,(i)}^{*} ds,\ \ \ \ \ G_{kj}^{(i)} = \int_{\Gamma_j^{(i)}} u_{k,(i)}^{*} ds \ \ \ \ \ \ \ \ \ (20) Hkj(i)=Γj(i)qk,(i)ds,     Gkj(i)=Γj(i)uk,(i)ds         (20)

公式 (20) 的边界元积分可分为两类,一类是源点落在被积元上的奇异积分。奇异积分数量不多,一般采用极坐标下的解析积分进行计算。另一类是源点不落在被积元上的非奇异积分。非奇异积分数量很多,计算量很大。当源点到被积元的距离与元的大小相比很小时,称为近奇异积分,是边界元积分的难点。高斯(Gauss) 积分是计算非奇异积分的常用方法。

间接边界积分方程

间接边界积分方程离散化

直接边界元与间接边界元对比

将拉普拉斯方程转化为边界积分方程有不同方法,常用的有以下两种。第一种方法是利用格林(Green)公式和拉普拉斯方程基本解的性质导出直接边界积方程并离散化。之所以称为直接边界积分方程,是因为积分方程所含未知量是与拉普拉斯方程待求函数直接相关的边界量,这些边界量往往具有明确的物理意义。离散直接边界积分方程的方法称为直接边界元素法。第二种方法是通过位势理论,即在边界上引入某个虚拟的密度函数去构造边界积分方程,这样得到一个关于密度函数的积分方程,称为间接边界积分方程。这种以间接边界积分方程为离散化对象的方法称为间接边界元素法。所谓间接是因为引入的密度函数并不是拉普拉斯方程所直接关心的未知函数,应当指出,与间接边界元素法相比,直接边界元素法能够十分方便地处理带混合边界条件的拉普拉斯方程。

  • 13
    点赞
  • 44
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
边界元法(Boundary Element Method, BEM)是一种常用于求解边界值问题的数值方法边界元法的核心思想是将求解区域分割为边界和内部两种类型。边界是指在问题的边界上离散取点,而内部是指在边界内部取点。然后,通过将问题转换为边界积分方程,利用边界与内部的构造和边界条件的约束,使用离散化的方法求解有限个界面积分方程,从而得到问题的数值解。 Matlab是一种功能强大的科学计算软件,广泛应用于工程与科学计算领域。在Matlab中,边界元法可以通过编写相应的程序进行实现。 首先,需要在Matlab中定义边界的构造和内部的取点。边界的构造通常利用基函数来实现,可以选择常用的线性或非线性基函数。内部的取点可以采用不同的方法,如均匀取点、Gauss-Legendre积分点等,以保证数值计算的精度。 然后,需要将问题转化为边界积分方程表示。根据具体问题的边界条件,可以得到相应的边界积分方程。在编写程序时,需要将边界积分方程离散化,将迭代求解问题转化为求解一个线性方程组。 最后,利用Matlab中的矩阵运算和求解线性方程组的函数,可以求解得到问题的数值解。根据具体问题的要求,可以通过调整边界和内部的离散点数量、改变基函数的选择和内部的取点方式等来控制求解的精度和计算效率。 总之,边界元法是一种常用且有效的数值方法,Matlab作为一种强大的科学计算软件,可以提供丰富的工具和函数来实现边界元法的求解过程。
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值