电磁场仿真-主题078-FEM-BEM耦合方法

kkchenjj

92人浏览 · 2026-03-25 21:26:53

kkchenjj · 2026-03-25 21:26:53 发布

主题078：FEM-BEM耦合方法

1. 引言

1.1 开域电磁场问题的挑战

在电磁场数值仿真中，开域问题（Open Domain Problems）是一类具有重要工程意义但求解困难的课题。这类问题的特点是计算区域延伸至无穷远，例如：

天线辐射问题：电磁波向自由空间无限传播
电磁散射问题：目标体对入射波的散射场在远场衰减
电磁兼容分析：设备辐射干扰的传播范围无法预先限定

传统的纯有限元法(FEM)在处理开域问题时面临根本困难：必须在有限计算域内截断无限区域，这引入了人工边界。为消除人工边界的反射效应，需要采用复杂的吸收边界条件(ABC)或完美匹配层(PML)技术，增加了计算复杂度和误差来源。

1.2 FEM-BEM耦合的优势

有限元-边界元耦合方法(FEM-BEM Coupling)巧妙地结合了两种方法的优势：

有限元法(FEM)的优势：

擅长处理非均匀、各向异性介质
能够精确模拟复杂几何结构
稀疏矩阵特性，内存效率高

边界元法(BEM)的优势：

自动满足无穷远辐射条件
只需离散边界，降低维度
适合处理无限域问题

FEM-BEM耦合的核心思想：
将计算域分为两部分：

内部区域（有界）：使用FEM处理复杂介质和结构
外部区域（无限）：使用BEM自动满足远场条件

两部分在耦合界面上通过连续性条件连接，形成统一的求解系统。
在这里插入图片描述

1.3 本教程目标

本教程将系统介绍FEM-BEM耦合方法的理论基础、数值实现和工程应用，包括：

FEM和BEM的基本原理回顾
耦合方法的理论推导
Python代码实现与数值实验
电磁散射和天线辐射问题的仿真案例
方法性能分析与比较

2. FEM-BEM耦合的基本原理

2.1 问题域的分解

考虑二维Helmholtz方程描述的开域电磁问题：

$Ω\nabla^2 u + k^2 u = f \quad \text{in } \Omega$

其中 $Ω\Omega$ 是无限域。FEM-BEM方法将其分解为：

内部域 $Ωi\Omega_i$ （有界）：
$Ωi\nabla^2 u_i + k_i^2 u_i = f \quad \text{in } \Omega_i$

外部域 $Ωe\Omega_e$ （无限）：
$Ωe\nabla^2 u_e + k_0^2 u_e = 0 \quad \text{in } \Omega_e$

其中 $Γ=∂Ωi=∂Ωe\Gamma = \partial\Omega_i = \partial\Omega_e$ 是耦合界面。

2.2 耦合条件

在界面 $Γ\Gamma$ 上，必须满足两个连续性条件：

场连续性（Dirichlet条件）：
$Γu_i = u_e \quad \text{on } \Gamma$

通量连续性（Neumann条件）：
$Γ\frac{\partial u_i}{\partial n} = \frac{\partial u_e}{\partial n} \quad \text{on } \Gamma$

这两个条件确保物理量在界面上的连续过渡。

2.3 变分 formulation

内部域采用FEM的弱形式：

$∫Ωi(∇ui⋅∇v−ki2uiv)dΩ−∫Γ∂ui∂nvdΓ=∫ΩifvdΩ\int_{\Omega_i} \left( \nabla u_i \cdot \nabla v - k_i^2 u_i v \right) d\Omega - \int_{\Gamma} \frac{\partial u_i}{\partial n} v d\Gamma = \int_{\Omega_i} f v d\Omega$

外部域采用BEM的边界积分方程：

$c(x)ue(x)+∫Γ∂G∂n(x,y)ue(y)dΓ(y)=∫ΓG(x,y)∂ue∂n(y)dΓ(y)c(x)u_e(x) + \int_{\Gamma} \frac{\partial G}{\partial n}(x,y) u_e(y) d\Gamma(y) = \int_{\Gamma} G(x,y) \frac{\partial u_e}{\partial n}(y) d\Gamma(y)$

其中 $G$ 是Helmholtz格林函数。

2.4 耦合系统的矩阵形式

离散化后，耦合系统可写为块矩阵形式：

$[AFEMCFICIFABEM][uFEMuBEM]=[fFEMfBEM]\begin{bmatrix} A_{FEM} & C_{FI} \\ C_{IF} & A_{BEM} \end{bmatrix} \begin{bmatrix} u_{FEM} \\ u_{BEM} \end{bmatrix} = \begin{bmatrix} f_{FEM} \\ f_{BEM} \end{bmatrix}$

其中：

$A_{FEM}$ ：FEM刚度矩阵（稀疏）
$A_{BEM}$ ：BEM系统矩阵（稠密）
$C_{FI}, C_{IF}$ ：耦合矩阵

3. 有限元法(FEM)基础

3.1 二维FEM求解器实现

class FEM2D:
    """二维有限元求解器 - 处理内部区域"""
    
    def __init__(self, nodes, elements, epsilon_r=1.0, mu_r=1.0):
        """
        参数:
            nodes: 节点坐标数组 (N, 2)
            elements: 三角形单元连接数组 (M, 3)
            epsilon_r: 相对介电常数
            mu_r: 相对磁导率
        """
        self.nodes = nodes
        self.elements = elements
        self.n_nodes = len(nodes)
        self.n_elements = len(elements)
        self.epsilon_r = epsilon_r
        self.mu_r = mu_r
        self.k = 2 * np.pi  # 波数

3.2 单元矩阵计算

对于三角形线性单元，单元刚度矩阵和质量矩阵的计算如下：

def compute_element_matrices(self, element_idx):
    """计算单元刚度矩阵和质量矩阵"""
    elem = self.elements[element_idx]
    x = self.nodes[elem, 0]
    y = self.nodes[elem, 1]
    
    # 计算单元面积
    area = 0.5 * abs((x[1] - x[0]) * (y[2] - y[0]) - 
                     (x[2] - x[0]) * (y[1] - y[0]))
    
    # 形状函数梯度
    b = np.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]])
    c = np.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]])
    
    # 刚度矩阵 (梯度项)
    K_elem = np.outer(b, b) + np.outer(c, c)
    K_elem = K_elem / (4 * area)
    
    # 质量矩阵 (势能项)
    M_elem = area / 12 * np.array([[2, 1, 1],
                                    [1, 2, 1],
                                    [1, 1, 2]])
    
    return K_elem, M_elem, area

理论说明：

对于线性三角形单元，形状函数为：
$Ni=ai+bix+ciy,i=1,2,3N_i = a_i + b_i x + c_i y, \quad i = 1, 2, 3$

其中系数由节点坐标确定：
$b1=y2−y3,c1=x3−x2b_1 = y_2 - y_3, \quad c_1 = x_3 - x_2$

单元刚度矩阵：
$Kij(e)=∫Ωe∇Ni⋅∇NjdΩ=14Ae(bibj+cicj)K_{ij}^{(e)} = \int_{\Omega_e} \nabla N_i \cdot \nabla N_j d\Omega = \frac{1}{4A_e}(b_i b_j + c_i c_j)$

单元质量矩阵（一致质量矩阵）：
$Mij(e)=∫ΩeNiNjdΩ=Ae12(1+δij)M_{ij}^{(e)} = \int_{\Omega_e} N_i N_j d\Omega = \frac{A_e}{12}(1 + \delta_{ij})$

3.3 全局矩阵组装

def assemble_system(self, k=None):
    """组装全局系统矩阵"""
    if k is None:
        k = self.k
    
    # 全局矩阵
    rows, cols, data_K, data_M = [], [], [], []
    
    for e in range(self.n_elements):
        K_elem, M_elem, _ = self.compute_element_matrices(e)
        elem = self.elements[e]
        
        for i in range(3):
            for j in range(3):
                rows.append(elem[i])
                cols.append(elem[j])
                data_K.append(K_elem[i, j])
                data_M.append(M_elem[i, j])
    
    K_global = sparse.coo_matrix((data_K, (rows, cols)), 
                                  shape=(self.n_nodes, self.n_nodes)).tocsr()
    M_global = sparse.coo_matrix((data_M, (rows, cols)), 
                                  shape=(self.n_nodes, self.n_nodes)).tocsr()
    
    # Helmholtz方程: -∇²u - k²u = f
    A = K_global - k**2 * M_global
    
    return A, K_global, M_global

3.4 边界条件处理

def apply_dirichlet_boundary(self, A, b, bc_nodes, bc_values):
    """应用Dirichlet边界条件"""
    A = A.tolil()
    
    for node, value in zip(bc_nodes, bc_values):
        A[node, :] = 0
        A[node, node] = 1
        b[node] = value
    
    return A.tocsr(), b

3.5 FEM求解器演示结果

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了FEM求解器的网格和计算结果：

左图：规则三角形网格，用于离散计算域
右图：Helmholtz方程的数值解实部，显示波动特性

4. 边界元法(BEM)基础

4.1 边界积分方程

对于外部Helmholtz问题，使用边界积分方程：

$\int_{\Gamma} q^*(x,y) u(y) d\Gamma(y) = \int_{\Gamma} u^*(x,y) q(y) d\Gamma(y)$

其中：

$u^*$ ：Helmholtz格林函数（基本解）
$q∗=∂u∗/∂nq^* = \partial u^*/\partial n$ ：格林函数法向导数
$c (x)$ ：几何系数（光滑边界上 $c = 0.5$ ）

4.2 二维Helmholtz格林函数

二维Helmholtz方程的格林函数为：

$\frac{i}{4} H_0^{(1)}(kr)$

其中 $H_0^{(1)}$ 是第一类零阶Hankel函数， $r = ∣ x - y ∣$ 。

def green_function(self, x, y, x0, y0):
    """二维Helmholtz格林函数"""
    r = np.sqrt((x - x0)**2 + (y - y0)**2)
    if r < 1e-10:
        return 0.0
    return 0.25j * hankel1(0, self.k * r)

Hankel函数的性质：

渐近行为： $H0(1)(kr)∼2πkrei(kr−π/4)H_0^{(1)}(kr) \sim \sqrt{\frac{2}{\pi kr}} e^{i(kr - \pi/4)}$ （当 $\to \infty$ ）
奇异性： $H0(1)(kr)∼2iπln⁡(kr)H_0^{(1)}(kr) \sim \frac{2i}{\pi} \ln(kr)$ （当 $\to 0$ ）

4.3 BEM求解器实现

class BEM2D:
    """二维边界元求解器 - 处理无限远场"""
    
    def __init__(self, boundary_nodes, k=2*np.pi):
        self.boundary_nodes = boundary_nodes
        self.n_bound = len(boundary_nodes)
        self.k = k

4.4 边界矩阵计算

def compute_boundary_matrices(self):
    """计算边界元系统矩阵"""
    n = self.n_bound
    H = np.zeros((n, n), dtype=complex)
    G = np.zeros((n, n), dtype=complex)
    
    # 计算边界法向量
    normals = self._compute_normals()
    
    for i in range(n):
        xi, yi = self.boundary_nodes[i]
        
        for j in range(n):
            xj, yj = self.boundary_nodes[j]
            
            if i == j:
                # 对角项处理 (奇异性)
                H[i, j] = 0.5
                G[i, j] = self._compute_singular_integral(i)
            else:
                # 非对角项
                r = np.sqrt((xi - xj)**2 + (yi - yj)**2)
                H[i, j] = -0.25j * self.k * hankel1(1, self.k * r)
                G[i, j] = 0.25j * hankel1(0, self.k * r)
    
    return H, G, normals

4.5 BEM求解器演示结果

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了BEM求解器的特性：

左图：圆形边界离散和法向量分布
右图：BEM系统矩阵 $H$ 的结构，显示稠密矩阵特性

5. 耦合方法的理论推导

5.1 变分 formulation

FEM-BEM耦合的核心是建立统一的变分问题。考虑内部域 $Ωi\Omega_i$ 和边界 $Γ\Gamma$ ：

内部域弱形式：
$∫Ωi(∇u⋅∇v−k2uv)dΩ=∫ΓqvdΓ\int_{\Omega_i} (\nabla u \cdot \nabla v - k^2 u v) d\Omega = \int_{\Gamma} q v d\Gamma$

边界积分方程：
$12u(x)+∫Γq∗udΓ=∫Γu∗qdΓ\frac{1}{2}u(x) + \int_{\Gamma} q^* u d\Gamma = \int_{\Gamma} u^* q d\Gamma$

5.2 耦合系统的建立

将FEM和BEM方程组合，得到耦合系统：

$[K−k2M−LT(12I+H)−G][uq]=[f0]\begin{bmatrix} K - k^2 M & -L^T \\ (\frac{1}{2}I + H) & -G \end{bmatrix} \begin{bmatrix} u \\ q \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}$

其中：

$K, M$ ：FEM刚度矩阵和质量矩阵
$H, G$ ：BEM系统矩阵
$L$ ：耦合矩阵

5.3 对称耦合 formulation

为获得对称系统，可采用Costabel的symmetric coupling方法：

$[AFEM−MT−MABEM][uλ]=[fg]\begin{bmatrix} A_{FEM} & -M^T \\ -M & A_{BEM} \end{bmatrix} \begin{bmatrix} u \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ g \end{bmatrix}$

对称系统的优势：

可使用共轭梯度等迭代求解器
更好的数值稳定性
便于误差分析

6. 二维Helmholtz问题的FEM-BEM耦合

6.1 耦合求解器实现

class FEMBEMCoupling:
    """FEM-BEM耦合求解器"""
    
    def __init__(self, fem_solver, bem_solver, coupling_nodes_fem, 
                 coupling_nodes_bem):
        self.fem = fem_solver
        self.bem = bem_solver
        self.coupling_fem = coupling_nodes_fem
        self.coupling_bem = coupling_nodes_bem
        self.n_coupling = len(coupling_nodes_fem)

6.2 耦合系统构建

def build_coupled_system(self, source_nodes=None, source_values=None):
    """构建耦合系统矩阵"""
    # FEM系统
    A_fem, K_fem, M_fem = self.fem.assemble_system()
    n_fem = self.fem.n_nodes
    
    # BEM系统
    H_bem, G_bem, normals = self.bem.compute_boundary_matrices()
    n_bem = self.bem.n_bound
    
    # 构建耦合矩阵
    n_total = n_fem + n_bem
    A_coupled = np.zeros((n_total, n_total), dtype=complex)
    b_coupled = np.zeros(n_total, dtype=complex)
    
    # FEM块
    A_coupled[:n_fem, :n_fem] = A_fem.toarray()
    
    # BEM块
    A_coupled[n_fem:, n_fem:] = H_bem
    
    # 耦合块
    for i, (fem_idx, bem_idx) in enumerate(zip(self.coupling_fem, 
                                                self.coupling_bem)):
        A_coupled[n_fem + bem_idx, fem_idx] = -1.0
        A_coupled[fem_idx, n_fem + bem_idx] = -1.0
    
    return A_coupled, b_coupled

6.3 系统求解

def solve(self, source_nodes=None, source_values=None):
    """求解耦合系统"""
    A, b = self.build_coupled_system(source_nodes, source_values)
    
    # 求解线性系统
    u = np.linalg.solve(A, b)
    
    # 分离FEM和BEM解
    n_fem = self.fem.n_nodes
    u_fem = u[:n_fem]
    u_bem = u[n_fem:]
    
    return u_fem, u_bem

7. 格林函数与基本解

7.1 Helmholtz格林函数

格林函数是点源在无限域中产生的场，满足：

$∇2G+k2G=−δ(x−x0)\nabla^2 G + k^2 G = -\delta(x - x_0)$

二维格林函数：
$\frac{i}{4} H_0^{(1)}(kr)$

三维格林函数：
$\frac{e^{ikr}}{4\pi r}$

7.2 格林函数的物理意义

格林函数表示位于 $x_0$ 的点源在 $x$ 处产生的响应：

实部：近场振荡特性
虚部：辐射损耗
幅度： $\propto 1/\sqrt{r}$ （二维）或 $1/ r$ （三维）

7.3 格林函数可视化

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了二维Helmholtz格林函数的特性：

左图：实部显示同心圆状的波动模式
中图：虚部表示辐射场的相位特性
右图：幅度随距离衰减，符合 $1/r1/\sqrt{r}$ 规律

8. 耦合界面处理技术

8.1 界面离散化

耦合界面需要同时在FEM和BEM中离散：

FEM侧：界面作为区域边界的一部分
BEM侧：界面作为唯一的离散对象

8.2 节点匹配策略

一对一匹配：
FEM和BEM在界面上使用相同的节点分布，简化耦合矩阵构建。

非匹配网格：
使用Mortar方法或插值技术处理不同离散密度的界面。

8.3 耦合界面可视化

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了FEM-BEM耦合界面的结构：

左图：FEM内部节点（蓝色）和BEM边界节点（红色）的分布
右图：耦合系统矩阵的结构，显示FEM块、BEM块和耦合块的分布

9. 电磁散射问题仿真

9.1 散射问题描述

平面波入射到介质散射体上，产生散射场：

$u_{total} = u_{incident} + u_{scattered}$

9.2 散射问题求解器

class ScatteringProblem:
    """电磁散射问题仿真"""
    
    def __init__(self, frequency=1e9, epsilon_r=4.0, mu_r=1.0):
        self.frequency = frequency
        self.c = 3e8
        self.wavelength = self.c / frequency
        self.k0 = 2 * np.pi / self.wavelength
        self.epsilon_r = epsilon_r
        self.mu_r = mu_r

9.3 圆形散射体网格生成

def create_circular_scatterer(self, radius, n_elements=40):
    """创建圆形散射体网格"""
    theta = np.linspace(0, 2*np.pi, n_elements, endpoint=False)
    boundary_nodes = np.column_stack([
        radius * np.cos(theta),
        radius * np.sin(theta)
    ])
    
    # 内部网格 (同心圆)
    n_rings = 5
    all_nodes = [boundary_nodes]
    
    for ring in range(1, n_rings + 1):
        r = radius * (n_rings - ring) / n_rings
        ring_nodes = np.column_stack([
            r * np.cos(theta),
            r * np.sin(theta)
        ])
        all_nodes.append(ring_nodes)
    
    nodes = np.vstack(all_nodes[::-1])
    
    # 创建三角形单元
    elements = []
    n_per_ring = n_elements
    for ring in range(n_rings):
        for i in range(n_per_ring):
            i1 = ring * n_per_ring + i
            i2 = ring * n_per_ring + (i + 1) % n_per_ring
            i3 = (ring + 1) * n_per_ring + i
            i4 = (ring + 1) * n_per_ring + (i + 1) % n_per_ring
            
            elements.append([i1, i2, i3])
            if ring < n_rings - 1:
                elements.append([i2, i4, i3])
    
    return nodes, np.array(elements), boundary_nodes

9.4 散射场计算结果

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了电磁散射问题的数值解：

左上：网格划分，显示同心圆结构的散射体离散
右上：总场分布，显示入射波与散射波的叠加
左下：散射场，显示散射体产生的辐射模式
右下：边界上的入射场分布

10. 天线辐射问题仿真

10.1 偶极子天线模型

class AntennaRadiation:
    """天线辐射问题FEM-BEM仿真"""
    
    def __init__(self, frequency=1e9):
        self.frequency = frequency
        self.c = 3e8
        self.wavelength = self.c / frequency
        self.k = 2 * np.pi / self.wavelength

10.2 天线几何建模

def create_dipole_antenna(self, length, n_segments=20):
    """创建偶极子天线模型"""
    z = np.linspace(-length/2, length/2, n_segments)
    x = np.zeros_like(z)
    y = np.zeros_like(z)
    
    antenna_nodes = np.column_stack([x, y, z])
    
    # 创建周围介质区域 (2D截面)
    r_domain = length * 2
    n_rings = 10
    theta = np.linspace(0, 2*np.pi, 40)
    
    all_nodes_2d = []
    for ring in range(n_rings + 1):
        r = r_domain * ring / n_rings
        ring_x = r * np.cos(theta)
        ring_y = r * np.sin(theta)
        all_nodes_2d.append(np.column_stack([ring_x, ring_y]))
    
    nodes_2d = np.vstack(all_nodes_2d)
    
    return antenna_nodes, nodes_2d

10.3 辐射方向图计算

def compute_radiation_pattern(self, antenna_nodes, current_distribution):
    """计算辐射方向图"""
    theta_obs = np.linspace(0, 2*np.pi, 360)
    r_far = 100 * self.wavelength
    
    E_theta = np.zeros(len(theta_obs), dtype=complex)
    
    for i, theta in enumerate(theta_obs):
        x_obs = r_far * np.cos(theta)
        y_obs = r_far * np.sin(theta)
        
        for j in range(len(antenna_nodes) - 1):
            x_mid = (antenna_nodes[j, 0] + antenna_nodes[j+1, 0]) / 2
            y_mid = (antenna_nodes[j, 1] + antenna_nodes[j+1, 1]) / 2
            z_mid = (antenna_nodes[j, 2] + antenna_nodes[j+1, 2]) / 2
            
            r = np.sqrt((x_obs - x_mid)**2 + (y_obs - y_mid)**2 + z_mid**2)
            
            if r > 0:
                E_theta[i] += current_distribution[j] * \
                              np.exp(-1j * self.k * r) / r
    
    return theta_obs, np.abs(E_theta)

10.4 天线辐射结果

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了半波偶极子天线的辐射特性：

左图：天线几何结构，沿z轴放置
中图：正弦电流分布，符合理论预期
右图：辐射方向图，显示典型的"8"字形图案

11. 快速多极子加速技术

11.1 BEM的计算瓶颈

传统BEM的主要瓶颈：

稠密矩阵：存储复杂度 $O(N^2)$
矩阵向量乘法：计算复杂度 $O(N^2)$
直接求解： $O(N^3)$

11.2 快速多极子方法(FMM)

FMM通过分层分组和近似计算，将复杂度降至 $\log N)$ 或 $O (N)$ 。

核心思想：

远场近似：远距离组的相互作用用多极展开近似
分层结构：使用八叉树（3D）或四叉树（2D）组织空间
转换操作：M2M、M2L、L2L操作传递场信息

11.3 FMM在FEM-BEM中的应用

class FastMultipoleBEM:
    """快速多极子加速的BEM"""
    
    def __init__(self, boundary_nodes, k, n_levels=4):
        self.boundary_nodes = boundary_nodes
        self.k = k
        self.n_levels = n_levels
        self.tree = self._build_tree()
    
    def _build_tree(self):
        """构建空间分层树"""
        # 实现四叉树/八叉树构建
        pass
    
    def compute_matrix_vector_product(self, x):
        """快速矩阵向量乘法"""
        # 上行遍历：计算多极矩
        # 下行遍历：计算局部展开
        # 直接计算：近场相互作用
        pass

12. 数值算例与验证

12.1 验证方法

解析解对比：
对于简单几何（球、圆柱），与Mie级数解析解对比。

收敛性分析：
$∥uh−uexact∥≤Chp\|u_h - u_{exact}\| \leq Ch^p$

其中 $h$ 是网格尺寸， $p$ 是收敛阶。

12.2 误差分析

def compute_error(u_numerical, u_exact, norm_type='L2'):
    """计算数值误差"""
    if norm_type == 'L2':
        error = np.sqrt(np.sum(np.abs(u_numerical - u_exact)**2))
    elif norm_type == 'H1':
        # 包含梯度误差
        error = compute_h1_error(u_numerical, u_exact)
    
    return error

13. 工程应用案例

13.1 天线设计与分析

FEM-BEM耦合在天线设计中的应用：

复杂介质天线：介质透镜、人工电磁材料天线
共形天线：安装在平台表面的天线
阵列天线：考虑互耦效应的阵列分析

13.2 电磁散射与隐身

雷达截面(RCS)计算：飞机、舰船等目标的RCS分析
隐身设计：吸波材料与外形优化
气象雷达：雨滴、冰雹等粒子的散射

13.3 电磁兼容(EMC)

辐射发射分析：设备的电磁辐射评估
敏感度分析：外部电磁干扰的影响
屏蔽效能计算：屏蔽结构的性能评估

14. 方法比较与选择

14.1 FEM vs BEM vs FEM-BEM

特性	FEM	BEM	FEM-BEM
适用问题	有界区域	开域问题	复杂介质+开域
网格要求	区域网格	边界网格	区域+边界
矩阵类型	稀疏	稠密	混合
远场计算	需ABC/PML	自动满足	自动满足
计算复杂度	$O(N^{1.5})$	$O(N^3)$	$O(N^2)$
内存需求	$O (N)$	$O(N^2)$	$O(N^{1.5})$

14.2 效率比较

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图展示了三种方法的性能对比：

左图：计算时间随问题规模的变化
右图：内存使用随问题规模的变化

14.3 综合对比

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

上图从多个维度比较了三种方法：

左上：不同问题类型的适用性
右上：计算复杂度对比
左下：内存使用对比
右下：应用领域适用性

15. 前沿发展

15.1 快速算法

快速多极子FMM：
将BEM复杂度从 $O(N^3)$ 降至 $\log N)$ 。

H矩阵方法：
利用矩阵块的低秩特性进行压缩存储。

快速直接求解器：
利用层次矩阵结构构建近似逆矩阵。

15.2 高阶方法

高阶FEM：
使用p型有限元提高精度，减少自由度。

高阶BEM：
采用NURBS等参数曲面精确描述几何。

等几何分析(IGA)：
CAD与CAE的无缝集成。

15.3 多物理场耦合

电磁-热耦合：
考虑电磁损耗引起的温升效应。

电磁-结构耦合：
分析电磁力对结构的影响。

电磁-流体耦合：
等离子体与电磁场的相互作用。

16. 总结

16.1 核心要点回顾

FEM-BEM耦合方法是解决开域电磁场问题的强大工具：

优势互补：结合FEM处理复杂介质的能力和BEM处理无限域的优势
精确建模：避免人工边界带来的近似误差
灵活应用：适用于天线、散射、EMC等多种工程问题

16.2 实现要点

网格生成：确保耦合界面处网格匹配
矩阵组装：正确处理耦合块矩阵
求解策略：利用矩阵结构选择高效求解器
后处理：提取远场辐射/散射特性

16.3 学习建议

理论基础：深入理解变分原理和边界积分方程
编程实践：从简单问题开始，逐步增加复杂度
文献阅读：关注IEEE TAP、JCP等顶级期刊
软件使用：学习商业软件（HFSS、CST）中的混合方法

17. 参考文献

Steinbach, O. (2008). Numerical Approximation Methods for Elliptic Boundary Value Problems: Finite and Boundary Elements. Springer.
McLean, W. (2000). Strongly Elliptic Systems and Boundary Integral Equations. Cambridge University Press.
Nédélec, J. C. (2001). Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems. Springer.
Sheng, X. Q. (2013). Essential Wavelets for Electromagnetic Applications. Artech House.
Jin, J. M. (2015). The Finite Element Method in Electromagnetics (3rd ed.). Wiley-IEEE Press.
Peterson, A. F., Ray, S. L., & Mittra, R. (1998). Computational Methods for Electromagnetics. IEEE Press.
Liu, Y. J. (2009). Fast Multipole Boundary Element Method: Theory and Applications in Engineering. Cambridge University Press.
Costabel, M. (1987). Symmetric methods for the coupling of finite elements and boundary elements. Boundary Elements IX, 411-420.

18. 附录

# -*- coding: utf-8 -*-
"""
主题078: FEM-BEM耦合方法
有限元-边界元耦合方法用于电磁场开域问题的数值仿真
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse import linalg as splinalg
from scipy.special import hankel1, jv
from scipy.integrate import quad
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class FEM2D:
    """二维有限元求解器 - 处理内部区域"""
    
    def __init__(self, nodes, elements, epsilon_r=1.0, mu_r=1.0):
        """
        参数:
            nodes: 节点坐标数组 (N, 2)
            elements: 三角形单元连接数组 (M, 3)
            epsilon_r: 相对介电常数
            mu_r: 相对磁导率
        """
        self.nodes = nodes
        self.elements = elements
        self.n_nodes = len(nodes)
        self.n_elements = len(elements)
        self.epsilon_r = epsilon_r
        self.mu_r = mu_r
        self.k = 2 * np.pi  # 波数 (假设频率为1，c=1)
        
    def compute_element_matrices(self, element_idx):
        """计算单元刚度矩阵和质量矩阵"""
        elem = self.elements[element_idx]
        x = self.nodes[elem, 0]
        y = self.nodes[elem, 1]
        
        # 计算单元面积
        area = 0.5 * abs((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]))
        
        # 形状函数梯度
        b = np.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]])
        c = np.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]])
        
        # 刚度矩阵 (梯度项)
        K_elem = np.outer(b, b) + np.outer(c, c)
        K_elem = K_elem / (4 * area)
        
        # 质量矩阵 (势能项)
        M_elem = area / 12 * np.array([[2, 1, 1],
                                        [1, 2, 1],
                                        [1, 1, 2]])
        
        return K_elem, M_elem, area
    
    def assemble_system(self, k=None):
        """组装全局系统矩阵"""
        if k is None:
            k = self.k
        
        # 全局矩阵
        rows, cols, data_K, data_M = [], [], [], []
        
        for e in range(self.n_elements):
            K_elem, M_elem, _ = self.compute_element_matrices(e)
            elem = self.elements[e]
            
            for i in range(3):
                for j in range(3):
                    rows.append(elem[i])
                    cols.append(elem[j])
                    data_K.append(K_elem[i, j])
                    data_M.append(M_elem[i, j])
        
        K_global = sparse.coo_matrix((data_K, (rows, cols)), 
                                      shape=(self.n_nodes, self.n_nodes)).tocsr()
        M_global = sparse.coo_matrix((data_M, (rows, cols)), 
                                      shape=(self.n_nodes, self.n_nodes)).tocsr()
        
        # Helmholtz方程: -∇²u - k²u = f
        A = K_global - k**2 * M_global
        
        return A, K_global, M_global
    
    def apply_dirichlet_boundary(self, A, b, bc_nodes, bc_values):
        """应用Dirichlet边界条件"""
        A = A.tolil()
        
        for node, value in zip(bc_nodes, bc_values):
            A[node, :] = 0
            A[node, node] = 1
            b[node] = value
        
        return A.tocsr(), b


class BEM2D:
    """二维边界元求解器 - 处理无限远场"""
    
    def __init__(self, boundary_nodes, k=2*np.pi):
        """
        参数:
            boundary_nodes: 边界节点坐标 (N, 2)
            k: 波数
        """
        self.boundary_nodes = boundary_nodes
        self.n_bound = len(boundary_nodes)
        self.k = k
        
    def green_function(self, x, y, x0, y0):
        """二维Helmholtz格林函数"""
        r = np.sqrt((x - x0)**2 + (y - y0)**2)
        if r < 1e-10:
            return 0.0
        return 0.25j * hankel1(0, self.k * r)
    
    def green_normal_derivative(self, x, y, x0, y0, nx, ny):
        """格林函数法向导数"""
        r = np.sqrt((x - x0)**2 + (y - y0)**2)
        if r < 1e-10:
            return 0.0
        dr_dn = ((x - x0) * nx + (y - y0) * ny) / r
        return -0.25j * self.k * hankel1(1, self.k * r) * dr_dn
    
    def compute_boundary_matrices(self):
        """计算边界元系统矩阵"""
        n = self.n_bound
        H = np.zeros((n, n), dtype=complex)
        G = np.zeros((n, n), dtype=complex)
        
        # 计算边界法向量
        normals = self._compute_normals()
        
        for i in range(n):
            xi, yi = self.boundary_nodes[i]
            
            for j in range(n):
                xj, yj = self.boundary_nodes[j]
                
                if i == j:
                    # 对角项处理 (奇异性)
                    H[i, j] = 0.5
                    # 计算G的对角项积分
                    G[i, j] = self._compute_singular_integral(i)
                else:
                    # 非对角项
                    r = np.sqrt((xi - xj)**2 + (yi - yj)**2)
                    H[i, j] = -0.25j * self.k * hankel1(1, self.k * r)
                    G[i, j] = 0.25j * hankel1(0, self.k * r)
        
        return H, G, normals
    
    def _compute_normals(self):
        """计算边界外法向量"""
        n = self.n_bound
        normals = np.zeros((n, 2))
        
        for i in range(n):
            i_prev = (i - 1) % n
            i_next = (i + 1) % n
            
            dx = self.boundary_nodes[i_next, 0] - self.boundary_nodes[i_prev, 0]
            dy = self.boundary_nodes[i_next, 1] - self.boundary_nodes[i_prev, 1]
            
            length = np.sqrt(dx**2 + dy**2)
            if length > 0:
                # 外法向量 (顺时针旋转90度)
                normals[i, 0] = dy / length
                normals[i, 1] = -dx / length
        
        return normals
    
    def _compute_singular_integral(self, idx):
        """计算奇异积分 (对角项)"""
        # 简化的奇异积分处理
        # 使用对数奇异性的解析结果
        return 1.0 / (2 * np.pi) * np.log(1.0 / self.k)


class FEMBEMCoupling:
    """FEM-BEM耦合求解器"""
    
    def __init__(self, fem_solver, bem_solver, coupling_nodes_fem, 
                 coupling_nodes_bem):
        """
        参数:
            fem_solver: FEM求解器实例
            bem_solver: BEM求解器实例
            coupling_nodes_fem: FEM侧耦合节点索引
            coupling_nodes_bem: BEM侧耦合节点索引
        """
        self.fem = fem_solver
        self.bem = bem_solver
        self.coupling_fem = coupling_nodes_fem
        self.coupling_bem = coupling_nodes_bem
        self.n_coupling = len(coupling_nodes_fem)
        
    def build_coupled_system(self, source_nodes=None, source_values=None):
        """构建耦合系统矩阵"""
        # FEM系统
        A_fem, K_fem, M_fem = self.fem.assemble_system()
        n_fem = self.fem.n_nodes
        
        # BEM系统
        H_bem, G_bem, normals = self.bem.compute_boundary_matrices()
        n_bem = self.bem.n_bound
        
        # 构建耦合矩阵
        # 系统形式: [A_fem   C1  ][u_fem]   [f_fem]
        #           [C2    A_bem][u_bem] = [f_bem]
        
        n_total = n_fem + n_bem
        A_coupled = np.zeros((n_total, n_total), dtype=complex)
        b_coupled = np.zeros(n_total, dtype=complex)
        
        # FEM块
        A_coupled[:n_fem, :n_fem] = A_fem.toarray()
        
        # BEM块
        A_coupled[n_fem:, n_fem:] = H_bem
        
        # 耦合块 (简化处理)
        # 在耦合边界上，FEM和BEM的解应该连续
        for i, (fem_idx, bem_idx) in enumerate(zip(self.coupling_fem, 
                                                    self.coupling_bem)):
            # FEM到BEM的耦合
            A_coupled[n_fem + bem_idx, fem_idx] = -1.0
            # BEM到FEM的耦合
            A_coupled[fem_idx, n_fem + bem_idx] = -1.0
        
        # 源项
        if source_nodes is not None and source_values is not None:
            b_coupled[source_nodes] = source_values
        
        return A_coupled, b_coupled
    
    def solve(self, source_nodes=None, source_values=None):
        """求解耦合系统"""
        A, b = self.build_coupled_system(source_nodes, source_values)
        
        # 求解线性系统
        u = np.linalg.solve(A, b)
        
        # 分离FEM和BEM解
        n_fem = self.fem.n_nodes
        u_fem = u[:n_fem]
        u_bem = u[n_fem:]
        
        return u_fem, u_bem


class ScatteringProblem:
    """电磁散射问题仿真"""
    
    def __init__(self, frequency=1e9, epsilon_r=4.0, mu_r=1.0):
        """
        参数:
            frequency: 频率 (Hz)
            epsilon_r: 散射体相对介电常数
            mu_r: 相对磁导率
        """
        self.frequency = frequency
        self.c = 3e8
        self.wavelength = self.c / frequency
        self.k0 = 2 * np.pi / self.wavelength
        self.epsilon_r = epsilon_r
        self.mu_r = mu_r
        
    def create_circular_scatterer(self, radius, n_elements=40):
        """创建圆形散射体网格"""
        # 边界节点
        theta = np.linspace(0, 2*np.pi, n_elements, endpoint=False)
        boundary_nodes = np.column_stack([
            radius * np.cos(theta),
            radius * np.sin(theta)
        ])
        
        # 内部网格 (简化：使用同心圆)
        n_rings = 5
        all_nodes = [boundary_nodes]
        all_elements = []
        
        for ring in range(1, n_rings + 1):
            r = radius * (n_rings - ring) / n_rings
            ring_nodes = np.column_stack([
                r * np.cos(theta),
                r * np.sin(theta)
            ])
            all_nodes.append(ring_nodes)
        
        # 合并节点
        nodes = np.vstack(all_nodes[::-1])
        
        # 创建三角形单元
        n_per_ring = n_elements
        for ring in range(n_rings):
            for i in range(n_per_ring):
                i1 = ring * n_per_ring + i
                i2 = ring * n_per_ring + (i + 1) % n_per_ring
                i3 = (ring + 1) * n_per_ring + i
                i4 = (ring + 1) * n_per_ring + (i + 1) % n_per_ring
                
                all_elements.append([i1, i2, i3])
                if ring < n_rings - 1:
                    all_elements.append([i2, i4, i3])
        
        elements = np.array(all_elements)
        
        # 耦合节点索引 (最外层)
        coupling_fem = np.arange(n_rings * n_per_ring, (n_rings + 1) * n_per_ring)
        coupling_bem = np.arange(n_elements)
        
        return nodes, elements, boundary_nodes, coupling_fem, coupling_bem
    
    def incident_field(self, x, y, direction=0):
        """入射平面波"""
        kx = self.k0 * np.cos(direction)
        ky = self.k0 * np.sin(direction)
        return np.exp(1j * (kx * x + ky * y))
    
    def solve_scattering(self, radius=0.5, n_elements=40):
        """求解散射问题"""
        # 创建网格
        nodes, elements, boundary_nodes, coupling_fem, coupling_bem = \
            self.create_circular_scatterer(radius, n_elements)
        
        # 创建求解器
        fem = FEM2D(nodes, elements, self.epsilon_r, self.mu_r)
        fem.k = self.k0 * np.sqrt(self.epsilon_r * self.mu_r)
        
        bem = BEM2D(boundary_nodes, self.k0)
        
        # 创建耦合求解器
        coupling = FEMBEMCoupling(fem, bem, coupling_fem, coupling_bem)
        
        # 计算入射场在边界上的值
        x_b = boundary_nodes[:, 0]
        y_b = boundary_nodes[:, 1]
        u_inc = self.incident_field(x_b, y_b)
        
        # 求解 (简化：直接求解)
        # 这里使用简化的耦合求解
        A_fem, _, _ = fem.assemble_system()
        
        # 右端项
        b_fem = np.zeros(fem.n_nodes, dtype=complex)
        
        # 在边界上施加入射场条件
        A_fem_lil = A_fem.tolil()
        for idx in coupling_fem:
            A_fem_lil[idx, :] = 0
            A_fem_lil[idx, idx] = 1
        A_fem = A_fem_lil.tocsr()
        
        for i, idx in enumerate(coupling_fem):
            b_fem[idx] = u_inc[i]
        
        # 求解
        u_total = splinalg.spsolve(A_fem, b_fem)
        
        # 计算散射场
        u_scattered = u_total.copy()
        for i, idx in enumerate(coupling_fem):
            u_scattered[idx] -= u_inc[i]
        
        return nodes, elements, boundary_nodes, u_total, u_scattered, u_inc


class AntennaRadiation:
    """天线辐射问题FEM-BEM仿真"""
    
    def __init__(self, frequency=1e9):
        self.frequency = frequency
        self.c = 3e8
        self.wavelength = self.c / frequency
        self.k = 2 * np.pi / self.wavelength
        
    def create_dipole_antenna(self, length, n_segments=20):
        """创建偶极子天线模型"""
        # 天线沿z轴放置
        z = np.linspace(-length/2, length/2, n_segments)
        x = np.zeros_like(z)
        y = np.zeros_like(z)
        
        antenna_nodes = np.column_stack([x, y, z])
        
        # 创建周围介质区域 (简化为2D截面)
        r_domain = length * 2
        n_rings = 10
        theta = np.linspace(0, 2*np.pi, 40)
        
        all_nodes_2d = []
        for ring in range(n_rings + 1):
            r = r_domain * ring / n_rings
            ring_x = r * np.cos(theta)
            ring_y = r * np.sin(theta)
            all_nodes_2d.append(np.column_stack([ring_x, ring_y]))
        
        nodes_2d = np.vstack(all_nodes_2d)
        
        # 创建三角形单元
        elements = []
        n_theta = len(theta)
        for ring in range(n_rings):
            for i in range(n_theta):
                i1 = ring * n_theta + i
                i2 = ring * n_theta + (i + 1) % n_theta
                i3 = (ring + 1) * n_theta + i
                i4 = (ring + 1) * n_theta + (i + 1) % n_theta
                
                elements.append([i1, i2, i3])
                elements.append([i2, i4, i3])
        
        return antenna_nodes, nodes_2d, np.array(elements)
    
    def compute_radiation_pattern(self, antenna_nodes, current_distribution):
        """计算辐射方向图"""
        # 远场观测角度
        theta_obs = np.linspace(0, 2*np.pi, 360)
        r_far = 100 * self.wavelength  # 远场距离
        
        E_theta = np.zeros(len(theta_obs), dtype=complex)
        
        for i, theta in enumerate(theta_obs):
            # 观测点
            x_obs = r_far * np.cos(theta)
            y_obs = r_far * np.sin(theta)
            
            # 计算辐射场 (简化)
            for j in range(len(antenna_nodes) - 1):
                # 天线单元中点
                x_mid = (antenna_nodes[j, 0] + antenna_nodes[j+1, 0]) / 2
                y_mid = (antenna_nodes[j, 1] + antenna_nodes[j+1, 1]) / 2
                z_mid = (antenna_nodes[j, 2] + antenna_nodes[j+1, 2]) / 2
                
                # 距离
                r = np.sqrt((x_obs - x_mid)**2 + (y_obs - y_mid)**2 + z_mid**2)
                
                # 辐射积分
                if r > 0:
                    E_theta[i] += current_distribution[j] * \
                                  np.exp(-1j * self.k * r) / r
        
        return theta_obs, np.abs(E_theta)


class CouplingEfficiency:
    """FEM-BEM耦合效率分析"""
    
    def __init__(self):
        self.results = {}
        
    def compare_methods(self, problem_size_range):
        """比较FEM、BEM和FEM-BEM的效率"""
        fem_times = []
        bem_times = []
        fem_bem_times = []
        
        for n in problem_size_range:
            # 简化的性能估计
            # FEM: O(n^2) for sparse, O(n^1.5) with good preconditioner
            fem_time = n**1.5 * 1e-6
            
            # BEM: O(n^3) for dense, O(n^2) with fast methods
            bem_time = n**2.5 * 1e-6
            
            # FEM-BEM: 介于两者之间
            fem_bem_time = n**2 * 1e-6
            
            fem_times.append(fem_time)
            bem_times.append(bem_time)
            fem_bem_times.append(fem_bem_time)
        
        self.results = {
            'problem_size': problem_size_range,
            'fem_time': fem_times,
            'bem_time': bem_times,
            'fem_bem_time': fem_bem_times
        }
        
        return self.results
    
    def memory_comparison(self, problem_size_range):
        """内存使用比较"""
        fem_memory = []
        bem_memory = []
        fem_bem_memory = []
        
        for n in problem_size_range:
            # FEM: 稀疏矩阵，O(n)
            fem_mem = n * 50  # 50非零元每行
            
            # BEM: 稠密矩阵，O(n^2)
            bem_mem = n**2
            
            # FEM-BEM: 混合
            fem_bem_mem = n * 50 + n**2 * 0.1
            
            fem_memory.append(fem_mem)
            bem_memory.append(bem_mem)
            fem_bem_memory.append(fem_bem_mem)
        
        return {
            'problem_size': problem_size_range,
            'fem_memory': fem_memory,
            'bem_memory': bem_memory,
            'fem_bem_memory': fem_bem_memory
        }


def demo_fem_solver():
    """演示FEM求解器"""
    print("=" * 60)
    print("演示1: 2D有限元求解器")
    print("=" * 60)
    
    # 创建简单网格 (正方形区域)
    n = 20
    x = np.linspace(-1, 1, n)
    y = np.linspace(-1, 1, n)
    X, Y = np.meshgrid(x, y)
    
    nodes = np.column_stack([X.flatten(), Y.flatten()])
    
    # 创建三角形单元
    elements = []
    for i in range(n-1):
        for j in range(n-1):
            n1 = i * n + j
            n2 = i * n + (j + 1)
            n3 = (i + 1) * n + j
            n4 = (i + 1) * n + (j + 1)
            
            elements.append([n1, n2, n3])
            elements.append([n2, n4, n3])
    
    elements = np.array(elements)
    
    # 创建FEM求解器
    fem = FEM2D(nodes, elements, epsilon_r=1.0)
    
    # 组装系统
    A, K, M = fem.assemble_system()
    
    print(f"网格节点数: {fem.n_nodes}")
    print(f"单元数: {fem.n_elements}")
    print(f"系统矩阵大小: {A.shape}")
    print(f"矩阵非零元素: {A.nnz}")
    
    # 求解简单问题
    b = np.zeros(fem.n_nodes, dtype=complex)
    
    # 在中心设置源
    center_idx = (n//2) * n + n//2
    b[center_idx] = 1.0
    
    # 边界条件 (Dirichlet)
    boundary_nodes = []
    for i in range(n):
        boundary_nodes.extend([i, i + n*(n-1), i*n, i*n + (n-1)])
    boundary_nodes = list(set(boundary_nodes))
    
    bc_values = np.zeros(len(boundary_nodes))
    A, b = fem.apply_dirichlet_boundary(A, b, boundary_nodes, bc_values)
    
    # 求解
    u = splinalg.spsolve(A, b)
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 网格
    ax = axes[0]
    ax.triplot(nodes[:, 0], nodes[:, 1], elements, 'b-', lw=0.5)
    ax.set_aspect('equal')
    ax.set_title('FEM Mesh')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    
    # 解
    ax = axes[1]
    u_plot = np.real(u).reshape(n, n)
    im = ax.contourf(X, Y, u_plot, levels=20, cmap='RdBu_r')
    ax.set_aspect('equal')
    ax.set_title('FEM Solution (Real part)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(im, ax=ax)
    
    plt.tight_layout()
    plt.savefig('fem_solver_demo.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 fem_solver_demo.png")
    plt.close()
    
    return fem, u


def demo_bem_solver():
    """演示BEM求解器"""
    print("\n" + "=" * 60)
    print("演示2: 2D边界元求解器")
    print("=" * 60)
    
    # 创建圆形边界
    n = 40
    theta = np.linspace(0, 2*np.pi, n, endpoint=False)
    radius = 1.0
    boundary_nodes = np.column_stack([
        radius * np.cos(theta),
        radius * np.sin(theta)
    ])
    
    # 创建BEM求解器
    k = 2 * np.pi  # 波数
    bem = BEM2D(boundary_nodes, k)
    
    # 计算边界矩阵
    H, G, normals = bem.compute_boundary_matrices()
    
    print(f"边界节点数: {bem.n_bound}")
    print(f"波数 k = {k:.4f}")
    print(f"波长 λ = {2*np.pi/k:.4f}")
    print(f"H矩阵条件数: {np.linalg.cond(H):.2e}")
    print(f"G矩阵条件数: {np.linalg.cond(G):.2e}")
    
    # 可视化边界和法向量
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax = axes[0]
    ax.plot(boundary_nodes[:, 0], boundary_nodes[:, 1], 'bo-', markersize=4)
    # 绘制法向量
    scale = 0.15
    for i in range(0, n, 4):
        ax.arrow(boundary_nodes[i, 0], boundary_nodes[i, 1],
                scale * normals[i, 0], scale * normals[i, 1],
                head_width=0.05, color='r')
    ax.set_aspect('equal')
    ax.set_title('BEM Boundary with Normals')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid(True)
    
    # 矩阵可视化
    ax = axes[1]
    im = ax.imshow(np.abs(H), cmap='viridis', aspect='auto')
    ax.set_title('|H| Matrix')
    ax.set_xlabel('Node index')
    ax.set_ylabel('Node index')
    plt.colorbar(im, ax=ax)
    
    plt.tight_layout()
    plt.savefig('bem_solver_demo.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 bem_solver_demo.png")
    plt.close()
    
    return bem, H, G


def demo_scattering():
    """演示散射问题"""
    print("\n" + "=" * 60)
    print("演示3: 电磁散射问题 (FEM-BEM)")
    print("=" * 60)
    
    # 创建散射问题
    scatterer = ScatteringProblem(frequency=1e9, epsilon_r=4.0)
    
    print(f"频率: {scatterer.frequency/1e9:.1f} GHz")
    print(f"波长: {scatterer.wavelength:.4f} m")
    print(f"波数: {scatterer.k0:.4f} rad/m")
    print(f"散射体介电常数: {scatterer.epsilon_r}")
    
    # 求解
    radius = 0.5  # 散射体半径
    nodes, elements, boundary_nodes, u_total, u_scattered, u_inc = \
        scatterer.solve_scattering(radius, n_elements=40)
    
    print(f"\n网格统计:")
    print(f"  节点数: {len(nodes)}")
    print(f"  单元数: {len(elements)}")
    print(f"  边界节点数: {len(boundary_nodes)}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 网格
    ax = axes[0, 0]
    ax.triplot(nodes[:, 0], nodes[:, 1], elements, 'b-', lw=0.5)
    ax.plot(boundary_nodes[:, 0], boundary_nodes[:, 1], 'r-', lw=2)
    ax.set_aspect('equal')
    ax.set_title('Scattering Mesh')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    
    # 总场
    ax = axes[0, 1]
    u_total_real = np.real(u_total)
    scatter = ax.scatter(nodes[:, 0], nodes[:, 1], c=u_total_real, 
                         cmap='RdBu_r', s=10)
    ax.set_aspect('equal')
    ax.set_title('Total Field (Real part)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(scatter, ax=ax)
    
    # 散射场
    ax = axes[1, 0]
    u_scat_real = np.real(u_scattered)
    scatter = ax.scatter(nodes[:, 0], nodes[:, 1], c=u_scat_real, 
                         cmap='RdBu_r', s=10)
    ax.set_aspect('equal')
    ax.set_title('Scattered Field (Real part)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(scatter, ax=ax)
    
    # 入射场
    ax = axes[1, 1]
    u_inc_real = np.real(u_inc)
    scatter = ax.scatter(boundary_nodes[:, 0], boundary_nodes[:, 1], 
                         c=u_inc_real, cmap='RdBu_r', s=30)
    ax.set_aspect('equal')
    ax.set_title('Incident Field on Boundary')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(scatter, ax=ax)
    
    plt.tight_layout()
    plt.savefig('scattering_fem_bem.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 scattering_fem_bem.png")
    plt.close()
    
    return scatterer, nodes, elements, u_total, u_scattered


def demo_radiation_pattern():
    """演示天线辐射方向图"""
    print("\n" + "=" * 60)
    print("演示4: 天线辐射方向图")
    print("=" * 60)
    
    # 创建天线模型
    antenna = AntennaRadiation(frequency=1e9)
    
    # 半波偶极子
    length = antenna.wavelength / 2
    antenna_nodes, nodes_2d, elements = antenna.create_dipole_antenna(length)
    
    print(f"天线长度: {length:.4f} m (半波长)")
    print(f"波长: {antenna.wavelength:.4f} m")
    print(f"频率: {antenna.frequency/1e9:.1f} GHz")
    
    # 假设正弦电流分布
    z_coords = antenna_nodes[:, 2]
    current_dist = np.sin(2 * np.pi / antenna.wavelength * 
                          (length/2 - np.abs(z_coords)))
    
    # 计算辐射方向图
    theta, E_theta = antenna.compute_radiation_pattern(antenna_nodes, 
                                                        current_dist)
    
    # 归一化
    E_theta_norm = np.abs(E_theta) / np.max(np.abs(E_theta))
    E_theta_db = 20 * np.log10(E_theta_norm + 1e-10)
    
    print(f"最大辐射方向: {theta[np.argmax(E_theta_norm)] * 180/np.pi:.1f}°")
    
    # 可视化
    fig, axes = plt.subplots(1, 3, figsize=(16, 5))
    
    # 天线几何
    ax = axes[0]
    ax.plot(antenna_nodes[:, 2], np.zeros(len(antenna_nodes)), 'b-o', 
            markersize=3)
    ax.set_title('Dipole Antenna Geometry')
    ax.set_xlabel('z (m)')
    ax.set_ylabel('x, y')
    ax.grid(True)
    ax.set_aspect('equal')
    
    # 电流分布
    ax = axes[1]
    ax.plot(z_coords, current_dist, 'b-', lw=2)
    ax.set_title('Current Distribution')
    ax.set_xlabel('z (m)')
    ax.set_ylabel('Current (A)')
    ax.grid(True)
    
    # 辐射方向图 (极坐标)
    ax = axes[2]
    ax = plt.subplot(1, 3, 3, projection='polar')
    ax.plot(theta, E_theta_db, 'b-', lw=2)
    ax.set_title('Radiation Pattern (dB)', pad=20)
    ax.set_ylim(-40, 5)
    ax.grid(True)
    
    plt.tight_layout()
    plt.savefig('antenna_radiation.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 antenna_radiation.png")
    plt.close()
    
    return antenna, theta, E_theta_db


def demo_efficiency_comparison():
    """演示效率比较"""
    print("\n" + "=" * 60)
    print("演示5: FEM-BEM效率比较")
    print("=" * 60)
    
    efficiency = CouplingEfficiency()
    
    # 问题规模范围
    problem_sizes = np.array([100, 200, 500, 1000, 2000, 5000])
    
    # 计算时间比较
    time_results = efficiency.compare_methods(problem_sizes)
    
    # 内存使用比较
    memory_results = efficiency.memory_comparison(problem_sizes)
    
    print("计算时间比较 (相对值):")
    print("-" * 50)
    print(f"{'Size':>10} {'FEM':>12} {'BEM':>12} {'FEM-BEM':>12}")
    print("-" * 50)
    for i, n in enumerate(problem_sizes):
        print(f"{n:>10} {time_results['fem_time'][i]:>12.2e} "
              f"{time_results['bem_time'][i]:>12.2e} "
              f"{time_results['fem_bem_time'][i]:>12.2e}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 计算时间
    ax = axes[0]
    ax.loglog(problem_sizes, time_results['fem_time'], 'o-', label='FEM', lw=2)
    ax.loglog(problem_sizes, time_results['bem_time'], 's-', label='BEM', lw=2)
    ax.loglog(problem_sizes, time_results['fem_bem_time'], '^-', 
              label='FEM-BEM', lw=2)
    ax.set_xlabel('Problem Size (DOF)')
    ax.set_ylabel('Computation Time (s)')
    ax.set_title('Computational Time Comparison')
    ax.legend()
    ax.grid(True, which='both', ls='--', alpha=0.5)
    
    # 内存使用
    ax = axes[1]
    ax.loglog(problem_sizes, memory_results['fem_memory'], 'o-', label='FEM', 
              lw=2)
    ax.loglog(problem_sizes, memory_results['bem_memory'], 's-', label='BEM', 
              lw=2)
    ax.loglog(problem_sizes, memory_results['fem_bem_memory'], '^-', 
              label='FEM-BEM', lw=2)
    ax.set_xlabel('Problem Size (DOF)')
    ax.set_ylabel('Memory Usage (elements)')
    ax.set_title('Memory Usage Comparison')
    ax.legend()
    ax.grid(True, which='both', ls='--', alpha=0.5)
    
    plt.tight_layout()
    plt.savefig('efficiency_comparison.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存到 efficiency_comparison.png")
    plt.close()
    
    return efficiency, time_results, memory_results


def demo_coupling_interface():
    """演示FEM-BEM耦合界面"""
    print("\n" + "=" * 60)
    print("演示6: FEM-BEM耦合界面处理")
    print("=" * 60)
    
    # 创建耦合界面示例
    n_interface = 32
    theta = np.linspace(0, 2*np.pi, n_interface, endpoint=False)
    radius = 1.0
    
    # FEM侧 (内部，更密的网格)
    n_fem_ring = 3
    fem_nodes = []
    for ring in range(n_fem_ring + 1):
        r = radius * (n_fem_ring - ring) / n_fem_ring
        ring_nodes = np.column_stack([
            r * np.cos(theta),
            r * np.sin(theta)
        ])
        fem_nodes.append(ring_nodes)
    fem_nodes = np.vstack(fem_nodes)
    
    # BEM侧 (边界)
    bem_nodes = np.column_stack([
        radius * np.cos(theta),
        radius * np.sin(theta)
    ])
    
    # 耦合映射
    coupling_map = list(range(n_interface))
    
    print(f"FEM节点数: {len(fem_nodes)}")
    print(f"BEM节点数: {len(bem_nodes)}")
    print(f"耦合界面节点数: {n_interface}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 耦合界面
    ax = axes[0]
    ax.scatter(fem_nodes[:, 0], fem_nodes[:, 1], c='blue', s=30, 
               label='FEM nodes', alpha=0.6)
    ax.scatter(bem_nodes[:, 0], bem_nodes[:, 1], c='red', s=50, 
               label='BEM nodes', marker='s')
    
    # 绘制耦合连线
    for i in range(0, n_interface, 4):
        ax.plot([fem_nodes[n_fem_ring*n_interface + i, 0], bem_nodes[i, 0]],
                [fem_nodes[n_fem_ring*n_interface + i, 1], bem_nodes[i, 1]],
                'g--', lw=1, alpha=0.5)
    
    circle = plt.Circle((0, 0), radius, fill=False, color='k', ls='--', lw=2)
    ax.add_patch(circle)
    ax.set_aspect('equal')
    ax.set_title('FEM-BEM Coupling Interface')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.legend()
    ax.grid(True)
    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    
    # 耦合矩阵示意图
    ax = axes[1]
    n_total = len(fem_nodes) + len(bem_nodes)
    coupling_matrix = np.zeros((n_total, n_total))
    
    # FEM块
    n_fem = len(fem_nodes)
    coupling_matrix[:n_fem, :n_fem] = 0.3
    
    # BEM块
    coupling_matrix[n_fem:, n_fem:] = 0.5
    
    # 耦合块
    for i in range(n_interface):
        fem_idx = n_fem - n_interface + i
        bem_idx = n_fem + i
        coupling_matrix[fem_idx, bem_idx] = 0.8
        coupling_matrix[bem_idx, fem_idx] = 0.8
    
    im = ax.imshow(coupling_matrix, cmap='Blues', aspect='auto')
    ax.axhline(n_fem - 0.5, color='r', lw=2)
    ax.axvline(n_fem - 0.5, color='r', lw=2)
    ax.set_title('Coupled System Matrix Structure')
    ax.set_xlabel('Matrix Column')
    ax.set_ylabel('Matrix Row')
    
    # 添加文本标注
    ax.text(n_fem//2, n_fem//2, 'FEM', ha='center', va='center', 
            fontsize=14, color='white', fontweight='bold')
    ax.text(n_fem + len(bem_nodes)//2, n_fem + len(bem_nodes)//2, 'BEM', 
            ha='center', va='center', fontsize=14, color='white', fontweight='bold')
    ax.text(n_fem//2, n_fem + len(bem_nodes)//2, 'Coupling', ha='center', 
            va='center', fontsize=12, color='red', fontweight='bold')
    ax.text(n_fem + len(bem_nodes)//2, n_fem//2, 'Coupling', ha='center', 
            va='center', fontsize=12, color='red', fontweight='bold')
    
    plt.colorbar(im, ax=ax)
    
    plt.tight_layout()
    plt.savefig('coupling_interface.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 coupling_interface.png")
    plt.close()
    
    return fem_nodes, bem_nodes, coupling_map


def demo_green_function():
    """演示格林函数"""
    print("\n" + "=" * 60)
    print("演示7: Helmholtz格林函数")
    print("=" * 60)
    
    # 创建网格
    x = np.linspace(-2, 2, 100)
    y = np.linspace(-2, 2, 100)
    X, Y = np.meshgrid(x, y)
    
    # 源点位置
    x0, y0 = 0.5, -0.3
    
    # 计算格林函数
    k = 2 * np.pi
    r = np.sqrt((X - x0)**2 + (Y - y0)**2)
    r = np.maximum(r, 0.01)  # 避免奇异性
    
    G = 0.25j * hankel1(0, k * r)
    
    print(f"源点位置: ({x0}, {y0})")
    print(f"波数 k = {k:.4f}")
    print(f"格林函数最大值: {np.max(np.abs(G)):.4f}")
    print(f"格林函数最小值: {np.min(np.abs(G)):.4f}")
    
    # 可视化
    fig, axes = plt.subplots(1, 3, figsize=(16, 5))
    
    # 实部
    ax = axes[0]
    im = ax.contourf(X, Y, np.real(G), levels=20, cmap='RdBu_r')
    ax.plot(x0, y0, 'ko', markersize=10)
    ax.set_aspect('equal')
    ax.set_title('Green Function (Real part)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(im, ax=ax)
    
    # 虚部
    ax = axes[1]
    im = ax.contourf(X, Y, np.imag(G), levels=20, cmap='RdBu_r')
    ax.plot(x0, y0, 'ko', markersize=10)
    ax.set_aspect('equal')
    ax.set_title('Green Function (Imaginary part)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(im, ax=ax)
    
    # 幅度
    ax = axes[2]
    im = ax.contourf(X, Y, np.abs(G), levels=20, cmap='viridis')
    ax.plot(x0, y0, 'ro', markersize=10)
    ax.set_aspect('equal')
    ax.set_title('|G| (Magnitude)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(im, ax=ax)
    
    plt.tight_layout()
    plt.savefig('green_function.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 green_function.png")
    plt.close()
    
    return X, Y, G


def demo_fem_bem_summary():
    """FEM-BEM方法总结"""
    print("\n" + "=" * 60)
    print("演示8: FEM-BEM方法综合对比")
    print("=" * 60)
    
    # 方法特性对比
    methods = ['FEM', 'BEM', 'FEM-BEM Coupling']
    characteristics = {
        '适用问题': [
            '有界区域、复杂介质',
            '开域问题、均匀介质',
            '复杂介质+开域边界'
        ],
        '网格要求': [
            '区域网格',
            '边界网格',
            '区域+边界网格'
        ],
        '矩阵类型': [
            '稀疏矩阵',
            '稠密矩阵',
            '稀疏+稠密块'
        ],
        '远场计算': [
            '需要吸收边界',
            '自动满足',
            '自动满足'
        ],
        '计算复杂度': [
            'O(N^1.5)',
            'O(N^3)',
            'O(N^2)'
        ]
    }
    
    print("\n方法特性对比:")
    print("-" * 80)
    for key, values in characteristics.items():
        print(f"\n{key}:")
        for method, value in zip(methods, values):
            print(f"  {method:20s}: {value}")
    
    # 可视化对比
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 适用问题类型
    ax = axes[0, 0]
    categories = ['有界区域', '开域问题', '复杂介质', '均匀介质']
    fem_score = [1, 0.3, 1, 0.5]
    bem_score = [0.3, 1, 0.4, 1]
    fem_bem_score = [0.8, 1, 1, 0.8]
    
    x = np.arange(len(categories))
    width = 0.25
    
    ax.bar(x - width, fem_score, width, label='FEM', color='blue', alpha=0.7)
    ax.bar(x, bem_score, width, label='BEM', color='red', alpha=0.7)
    ax.bar(x + width, fem_bem_score, width, label='FEM-BEM', color='green', 
           alpha=0.7)
    
    ax.set_ylabel('Suitability Score')
    ax.set_title('Method Suitability for Different Problems')
    ax.set_xticks(x)
    ax.set_xticklabels(categories)
    ax.legend()
    ax.set_ylim(0, 1.2)
    ax.grid(True, axis='y', alpha=0.3)
    
    # 计算复杂度对比
    ax = axes[0, 1]
    n_values = np.logspace(1, 4, 50)
    fem_complexity = n_values**1.5 * 1e-8
    bem_complexity = n_values**3 * 1e-12
    fem_bem_complexity = n_values**2 * 1e-10
    
    ax.loglog(n_values, fem_complexity, 'b-', lw=2, label='FEM O(N^1.5)')
    ax.loglog(n_values, bem_complexity, 'r-', lw=2, label='BEM O(N^3)')
    ax.loglog(n_values, fem_bem_complexity, 'g-', lw=2, label='FEM-BEM O(N^2)')
    ax.set_xlabel('Degrees of Freedom (N)')
    ax.set_ylabel('Relative Computation Time')
    ax.set_title('Computational Complexity Comparison')
    ax.legend()
    ax.grid(True, which='both', ls='--', alpha=0.5)
    
    # 内存使用
    ax = axes[1, 0]
    fem_memory = n_values * 50  # 稀疏矩阵
    bem_memory = n_values**2    # 稠密矩阵
    fem_bem_memory = n_values * 50 + (n_values*0.1)**2
    
    ax.loglog(n_values, fem_memory, 'b-', lw=2, label='FEM')
    ax.loglog(n_values, bem_memory, 'r-', lw=2, label='BEM')
    ax.loglog(n_values, fem_bem_memory, 'g-', lw=2, label='FEM-BEM')
    ax.set_xlabel('Degrees of Freedom (N)')
    ax.set_ylabel('Memory Usage (elements)')
    ax.set_title('Memory Usage Comparison')
    ax.legend()
    ax.grid(True, which='both', ls='--', alpha=0.5)
    
    # 应用范围
    ax = axes[1, 1]
    applications = ['天线设计', '散射分析', '微波器件', 'EMC分析', '生物电磁']
    fem_applicability = [0.7, 0.5, 0.9, 0.8, 0.9]
    bem_applicability = [0.9, 0.9, 0.6, 0.7, 0.6]
    fem_bem_applicability = [0.95, 0.95, 0.85, 0.85, 0.85]
    
    y_pos = np.arange(len(applications))
    
    ax.barh(y_pos - 0.25, fem_applicability, 0.25, label='FEM', 
            color='blue', alpha=0.7)
    ax.barh(y_pos, bem_applicability, 0.25, label='BEM', 
            color='red', alpha=0.7)
    ax.barh(y_pos + 0.25, fem_bem_applicability, 0.25, label='FEM-BEM', 
            color='green', alpha=0.7)
    
    ax.set_xlabel('Applicability Score')
    ax.set_title('Application Domain Suitability')
    ax.set_yticks(y_pos)
    ax.set_yticklabels(applications)
    ax.legend()
    ax.set_xlim(0, 1.1)
    ax.grid(True, axis='x', alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fem_bem_summary.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存到 fem_bem_summary.png")
    plt.close()


def main():
    """主函数"""
    print("=" * 70)
    print("主题078: FEM-BEM耦合方法")
    print("有限元-边界元耦合方法用于电磁场开域问题仿真")
    print("=" * 70)
    
    # 运行所有演示
    demo_fem_solver()
    demo_bem_solver()
    demo_scattering()
    demo_radiation_pattern()
    demo_efficiency_comparison()
    demo_coupling_interface()
    demo_green_function()
    demo_fem_bem_summary()
    
    print("\n" + "=" * 70)
    print("所有演示完成!")
    print("=" * 70)


if __name__ == "__main__":
    main()

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

2026论文写作工具红黑榜：AI论文平台怎么选？一篇讲透：

2026年论文写作工具竞争白热化，红榜首选千笔AI、ThouPen、豆包，全面适配国内学术规范与格式要求；黑榜需警惕低质免费工具、无真实引用来源平台及过度依赖全文生成的工具。选择时

AtomGit开源社区

RabbitMQ在大数据领域的集群搭建指南

在大数据处理场景中，消息队列作为分布式系统的核心组件，需要支撑每秒数万级的消息吞吐量、毫秒级的延迟要求以及99.99%的可用性。RabbitMQ作为实现AMQP协议的开源消息中间件，其集群架构能有效解决单点故障、性能瓶颈等问题。本文聚焦RabbitMQ 3.10+版本，覆盖从基础概念到生产级集群部署的完整技术栈，包含架构设计、节点配置、数据同步、负载均衡、监控告警等核心内容。核心概念：解析Rabb

AtomGit开源社区

godot2D游戏教程系列二（22）

我是根据b站上的视频进行学习，并且总结写下笔记，然后做此分享。笔记非常详细以至于你可以通过查看文章进行快速学习并且制作出游戏出来。当然视频学习的链接我也放在下面了：视频学习：https://www.bilibili.com/video/BV1Nd25BNEA6?音频/美术资源包下载： https://pan.baidu.com/s/1mY05SbG3XFUpn-qiKxDRhw?pwd=b6bg。