文章目录
前言
边界元法是在经典积分方程法和有限元方法基础上发展起来的一种偏微分方程数值解法。 Brebbia 于 1978 年用加权余量法统一边界积分方程的不同推导,并首次以“边界元法”这一名称统一称呼此前边界积分方程离散化的多种不同方法。
边界元素法的基本思路是:将计算区域内的微分方程变换为区域边界上的积分方程;然后将区域边界剖分为有限大小的边界元素,从而将边界积分方程离散化为线性代数方程组,并经求解线性方程组而得到边界节点物理量。
边界元素法与有限差分和有限元素法相比,具有以下显著优点:
- 计算精度高:边界元素法直接求解边界广义场源的分布。边界上任一点场量将通过线性叠加各离散广义场源的作用而得到,不必经微分计算间接得到。
- 计算速度快:由于仅对边界离散化,它降低了问题的维数,不必如有限差分及有限元法那样需引进大量区域变量,使变量数极大减少。
- 处理复杂边界能力强:由于只离散边界,大大降低了几何处理复杂性,从而增强了处理复杂边界的能力。
不过,由于边界元计算形成的线性方程组系数矩阵一般是稠密、不对称的,随着集成密度指数级增加,问题规模不断扩大,给方程求解带来越来越大的困难,特别是对三维问题。
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(∂x∂Q−∂y∂P)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[(∂x∂u)2+(∂y∂u)2]dxdy=−∬suΔudxdy+∮cu∂n∂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=∂x∂ucosθ+∂y∂usinθ
可知有
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)2dy=(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}
∮cu∂n∂uds=∮cu(∂x∂u(dx)2+(dy)2dy+∂y∂u(dx)2+(dy)2dy)(dx)2+(dy)2=∮cu∂x∂udy−u∂y∂udx
由 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(∂x∂Q−∂y∂P)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}
∮cu∂x∂udy−u∂y∂udx=∬S[∂x∂(u∂x∂u)−∂y∂(−u∂y∂u)]dxdy=∬S[∂x∂(u∂x∂u)+∂y∂(u∂y∂u)]dxdy=∬S[∂x∂(∂x∂u)u+(∂x∂u)2+∂y∂(∂y∂u)u+(∂y∂u)2]dxdy=∬S[(∂x∂u)2+(∂y∂u)2]dxdy+∬S[∂x∂(∂x∂u)u+∂y∂(∂y∂u)u]dxdy=∬S[(∂x∂u)2+(∂y∂u)2]dxdy+∬Su[∂x∂(∂x∂u)+∂y∂(∂y∂u)]dxdy=∬S[(∂x∂u)2+(∂y∂u)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}
∮cu∂n∂uds=∬S[(∂x∂u)2+(∂y∂u)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∂uv∂n∂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Δu−uΔ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∂u∂n∂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∂uv−∂n∂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}}
=∮Cv∂x∂u(dx)2+(dy)2dy−v∂y∂u(dx)2+(dy)2dx
−
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}}
−u∂x∂v(dx)2+(dy)2dy
+
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}
+u∂y∂v(dx)2+(dy)2dx(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
=∮Cv∂x∂udy−v∂y∂udx−u∂x∂vdy+u∂y∂vdx
=
∮
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(u∂y∂v−v∂y∂u)dx+(v∂x∂u−u∂x∂v)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(∂x∂Q−∂y∂P)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=(u∂y∂v−v∂y∂u)Q=(v∂x∂u−u∂x∂v)∂x∂Q=∂x∂(v∂x∂u−u∂x∂v)=∂x∂v∂x∂u+v∂x2∂2u−∂x∂v∂x∂u−u∂x2∂2v=v∂x2∂2u−u∂x2∂2v
同理
∂ 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}} ∂y∂P=u∂y2∂2v−v∂y2∂2u
故有:
∬ 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(∂x∂Q−∂y∂P)dxdy=∬D(v∂x2∂2u−u∂x2∂2v−u∂y2∂2v+v∂y2∂2u)dxdy=∬DvΔu−uΔ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π1∮C(u∂n∂lnr−lnr∂n∂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(u∂n∂lnr−lnr∂n∂u)ds=∬DvΔu−uΔ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(u∂n∂lnr−lnr∂n∂u)ds=∬DvΔu−uΔ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(u∂n∂lnr−lnr∂n∂u)ds=∮∂D(u∂n∂lnr−lnr∂n∂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}
∮∂Dlnr∂n∂uds=lnt∮∂D∂n∂uds=lnt∬DΔuds=0∮∂Du∂n∂lnrds=∮∂Du∂r∂lnrds=∮∂Dut1ds=t1∮∂Duds=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)
t→0lim∮C(u∂n∂lnr−lnr∂n∂u)ds=t→0lim2π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π1∮C(u∂n∂lnr−lnr∂n∂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 ∫Ω∇⋅(c∇u)vdxdy=∫∂Ω(c∇u⋅n)vds−∫Ωc∇u⋅∇vdxdy
∬ Ω 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 ∬Ωf∂x∂gdxdy=∮∂Ω(fg)dy−∬Ωg∂x∂fdxdy
∬ Ω 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 ∬Ωf∂y∂gdxdy=−∮∂Ω(fg)dx−∬Ωg∂y∂fdxdy
证明:(由Green公式) ∬ Ω ∂ F ∂ x d x d y = ∮ ∂ Ω F d y \iint_{\Omega} \frac{\partial F}{\partial x} dxdy = \oint_{\partial \Omega} Fdy ∬Ω∂x∂Fdxdy=∮∂Ω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 ∬Ω(f∂x∂g+g∂x∂f)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 ∬Ωf∂x∂gdxdy=∮∂Ω(fg)dy−∬Ωg∂x∂fdxdy
注: ∫ u v ′ d x = u v − ∫ v u ′ d x \int u v'dx = uv - \int vu'dx ∫uv′dx=uv−∫vu′dx
积分方程
以二维的 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)
∫Ω(∇2u−f)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) −∫Ω∇u⋅∇vdxdy+∫∂Ω(∇u⋅n)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) ∫Ω(∇2vu−fv)dxdy+∫∂Ω(∇u⋅n)vds−∫∂Ω(∇v⋅n)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Ω=∫∂Ω(∇u⋅n)vds−∫∂Ω(∇v⋅n)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)u−fv}dxdy+∫Γ1qvds+∫Γ2q0vds−∫Γ1u0∂n∂vds−∫Γ2u∂n∂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)
∫Ω(−∂x∂v∂x∂u−∂y∂v∂y∂u−fv)dxdy+∫Γ∂n∂vuds+∫Γ1qvds+∫Γ2q0vds−∫Γ1u0∂n∂vds−∫Γ2u∂n∂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)
∫Ω(−∂x∂v∂x∂u−∂y∂v∂y∂u−fv)dxdy+∫Γ1∂n∂vuds+∫Γ1qvds+∫Γ2q0vds−∫Γ1u0∂n∂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)v−fv}dxdy−∫Γvqds+∫Γ1∂n∂vuds+∫Γ1qvds+∫Γ2q0vds−∫Γ1u0∂n∂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)v−fv}dxdy−∫Γ2vqds+∫Γ1∂n∂vuds+∫Γ2q0vds−∫Γ1u0∂n∂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)v−fv}dxdy−∫Γ2(q−q0)vds+∫Γ1(u−u0)∂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Ω=∫∂Ω(∇u⋅n)vds−∫∂Ω(∇v⋅n)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∂u∗−u∗n∂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πr21⋅n∂r=−4πr3(rn) 式 (15) 变为:
∫ Γ u q ∗ − u ∗ q d s = − c s u s ( 16 ) \int_{\Gamma} uq^* - u^*qds = -c_s u_s \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (16) ∫Γuq∗−u∗qds=−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)+∫Γiqs∗u(i)ds=∫Γius∗q(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=1∑NiHkj(i)uj(i)=j=1∑NiGkj(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)公式和拉普拉斯方程基本解的性质导出直接边界积方程并离散化。之所以称为直接边界积分方程,是因为积分方程所含未知量是与拉普拉斯方程待求函数直接相关的边界量,这些边界量往往具有明确的物理意义。离散直接边界积分方程的方法称为直接边界元素法。第二种方法是通过位势理论,即在边界上引入某个虚拟的密度函数去构造边界积分方程,这样得到一个关于密度函数的积分方程,称为间接边界积分方程。这种以间接边界积分方程为离散化对象的方法称为间接边界元素法。所谓间接是因为引入的密度函数并不是拉普拉斯方程所直接关心的未知函数,应当指出,与间接边界元素法相比,直接边界元素法能够十分方便地处理带混合边界条件的拉普拉斯方程。