悬臂梁的博客

8节点六面体单元的构造和数值积分(1)：完全积分

发表于 2024-07-21 更新于 2024-07-28 分类于 2.单元构造， 1.8节点六面体单元

前言

使用有限元方法对三维实体模型进行分析时，同样网格尺寸条件下，相比于四面体单元，六面体单元网格更少、效率和精度更高。

而不管是四面体单元还是六面体单元，都存在一阶和二阶之分，一阶单元精度不如二阶单元，但效率比二阶单元高。

在这篇文章和接下来几篇文章，我会介绍一下一阶六面体单元，即8节点六面体单元的构造及其积分方法，包括：

完全积分方法
选择减缩积分方法
附加9自由度的非协调模式积分方法
附加13自由度的非协调模式积分方法
减缩积分方法

单元构造

刚度矩阵

三维实体单元的单元刚度矩阵为： \[ \mathbf{k}_e = \int_{V_e}{\mathbf{B}^T \mathbf{DB}} \mathrm{d}V \tag{1} \]

其中，$\mathbf{D}$ 为弹性矩阵，由材料的弹性模量 $E$ 和泊松比 $\nu$ 确定， \[ \mathbf{D} = \frac{E(1 - \nu)}{(1 + \nu)(1 - 2 \nu)} \begin{bmatrix} 1 & \frac{\nu}{1 - \nu} & \frac{\nu}{1 - \nu} & 0 & 0 & 0 \\ \frac{\nu}{1 - \nu} & 1 & \frac{\nu}{1 - \nu} & 0 & 0 & 0 \\ \frac{\nu}{1 - \nu} & \frac{\nu}{1 - \nu} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1 - 2 \nu}{2(1 - \nu)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1 - 2 \nu}{2(1 - \nu)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1 - 2 \nu}{2(1 - \nu)} \\ \end{bmatrix} \tag{2} \]

$\mathbf{B}$ 为应变矩阵， \[ \mathbf{B} = \mathbf{LN} \tag{3} \]

其中，$\mathbf{L}$ 为微分算子， \[ \mathbf{L} = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \\ 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x} \end{bmatrix} \tag{4} \]

$\mathbf{N}$ 为单元形函数。在实际工程计算中，有限元网格往往是不规则的，直接计算形函数并不方便。因此，引入自然坐标 $(\xi, \eta, \zeta)$，以简化形函数的计算。 8节点六面体单元在总体（笛卡尔）坐标系与局部（自然）坐标系下的示意图如下所示。

节点按逆时针方向进行编号。

引入自然坐标后，自然坐标系下的形函数为： \[ N_i = \frac{1}{8}(1+\xi \xi_i)(1+\eta \eta_i)(1+\zeta \zeta_i) \tag{5} \]

将其展开，则为： \[ \begin{matrix} N_1 = \frac{1}{8}(1 - \xi)(1 - \eta)(1 - \zeta) \\ N_2 = \frac{1}{8}(1 + \xi)(1 - \eta)(1 - \zeta) \\ N_3 = \frac{1}{8}(1 + \xi)(1 + \eta)(1 - \zeta) \\ N_4 = \frac{1}{8}(1 - \xi)(1 + \eta)(1 - \zeta) \\ N_5 = \frac{1}{8}(1 - \xi)(1 - \eta)(1 + \zeta) \\ N_6 = \frac{1}{8}(1 + \xi)(1 - \eta)(1 + \zeta) \\ N_7 = \frac{1}{8}(1 + \xi)(1 + \eta)(1 + \zeta) \\ N_8 = \frac{1}{8}(1 - \xi)(1 + \eta)(1 + \zeta) \end{matrix} \tag{6} \]

上述形函数定义在自然坐标系 $(\xi, \eta, \zeta)$ 下，而应变矩阵需要的是形函数对总体坐标系 $(x, y, z)$ 求偏导，这就需要用到链式法则： \[ \begin{matrix} \frac{\partial N_i}{\partial \xi} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \xi} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial \xi} + \frac{\partial N_i}{\partial z} \frac{\partial z}{\partial \xi}\\ \frac{\partial N_i}{\partial \eta} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \eta} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial \eta} + \frac{\partial N_i}{\partial z} \frac{\partial z}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \zeta} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial \zeta} + \frac{\partial N_i}{\partial z} \frac{\partial z}{\partial \zeta} \tag{7} \end{matrix} \]

式(7)可写为矩阵形式： \[ \begin{Bmatrix} \partial N_i / \partial \xi \\ \partial N_i / \partial \eta \\ \partial N_i / \partial \zeta \end{Bmatrix} = \mathbf{J} \begin{Bmatrix} \partial N_i / \partial x \\ \partial N_i / \partial y \\ \partial N_i / \partial z \\ \end{Bmatrix} \tag{8} \] 其中，$\mathbf{J}$ 称为雅可比矩阵，定义为： \[ \mathbf{J} = \begin{bmatrix} \partial x / \partial \xi & \partial y / \partial \xi & \partial z / \partial \xi \\ \partial x / \partial \eta & \partial y / \partial \eta & \partial z / \partial \eta \\ \partial x / \partial \zeta & \partial y / \partial \zeta & \partial z / \partial \zeta \end{bmatrix} \tag{9} \]

所以，公式(8) 可改写为： \[ \begin{Bmatrix} \partial N_i / \partial x \\ \partial N_i / \partial y \\ \partial N_i / \partial z \\ \end{Bmatrix} = \mathbf{J}^{-1} \begin{Bmatrix} \partial N_i / \partial \xi \\ \partial N_i / \partial \eta \\ \partial N_i / \partial \zeta \end{Bmatrix} \tag{10} \]

公式(10) 可用于计算应变矩阵 $\mathbf{B}$。

计算出应变矩阵 $\mathbf{B}$ 后，公式(1) 可改写为： \[ \mathbf{k}_e = \int_{V_e} \mathbf{B}^T \mathbf{DB} \mathrm{d}V = \int_{-1}^{+1} \int_{-1}^{+1} \int_{-1}^{+1} \mathbf{B}^T \mathbf{DB} \mathrm{det}[\mathbf{J}] \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta \tag{11} \]

直接计算公式(11)是非常困难的，因此，一般使用数值方法对其进行计算，最常用的是高斯积分方法： \[ \mathbf{I} = \int_{-1}^{+1} \int_{-1}^{+1} \int_{-1}^{+1} f(\xi, \eta, \zeta) \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta = \sum_{i=1}^n \sum_{i=1}^m \sum_{i=1}^l w_i w_j w_k f(\xi_i, \eta_j, \eta_k) \tag{12} \]

质量矩阵

三维实体单元的质量矩阵为： \[ \mathbf{m}_e = \int_{V_e} \rho \mathbf{N}^T \mathbf{N} \mathrm{d}V \tag{13} \]

可参考式(11)，使用形函数进行计算。

载荷列阵

三维实体单元的载荷列阵为： \[ \mathbf{F}_b = \int_{V_e}\mathbf{N}^T \mathbf{f}_b \mathrm{d}V \tag{14} \] 和 \[ \mathbf{F}_s = \int_{S_e}\mathbf{N}^T \mathbf{f}_s \mathrm{d}S \tag{15} \] 其中， $\mathbf{F}_b$ 和 $\mathbf{F}_s$ 分别为体力和面力的等效节点载荷。

同样的，可以参考公式(11)，使用形函数进行计算。

完全积分

对于8节点六面体单元，如果使用完全积分方法，会产生剪切自锁现象，导致单元刚度增大，结构位移小于理论值。所以ANSYS和ABAQUS都没有集成完全积分方法（ABAQUS的C3D8单元使用的选择减缩积分方法，我会在下一篇文章中介绍）。

但是，完全积分方法比较基础，从其入门还是很有必要。

一维单元的高斯积分点和权系数如下图所示。

而三维六面体单元的高斯积分点和权系数，可在一维单元高斯积分点和权系数的基础上扩展而来，即在3个方向上均使用一维的高斯积分点和权系数。

所谓完全积分，指的是在3个方向上，都使用2个积分点进行积分，总共 2x2x2=8个积分点。

以下代码是完全积分方法的核心代码，截取自参考文献[1]：

其中，nip为8，der为形函数对自然坐标的偏导数，jac（前、后）分别为雅可比矩阵和雅可比矩阵的逆，deriv为形函数对总体坐标系的偏导数。截图中使用到的其它子例程，可以下载文献[1]的随书代码进行查看。需要注意的是，文献[1]中节点编号方式与本文所用的编号方式不同，而本文的节点编号方式是主流方式。

验证

参考文献[1]，自编有限元程序。为验证程序的正确性，将其与商业软件进行对比。

由于ANSYS和ABAQUS未集成完全积分方法，因此，使用Nastran进行验证。

Nastran默认也不会使用完全积分方法，需要修改PSOLID卡片，将IN设置为TWO，将ISOP设置为FULL。

1	PSOLID 1 1 TWO FULL

验证模型为长 1000mm 的悬臂梁，梁截面为边长 100mm 的正方形，一端固定，另一端承受 1000N 的集中力。材料为钢，弹性模量 200000MPa，泊松比 0.3，如下图所示：

静力学计算结果如下图所示：

可以看出，自编程序的结果（右）在z方向上与Nastran（左）完全一致，x、y方向有细微差别。这是由于在约束的处理上，自编代码使用了罚函数方法，而Nastran使用了删除行列的方法（文献[1]也是使用的删除行列的方法，若直接使用文献[1]的代码，其解过将与Nastran完全一致）。

去掉集中力，对模型进行模态分析，分别使用集中质量矩阵和一致质量矩阵，计算结果如下：

可以看出，无论是集中质量矩阵（上）还是一致质量矩阵（下），自编程序（右）与Nastran（左）计算出的模态频率均一致。

参考文献

[1] Smith I M, Griffiths D V, Margetts L. Programming the finite element method[M]. John Wiley & Sons, 2013.

[2] Liu G R, Quek S S. The finite element method: a practical course[M]. Butterworth-Heinemann, 2013.

[3] 王勖成. 有限单元法[M]. 清华大学出版社有限公司, 2003.

[4] ANSYS Mechanical APDL Theory Reference

特征问题的解法——向量迭代法(2)

发表于 2024-04-14 更新于 2024-07-28 分类于 3.数值计算， 3.特征问题

前言

上一篇文章中介绍了计算特征值的逆迭代算法，我们知道基本逆迭代算法收敛于 $\lambda_1$ 和 $\Phi_1$，基本正迭代算法收敛于 $\lambda_n$ 和 $\Phi_n$。现在假定我们已经计算出了其中一个特征值对，如 $(\lambda_k, \Phi_k)$，该怎么求解出其它的特征值对呢？

为了确保迭代过程不再收敛于 $(\lambda_k, \Phi_k)$，需要收缩 (Deflation)矩阵或迭代向量。

收缩

矩阵收缩

矩阵收缩的基本思想是，求一个正交矩阵 $\mathbf{H}$，使得 $\mathbf{HKH^{-1}}$ 为一个形如 \[ \left[ \begin{array}{c|cc} \lambda_1 & 0 & \cdots & 0 \\ \hline 0 & \\ \vdots & & \mathbf{K_1}\\ 0 \end{array} \tag{1} \right] \] 的矩阵。由于 $\mathbf{K}$ 和 $\mathbf{HKH^{-1}}$ 是相似的，故它们有相同的特征多项式。因此，若 $\mathbf{HKH^{-1}}$ 形如 (1)，则 \[ \rm{det} (\mathbf{K - \lambda I}) = \rm{det} (\mathbf{HKH^{-1} -\lambda I}) = \mathbf{(\lambda_1 - \lambda)} \rm{det} (\mathbf{K_1 - \lambda I}) \] 由此可得，$\mathbf{K}$ 的其余 $n-1$ 个特征值是 $\mathbf{K_1}$ 的特征值。

而一旦使用 $\mathbf{K_1}$ 算出了第二个所求特征对后，就可以对 $\mathbf{K_1}$ 重复该收缩过程，直到算出所有需要求解的特征值和特征向量。

应指出，矩阵 $\mathbf{H}$ 不是唯一的，因此可以使用各种方法构建一个适当的变换矩阵。另外，由于刚度矩阵 $\mathbf{K}$ 是带状的，所以，该变换不应破坏其带状结构。

向量收缩与Gram-Schmidt方法

为了得到其它特征对，我们也可以选择不收缩矩阵，而是收缩迭代向量。

向量收缩的基本思想是，逆迭代或正迭代过程中，为了使迭代向量收敛于所求的特征向量，该迭代向量应不与特征向量正交。反之，如果迭代向量正交于已计算出的特征向量，那么我们就消除了迭代向量收敛于这些特征向量的可能性，即，迭代向量将收敛于其它特征向量。

Gram-Schmidt方法是一个被广泛应用的向量正交化方法。该方法可用于求解广义特征值问题 $\mathbf{K \Phi = \lambda M \Phi}$。

假设我们已经利用逆迭代算法求解出特征向量 $\Phi_1,\Phi_2,\cdots,\Phi_m$，我们需要 $\mathbf{x_1}$ 与这些特征向量 $\mathbf{M}$ 正交化。

在Gram-Schmidt正交化中，特征向量 $\Phi_1,\Phi_2,\cdots,\Phi_m$ 的 $\mathbf{M}$ 正交向量 $\mathbf{\widetilde{x}_1}$ 的计算方法为： \[ \begin{align} \mathbf{\widetilde{x}_1} &= \mathbf{x}_1 - \sum_{i-1}^{m}{\alpha_i \Phi_i} \tag{2} \\ \alpha_i &= \mathbf{\Phi_i^T M x_1}; \quad i=1,\cdots,m \tag{3} \\ \end{align} \]

在逆迭代算法中，我们现在就可以使用 $\mathbf{\widetilde{x}_1}$ 作为初始迭代向量。然后，只要 $\mathbf{x_1^T M \Phi_{m+1}} \neq 0$，迭代向量就将收敛于 $\Phi_{m+1}$ 和 ${\lambda_{m+1}}$。

代码实现和测试

在上一篇文章中Python代码的基础上进行修改，实现向量收缩的代码。

代码实现和测试

import scipy as sp
import numpy as np


def inverse_iteration_method1(k: np.matrix, m: np.matrix):
    n = k.shape[0]
    eig_val = np.zeros(n)
    eig_vec = np.matrix(np.zeros((n, n)))

    eps = 1e-6
    for i in range(n):
        x = np.matrix(np.ones(n).reshape(n, 1))
        x_pre = x
        while True:
            for j in range(i):
                alpha = eig_vec[:, j].transpose() * m * x
                x -= alpha[0, 0] * eig_vec[:, j]
            x = np.linalg.solve(k, m * x)
            norm = x.transpose() * m * x
            x /= np.sqrt(norm)
            if np.all(abs(x - x_pre) < eps):
                eig_val[i] = 1/np.sqrt(norm[0, 0])
                eig_vec[:, i] = x
                break
            x_pre = x
    return eig_val, eig_vec


def inverse_iteration_method2(k: np.matrix, m: np.matrix):
    n = k.shape[0]
    eig_val = np.zeros(n)
    eig_vec = np.matrix(np.zeros((n, n)))

    eps = 1e-6
    rho_pre = 0.0
    for i in range(n):
        x = np.matrix(np.ones(n).reshape(n, 1))
        y = m * x
        while True:
            x = np.linalg.solve(k, y)
            for j in range(i):
                alpha = eig_vec[:, j].transpose() * m * x
                x -= alpha[0, 0] * eig_vec[:, j]
            y_bar = m * x
            rho = x.transpose() * y / (x.transpose() * y_bar)
            if abs(rho - rho_pre) / rho < eps:
                eig_val[i] = rho[0, 0]
                eig_vec[:, i] = x / np.sqrt(x.transpose() * y_bar)
                break
            y = y_bar / np.sqrt(x.transpose() * y_bar)
            rho_pre = rho
    return eig_val, eig_vec


def test():
    k = np.matrix(np.array([
        [5, -4, 1, 0],
        [-4, 6, -4, 1],
        [1, -4, 6, -4],
        [0, 1, -4, 5]
    ]))
    m = np.matrix(np.array([
        [2, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ]))
    n = k.shape[0]
    eig_vals, eig_vecs = sp.linalg.eig(k, m)[:2]
    eig_vals = eig_vals.real
    argsort = np.argsort(eig_vals)
    eig_vals = [eig_vals[argsort[i]] for i in range(n)]
    eig_vecs_copy = eig_vecs.copy()
    for i in range(n):
        eig_vecs[:, i] = eig_vecs_copy[:, argsort[i]]
    print('SciPy库：')
    print(eig_vals)
    print(eig_vecs)

    eig_vals, eig_vecs = inverse_iteration_method1(k, m)
    print('方法1：')
    print(f'特征值：\n{eig_vals}')
    print(f'特征向量：\n{eig_vecs}')

    eig_vals, eig_vecs = inverse_iteration_method2(k, m)
    print('方法2：')
    print(f'特征值：\n{eig_vals}')
    print(f'特征向量：\n{eig_vecs}')

程序输出如下：

SciPy库：
[0.09653732854936428, 1.3914654511583402, 4.373549554582955, 10.638447665709336]
[[-0.38576651  0.50215561 -0.55097093  0.11195372]
 [-0.61138815  0.14033458  0.52342732 -0.26606756]
 [-0.59120182 -0.55197537  0.02916661  0.75798869]
 [-0.35758794 -0.65074507 -0.64931054 -0.58491672]]
方法1：
特征值：
[ 0.09653733  1.39146545  4.37354955 10.63844769]
特征向量：
[[ 0.31262955  0.44526604  0.43866944 -0.10756245]
 [ 0.49547586  0.12443555 -0.41674032  0.25563076]
 [ 0.4791166  -0.48944215 -0.02322167 -0.72825456]
 [ 0.28979328 -0.5770219   0.51696616  0.56197132]]
方法2：
特征值：
[ 0.09653733  1.39146547  4.37354974 10.63844707]
特征向量：
[[ 0.31263514  0.44529746  0.43861285 -0.10764853]
 [ 0.49547613  0.12439816 -0.41670087  0.25571251]
 [ 0.47911154 -0.48944852 -0.023328   -0.72825001]
 [ 0.28978869 -0.57698412  0.51712097  0.56186986]]

可以看出，程序正确计算出了矩阵的全部特征值和特征向量。（特征向量因归一化方式不同，存在数字上的差异。若以相同的规则缩放，结果将相同）

参考文献

[1]. Bathe K J. Finite element procedures[M]. Klaus-Jurgen Bathe, 2006.

[2]. 史蒂文・J.利昂著STEVENJ.LEON.线性代数[M].机械工业出版社,2015.

[3]. 朱伯芳.有限单元法原理与应用-第3版[M].中国水利水电出版社,2009.

特征问题的解法——向量迭代法(1)

发表于 2024-04-06 更新于 2024-07-28 分类于 3.数值计算， 3.特征问题

前言

在工程中，我们经常需要对模型进行模态分析，以确定模型的模态频率和振型。

模态分析本质上是求解特征值问题 \[ {\mathbf K \phi} = \lambda {\mathbf M \phi} \tag{1} \] 特别是求解最小特征值 $\lambda_1, \dots, \lambda_p$ 及其特征向量 $\phi_1, \dots, \phi_p$ 。

常用的特征值求解算法可分为四类，对应于求解算法利用的基本性质。

第一类是向量迭代法(Vector Iteration Method)，又称为幂法(Power Iteration Method)，利用的基本性质是 \[ {\mathbf K \phi_i} = \lambda_i{\mathbf M \phi_i} \tag{2} \]

第二类为变换法，利用的基本性质是 \[ {\mathbf{\Phi^T K \Phi}} = {\mathbf{\Lambda}} \tag{3} \] \[ {\mathbf{\Phi^T M \Phi}} = {\mathbf{I}} \tag{4} \] 其中， ${\rm \Phi} = \left[ {\rm \phi_i, \dots, \phi_n} \right]$ 和 ${\rm \Lambda} = {\rm diag(\lambda_i)}, i=1,\dots,n$ 。

第三类为多项式迭代法，利用的基本性质是 \[ p(\lambda_i) = 0 \tag{5} \] 其中， \[ p(\lambda) = {\rm{det}(\mathbf{K} - \lambda \mathbf{M})} \tag{6} \]

第四类求解算法利用特征多项式的 Sturm 序列性质 \[ p(\lambda) = {\rm{det}(\mathbf{K} - \lambda \mathbf{M})} \tag{7} \] \[ p^{(r)}(\lambda^{(r)}) = {\rm{det}(\mathbf{K}^{(r)} - \lambda^{(r)}\mathbf{M}^{(r)})}; \quad r=1,\dots,n-1 \tag{8} \] 其中， $p^{(r)}(\lambda^{(r)})$ 是对应于${\rm K \phi = \lambda M \phi}$ 的第 $r$ 个相伴约束问题的特征多项式。

在四类求解方法中，每一类都提出了许多算法。

在学习具体算法前，需要认识到，求解算法应具有迭代性质，因为求解特征问题 ${\rm K \phi = \lambda M \phi}$ 就等价于计算多项式 $p(\lambda)$ 的根，它的阶等于 ${\rm K}$ 和 ${\rm M}$ 的阶。

接下来开始介绍第一类方法：向量迭代法。而在介绍向量迭代法之前，先尝试理解其基本思想。

向量迭代法

基本思想

假设向量 ${\rm A}$ 有 $n$ 个线性无关的特征向量 $\mathit{x_1, \dots, x_n}$ ，且相应的特征值满足 \[ |\lambda_1| > |\lambda_2| \geq \dots \geq |\lambda_n| \]

给定一个任意向量 $\mathit{v}_0$ ，假设 \[ \mathit{v}_0 = \alpha_1 \mathit{x}_1 + \dots + \alpha_n \mathit{x}_n \]

将 ${\rm A}$ 作用于该向量，有 \[ \mathbf{A} \mathit{v}_0 = \alpha_1 \lambda_1 \mathit{x}_1 + \alpha_2 \lambda_2 \mathit{x}_2 + \dots + \alpha_n \lambda_n \mathit{x}_n \\ \mathbf{A}^2 \mathit{v}_0 = \alpha_1 \lambda_1^2 \mathit{x}_1 + \alpha_2 \lambda_2^2 \mathit{x}_2 + \dots + \alpha_n \lambda_n^2 \mathit{x}_n \]

且一般地， \[ \mathbf{A}^k \mathit{v}_0 = \alpha_1 \lambda_1^k \mathit{x}_1 + \alpha_2 \lambda_2^k \mathit{x}_2 + \dots + \alpha_n \lambda_n^k \mathit{x}_n \tag{9} \]

若定义 \[ \mathit{v}_k = \mathbf{A}^k \mathit{v}_0, \quad k=1, 2, \dots \]

则 \[ \frac{1}{\lambda_1^k} \mathit{v}_k = \alpha_1 \mathit{x}_1 + \alpha_2 (\frac{\lambda_2}{\lambda_1})^k \mathit{x}_2 + \dots + \alpha_2 (\frac{\lambda_n}{\lambda_1})^k \mathit{x}_n \tag{10} \]

由于 \[ \left| \frac{\lambda_i}{\lambda_1} \right| < 1, \quad, i=2,3,\dots,n \]

由此得到： \[ \frac{1}{\lambda_1^k} \mathit{v}_k \rightarrow \alpha_1 \mathit{x}_1, \quad k \rightarrow \infty \]

因此，若 $\alpha_1 \neq 0$ ，则序列 $\{ (1/\lambda_1^k) \mathit{v}_k \}$ 收敛到 $\mathbf{A}$ 的特征向量 $\alpha_1 \mathit{x}_1$。

当然，由于 $\lambda_1$ 是未知的，所以无法计算 $\{ (1/\lambda_1^k) \mathit{v}_k \}$ 。但好在不需要将序列 $\{ \mathit{v}_k \}$ 用 $1/\lambda_1^k$ 进行缩放。

这就是正迭代法的思想，能计算出最大的特征值。

而模态分析需要计算的是最小特征值，需要使用的是逆迭代法，即需要将 $\mathbf{A}^{-1}$ 作用在向量上。

逆迭代和正迭代

向量迭代法所考虑的基本关系式是 \[ {\mathbf K \phi} = \lambda {\mathbf M \phi} \tag{1} \]

选择 $\phi$ 的一个向量 $\mathbf{x}_1$ ，对 $\lambda$ 设定一个值，令 $\lambda = 1$ 。于是，可以计算公式 (1) 的右手边，即可以计算 \[ \mathbf{R}_1 = (1)\mathbf{M x_1} \tag{11} \] 由于 $\mathbf{x}_1$ 是任意向量，一般不满足 $\mathbf{K x_1 = R_1}$ 。考虑静平衡方程 $$

; $$ 其中， $\mathbf{x_2}$ 是位移解，对应于作用力 $\mathbf{R}$ 。由于我们是使用迭代法求解特征向量，所以我们可以直观地认为，$\mathbf{x_2}$ 是比 $\mathbf{x_1}$ 更好的近似特征向量。通过反复进行此迭代，可以得到越来越好的近似特征向量。

上述过程即为逆迭代的基础。

而在正迭代中，迭代过程是反向的，即在第一步中，计算 $\mathbf{R_1 = K x_1}$ ，然后通过求解 $\mathbf{M x_2 = R_1}$ ，得到改进的近似特征向量 $\mathbf{x_2}$ 。

基本算法描述

首先考虑逆迭代法中所用的基本方程。

在求解中，假设初始迭代向量为 $\mathbf{x}_1$，在每个迭代步骤 $k = 1, 2,\dots$ 时，计算 \[ \begin{align} \mathbf{K \bar{x}}_{k+1} &= \mathbf{M x}_k \tag{13} \\ \mathbf{x}_{k+1} &= \frac{\mathbf{\bar{x}}_{k+1}}{(\mathbf{\bar{x}}_{k+1}^T \mathbf{M} \mathbf{\bar{x}}_{k+1})^{1/2}} \tag{14} \end{align} \] 其中，只要 $\mathbf{x}_1$ 和 $\phi_i$ 不与 $\mathbf{M}$ 正交，即 $\mathbf{x_1^T M \phi_1} \neq 0$ ，就有 \[ 当 k \rightarrow \infty, \quad \mathbf{x}_{k+1} \rightarrow \phi_1 \]

迭代法的基本步骤是求解方程 (13)，每求解一次，我们就能够得到比前一次迭代向量 $\mathbf{x}_{k}$ 更接近特征向量的新向量 $\mathbf{x}_{k+1}$ 。式 (14) 的计算对新的迭代向量进行缩放，使其与 $\mathbf{M}$ 的加权长度为 $1$ ，即令 $\mathbf{x}_{k+1}$ 满足质量正交条件 \[ \mathbf{x_{k+1}^{T} M x_{k+1}} = 1 \tag{15} \] 当然，也可以将向量 $\mathbf{x}_{k+1}$ 缩放到其他长度。

更有效的算法

式 (13)、(14) 是基本的逆迭代算法，但在实际实现中，以下方法更加有效。

假设 $\mathbf{y}_1 = \mathbf{M x_1}$ ，计算 $k = 1, 2, \dots$ ， \[ \begin{align} \mathbf{K \bar{x}}_{k+1} &= \mathbf{y}_k \tag{16} \\ \mathbf{\bar{y}}_{k+1} &= \mathbf{M \bar{x}}_{k+1} \tag{17} \\ \rho(\mathbf{\bar{x}}_{k+1}) &= \frac{\mathbf{\bar{x}}_{k+1}^T \mathbf{y}_k}{\mathbf{\bar{x}}_{k+1}^T \mathbf{\bar{y}}_{k+1}} \tag{18} \\ \mathbf{y}_{k+1} &= \frac{\mathbf{\bar{y}}_{k+1}}{(\mathbf{\bar{x}_{k+1}^T \mathbf{\bar{y}}_{k+1}})^{1/2}} \tag{19} \end{align} \] 其中，只要 $\mathbf{y}_1^T \phi_1 \neq 0$ ，有 \[ 当 k \rightarrow \infty 时，\mathbf{y}_{k+1} \rightarrow \mathbf{M \phi_1} 和 \phi(\mathbf{\bar{x}}_{k+1}) \rightarrow \lambda_1 \] 本算法对 $\mathbf{y}_k$ 进行迭代而非对 $\mathbf{x}_k$ 进行迭代。式 (18) 中，我们得到了特征值 $\lambda_1$ 的近似值，其以瑞利商 $\rho(\mathbf{\bar{x}}_{k+1})$ 的形式给出。 $\lambda_1$ 的近似值可以用于评估迭代的收敛情况。若当前的 $\lambda_1$ 的近似值用 $\lambda_1^{(k+1)}$ 表示，即 $\lambda_1^{(k+1)} = \rho(\mathbf{\bar{x}}_{k+1})$，那么我们通过式 (20) 判断迭代是否收敛： \[ \frac{\left| \lambda_1^{(k+1)} - \lambda_1^{(k)} \right|}{\lambda_1^{(k+1)}} \leq tol \tag{20} \] 当特征值 $\lambda_1$ 需要 $2s$ 的精度时，$tol$ 应该是 $10^{-2s}$ 或更小。

代码实现和测试

实现

使用 Python 语言对基本算法和更高效的算法进行实现，代码如下。

import numpy as np


def inverse_iteration_method1(k: np.ndarray, m: np.ndarray):
    n = k.shape[0]
    x = np.ones((n, 1))
    x_pre = x
    i_iter = 0
    while True:
        i_iter += 1
        x = np.linalg.solve(k, m @ x)
        norm = x.transpose() @ m @ x
        x /= np.sqrt(norm)
        if np.all(abs(x - x_pre) < 1e-6):
            break
        x_pre = x
    print(f'方法1：迭代{i_iter}次')
    print(f'特征值：\n{1/np.sqrt(norm)}')
    print(f'特征向量：\n{x}')


def inverse_iteration_method2(k: np.ndarray, m: np.ndarray):
    n = k.shape[0]
    x = np.ones((n, 1))
    y = m @ x
    rho_pre = 0.0
    n_iter = 0
    while True:
        n_iter += 1
        x = np.linalg.solve(k, y)
        y_bar = m @ x
        rho = x.transpose() @ y / (x.transpose() @ y_bar)
        if abs(rho - rho_pre) / rho < 1e-6:
            break
        y = y_bar / np.sqrt(x.transpose() @ y_bar)
        rho_pre = rho
    print(f'方法2：迭代{n_iter}次')
    print(f'特征值：\n{rho}')
    print(f'特征向量：\n{x / np.sqrt(x.transpose() @ y_bar)}')

测试

使用文献 [1] 中的例子进行测试。

def test1():
    k = np.array([
        [2, -1, 0, 0],
        [-1, 2, -1, 0],
        [0, -1, 2, -1],
        [0, 0, -1, 1]
    ])
    m = np.array([
        [0, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 1]
    ])
    inverse_iteration_method1(k, m)


def test2():
    k = np.array([
        [2, -1, 0, 0],
        [-1, 2, -1, 0],
        [0, -1, 2, -1],
        [0, 0, -1, 1]
    ])
    m = np.array([
        [0, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 1]
    ])
    inverse_iteration_method2(k, m)

程序输出如下。

方法1：迭代8次
特征值：
[[0.14644661]]
特征向量：
[[0.25000003]
 [0.50000006]
 [0.60355338]
 [0.70710669]]
方法2：迭代5次
特征值：
[[0.14644661]]
特征向量：
[[0.25000638]
 [0.50001275]
 [0.60355075]
 [0.70708874]]

最小特征值和特征向量的精确解为 $\lambda_1 = \frac{1}{2} - \frac{\sqrt{2}}{4}$ 和 $\phi_1 = [\frac{1}{4} \quad \frac{1}{2} \quad \frac{1+\sqrt{2}}{4} \quad \frac{\sqrt{2}}{2}]^T$，可以看出，两种逆迭代算法正确算出了最小特征值和特征向量，而后一种算法所需的迭代次数更少。

计算出最小特征值后，剩余特征值的计算方法将在下一篇文章中介绍。

参考文献

[1]. Bathe K J. Finite element procedures[M]. Klaus-Jurgen Bathe, 2006.

[2]. 史蒂文・J.利昂著STEVENJ.LEON.线性代数[M].机械工业出版社,2015.

线性方程组的解法——Cholesky分解

发表于 2024-04-06 更新于 2024-07-28 分类于 3.数值计算， 1.静力学方程， 1.直接解法

Cholesky分解

在科学和工程计算中，经常需要求解形如 $\mathbf{A} \mathit{x} = \mathit{b}$ 的线性方程组，其中 $\mathbf{A}$ 为 $n \times m$ 矩阵，称为系数矩阵， $\mathit{b}$ 为 $n$ 维列向量，称为右端向量， $\mathit{x}$ 为待求解的 $m$ 维列向量，称为解向量。

而科学和工程的实际计算中，经常遇到系数矩阵 $\mathbf{A}$ 为对称正定矩阵的情况。若 \[ \mathbf{A}=\begin{bmatrix} a_{11} \\ a_{21} & a_{22} & & 对称\\ a_{31} & a_{32} & a_{33} \\ \vdots & \vdots & \vdots & \ddots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{bmatrix} \] 为正定阵，则有如下三角阵 \[ \mathbf{L} = \begin{bmatrix} l_{11} \\ l_{21} & l_{22} & & \mathbf{0}\\ l_{31} & l_{32} & l_{33} \\ \vdots & \vdots & \vdots & \ddots \\ l_{n1} & l_{n2} & l_{n3} & \dots & l_{nn} \end{bmatrix} \\ \] 使 $\mathbf{A} = \mathbf{L \cdot L^T}$ 成立。若 $\mathbf{L}$ 的主对角线元素取正值，则这种分解是唯一的。

将矩阵关系式 $\mathbf{A} = \mathbf{L \cdot L^T}$ 直接展开，有 \[ \begin{align*} a_{11} &= l_{11}^{2} \\ a_{21} &= l_{21}l_{11},\quad a_{22} = l_{21}^{2}+l_{22}^{2}\\ a_{31} &= l_{31}l_{11},\quad a_{32} = l_{31}l_{21}+l_{32}l_{22},\quad a_{33}=l_{31}^{2}+l_{32}^{2}+l_{33}^{2}\\ \dots \end{align*}\\ \]

据此可逐行求出矩阵 $\mathbf{L}$ 的元素 $l_{11} \rightarrow l_{21} \rightarrow l_{22} \rightarrow l_{31} \rightarrow l_{32} \rightarrow \dots$，计算公式为 \[ \begin{cases} l_{ij} &= (a_{ij} - \sum\limits_{k=1}^{j-1}l_{ik}l_{jk}) / l_{jj}, \quad & j = 1, 2, \dots, i-1 \\ l_{ii} &= (a_{ii} - \sum\limits_{k=1}^{i-1}l_{ik}^2)^\frac{1}{2}, \quad & i = 1, 2, \dots, n \\ \end{cases} \\ \]

基于矩阵分解式 $\mathbf{A} = \mathbf{L \cdot L^T}$，对称正定方程组 $\mathbf{A} \mathit{x} = \mathit{b}$ 可归结为两个三角方程组 $\mathbf{L} \mathit{y} = \mathit{b}$ 和 $\mathbf{L}^T \mathit{x} = \mathit{y}$ 来求解。

由 $\mathbf{L} \mathit{y} = \mathit{b}$ 即 \[ \begin{cases} l_{11}y_{1} &= b_1 \\ l_{21}y_{2} + l_{22}y_{2} &= b_2 \\ \dots \dots \dots \\ l_{n1}y_{1} + l_{n2}y_{2} + \dots + l_{nn}y_{n} &= b_n \end{cases} \\ \] 可顺序计算出 $y_1 \rightarrow y_2 \rightarrow \dots \rightarrow y_n$ ： \[ y_i = (b_i - \sum\limits_{k=1}^{i-1}l_{ik}y_{k})/l_{ii}, \quad i = 1, 2, \dots,n \\ \]

而由 $\mathbf{L}^T \mathit{x} = \mathit{y}$ 即 \[ \begin{cases} \begin{alignat*}{2} l_{11}x_1 + l_{21}x_2 + \dots + l_{n1}x_n &= y_1 \\ l_{22}x_2 + \dots + l_{n2}x_n &= y_2 \\ \dots \dots \dots \\ l_{nn}x_n &= y_n \\ \end{alignat*} \end{cases} \\ \] 可逆序求得 $x_n \rightarrow x_{n-1} \rightarrow \dots \rightarrow x_1$： \[ x_i = (y_i - \sum\limits_{k=i+1}^nl_{ki}x_{k})/l_{ii}, \quad i = n, n-1, \dots, 1 \\ \]

由于矩阵分解时公式含有开方运算，所以该算法称为平方根法，又叫Cholesky分解法。

代码实现（Fortran版）

根据上述公式，编写程序即可对方程进行求解：

subroutine cholesky_full(n, a, y)
    implicit none
    
    integer, intent(in) :: n
    real, intent(inout) :: a(n, n), y(n)
    
    integer :: i, j, k
    real :: temp
    
    ! 分解矩阵，生成下三角阵L
    ! 工程问题中的很多矩阵非常庞大，所以，计算过程中的数据应该直接存放在原始数组a中，
    ! 而不是新创建一个数组
    do i = 1, n
        ! 公式中，j的取值范围为1到j-1，此处换成1到j，可以将分解式统一起来，省去一次判断。
        ! 因为j=i时，j循环虽然会执行错误的操作、生成错误的a(i,j)结果，
        ! 但a(i,j)马上就会被最外层的i循环生成的正确数据替换
        do j = 1, i
            temp = a(i, j)
            do k = 1, j-1
                temp = temp - a(i, k) * a(j, k)
            end do
            a(i, j) = temp / a(j, j)
            ! a(j, i) = 0.  ! 对角线上方d的元素赋0，可有可无
        end do
        a(i, i) = sqrt(temp)
    end do
    
    ! 根据Ly=b求解出y
    do i = 1, n
        temp = y(i)
        do j = 1, i-1
            temp = temp - a(i, j) * y(j)
        end do
        y(i) = temp / a(i, i)
    end do
    
    ! 求解出x
    do i = n, 1, -1
        temp = y(i) / a(i, i)
        y(i) = temp
        ! 公式中k的范围为i+1到n，此处为1到i-1，因为下方a(i,k)的下标和公式中交换了顺序
        do k = 1, i-1
            y(k) = y(k) - temp * a(i, k)
        end do
    end do
end subroutine

以上代码的Cholesky分解部分与前文公式基本上一致，很好理解，但引入了一个临时变量temp，用于存储数据。而如果我们将j、k两层循环交换一下位置，再稍微调整一下循环计数器的取值范围，就可以不借助临时变量直接完成分解操作。代码如下：

do i = 1, n
    do k = 1, i - 1
        a(i, k) = a(i, k) / a(k, k)
        do j = k + 1, i
            a(i, j) = a(i, j) - a(i, k) * a(j, k)
        end do
    end do
    a(i, i) = sqrt(a(i, i))
end do

参考文献

[1].王能超. 高等学校教材, 数值分析简明教程, （第2版）[M]. 2003.

[2].吴建平, 王正华, 李晓梅. 稀疏线性方程组的高效求解与并行计算[M]. 湖南科学技术出版社, 2004.

各种命令的笔记

发表于 2023-11-19 更新于 2024-07-28 分类于 4.其它

Linux

解压缩

zip格式：

1	unzip filename.zip

tar格式：

1	tar -xvf filename.tar.xxx # 其中，xxx是gz、bz2、xz、z等后缀（1.15版之后可自动识别压缩格式）

磁盘

查看磁盘空间：

df -hl

Linux子系统（WSL）

移除子系统中的Windows环境变量（仅当次起效）:

1	PATH=$(echo "$PATH" \| sed -e 's/:\/mnt.*//g') # strip out problematic Windows %PATH% imported var

Hexo

新建文章

1	hexo new [layout] <title>

清理缓存和生成的静态文件；生成静态文件；启动服务，端口为5000

1	hexo clean; hexo g; hexo s -p 5000

OCC+VTK显示CAD模型

发表于 2023-10-15 更新于 2024-07-28 分类于 1.前后处理

前言

早些时候，这篇文章已经被我发布在了知乎上：“VTK+OCC显示CAD模型”，今天将文章重新编辑，发布在Github上，主要是为了测试Hexo的用法。

测试环境

系统：Win11
IDE：Visual Studio Community 2019
VTK：VTK 9.2
OCC：OCCT 7.6

正文

VTK是一款十分优秀的可视化套件，开源且功能强大，基本上可以满足有限元领域的全部可视化需求。遗憾的是，VTK不支持CAD模型（如igs、stp格式的模型）的显示。

在网上搜索后可以发现，在不花钱的情况下，想要显示和处理CAD模型，基本上都得使用OpenCasCade，即OCC。OCC有自己的可视化系统，也可以集成在Qt中。但对我而已，OCC自己的可视化系统还是太复杂了。

好在OCC在6.8版本开发了VIS（VTK Integration Services）功能，之后的版本就可以使用VTK进行模型的可视化了。

为了使用VIS功能，编译OCC的时候需要选择USE_VTK的选项，编译完成后，将生成TKIVtk、TKIVtkDraw的动态库和静态库。如果编译路径下有这两个库，说明VIS的功能是编译成功了。

编写一个最小案例看看显示效果。

CMakeLists.txt文件：

cmake_minimum_required(VERSION 3.15)

project(occvtk LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

find_package(VTK REQUIRED)
find_package(OpenCASCADE REQUIRED)
include_directories(${OpenCASCADE_INCLUDE_DIR})
link_directories(${OpenCASCADE_LIBRARY_DIR})

add_executable(test test.cpp)
target_link_libraries(test
    ${VTK_LIBRARIES}
    ${OpenCASCADE_LIBRARIES}
)
vtk_module_autoinit(
    TARGETS test
    MODULES ${VTK_LIBRARIES}
)

在此CMakeLists文件中，没有写死VTK和OCC的路径，而是使用find_package命令查找。如果没有找到，可以在命令行执行cmake命令时使用-D参数将相关路径传入:

1	cmake .. -DVTK_DIR:PATH=/home/me/vtk_build -DVTK_DIR -DOpenCASCADE_DIR:PATH=/home/me/occ_build

test.cpp文件:

#include <STEPControl_Reader.hxx>
#include <Standard_Integer.hxx>
#include <TopoDS_Shape.hxx>
#include <IFSelect_ReturnStatus.hxx>
#include <IFSelect_PrintCount.hxx>
#include <IVtkTools_ShapeDataSource.hxx>
#include <IVtkOCC_ShapeMesher.hxx>
#include <IVtkTools_DisplayModeFilter.hxx>
#include <vtkType.h>
#include <vtkAutoInit.h>
#include <vtkRenderWindow.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPolyDataMapper.h>
#include <vtkInteractorStyleTrackballCamera.h>

VTK_MODULE_INIT(vtkRenderingOpenGL2)
VTK_MODULE_INIT(vtkInteractionStyle)

int main()
{
    STEPControl_Reader reader;
    IFSelect_ReturnStatus stat = reader.ReadFile("assembly_solid.stp");
    IFSelect_PrintCount mode = IFSelect_CountByItem;
    Standard_Integer NbRoots = reader.NbRootsForTransfer();
    Standard_Integer num = reader.TransferRoots();
    Standard_Integer NbTrans = reader.TransferRoots();
    TopoDS_Shape result = reader.OneShape();
    // TopoDS_Shape shape = reader.Shape();

    vtkNew<IVtkTools_ShapeDataSource> occSource;
    //occSource->SetShape(new IVtkOCC_Shape(shape));
    occSource->SetShape(new IVtkOCC_Shape(result));

    vtkNew<IVtkTools_DisplayModeFilter> filter;
    filter->AddInputConnection(occSource->GetOutputPort());
    filter->SetDisplayMode(DM_Shading);

    vtkNew<vtkPolyDataMapper> mapper;
    mapper->SetInputConnection(filter->GetOutputPort());

    vtkNew<vtkActor> actor;
    actor->SetMapper(mapper);

    vtkNew<vtkRenderer> ren;
    ren->AddActor(actor);

    vtkNew<vtkRenderWindow> renWin;
    renWin->AddRenderer(ren);
    renWin->SetSize(960, 800);

    vtkNew<vtkInteractorStyleTrackballCamera> istyle;
    vtkNew<vtkRenderWindowInteractor> iren;

    iren->SetRenderWindow(renWin);
    iren->SetInteractorStyle(istyle);

    renWin->Render();
    iren->Start();

    return 0;
}

此测试代码将模型文件“assembly_solid.stp”写死在了源代码中，编译完成后，需要确保可执行文件可以找到模型文件。

显示效果如下：