X.2 销接结构2
刚度矩阵的推导：续
目前我们讨论的是两个杆件构成的模型。
当我们拓展到多个杆件时：
我们可以对单个杆件做分析，以下图为例。当前杆件的节点在全局坐标系将是第i、第j个节点，在局部坐标系将是第1和第2个节点。
我们可以参照在两杆桁架下的操作，把单元刚度矩阵“拼接”得到全局刚度矩阵，抓住节点对应位移的特征：
拼接得到的全局刚度矩阵如下：
在有限元模型中，节点一般按一定顺序编号，相邻单元的节点编号往往相对接近。
当把单元刚度矩阵拼装进全局刚度矩阵时，刚度的非零项只会出现在与这个单元的两个端点节点相关的行列中，要和端点的位移相对应，如上图所示。
这样，每个单元刚度矩阵只会影响这个单元所连接的几个节点对应的有限个自由度（也可称为每个自由度只和有限数量的自由度耦合）。
在二维/三维有限元网格中，单个节点的连接单元是有限的，这是的全局刚度矩阵中，非零元素相对于矩阵总大小而言占比例可能较小，因此刚度矩阵具有稀疏性；
此外，单个节点连接的几个有限自由度，在编号上通常彼此接近，从而在整体刚度矩阵中，形成集中于对角线附近、密集的非零块。
随着节点编号增加、结构空间分布的延展，整体刚度矩阵的非零项就呈现带状分布了。
例如，对上述具有10个节点的桁架结构，对节点1，它只和节点2、3、4有连接，因此这几个节点对应的列，刚度\(\times\)位移的值不为0。
讨论对应节点1 x方向位移的行：
\[0=K_{11}U_1+K_{12}V_1+K_{13}U_3+K_{14}V_2+... \]这里，\(K_{12}\)表示当节点1在y方向位移\(V_1\)变化时，对节点1在x方向反力产生的影响刚度；\(K_{13}\)表示当节点2在x方向位移\(U_2\)（第3个自由度）变化时，在节点1的x方向位移（第1个自由度）上引起的反力（刚度影响）。
意思是，假设所有自由度中，只有第3个自由度（节点2在x方向位移）向其指向位移了1个单位长度，其他自由度位移为0，此时第1个自由度（节点1在x方向位移）感受到反力的大小。
依此类推，把节点4的y方向位移变化，在节点1在x方向位移引起的反力的刚度写出来之后，关于其他节点的位移就为0了，因为没有直接连接。
这样，更多自由度的全局刚度矩阵呈现带状的原因就很清楚了：当前节点能连接的其他节点、能互相影响的自由度数量有限。
我们看到上面的图，定义了一条“天际线”（skyline，喜欢我城市天际线吗？），天际线以内的区域，包含天际线本身，包括了所有非零元素，天际线以外的区域为0，这利用了刚度矩阵的对称性，即矩阵中关于主对角线对称位置的元素相等。
从物理上考虑，根据功的互等定理，节点i的单位位移引起的节点j的力，等于节点j的单位位移引起的节点i的力。
或者说，刚度矩阵的本质是关于能量的二阶导数矩阵，而能量对位移的二阶偏导是对称的：
\[\begin{equation} \frac{\partial^2U}{\partial u_i\partial u_j}=\frac{\partial^2U}{\partial u_j\partial u_i} \end{equation} \]刚度矩阵代码化
存储
刚度矩阵是对称的，且存在变带宽结构，即每一行非零元素的数量（带宽） 不同。
我们引入一维存储方法，只需存储矩阵下/上三角的非零部分与必要位置。
注意到这里有两个列表。第一个列表是直接存储刚度数值的，第二个列表是ID数组，用于记录矩阵中每一行对角元素（或其余自定义的关键位置元素）在一维数组中的存放位置。通过这种标记，可以完整重建二维矩阵。
import numpy as np
K = np.array([
[10, 2, 0, 0, 0], # 第1行只与第2行有非零耦合
[ 2, 15, 5, 0, 0], # 第2行与第1、3行有联系
[ 0, 5, 20, 4, 0], # 第3行与第2、4行有联系
[ 0, 0, 4, 25, 3], # 第4行与第3、5行有联系
[ 0, 0, 0, 3, 30], # 第5行只与第4行有联系
], dtype=float)
# -----------------------------------------------
# 第一步：确定每行的非零带宽
# numpy中的array.shape返回一个表示数组维度大小的元组，这里返回(m, n)，表示行数和列数
n = K.shape[0]
band_width_per_row = []
for i in range(n): # 行遍历
# 找到该行中从左侧到对角线间的第一个非零元素位置（包括对角线）
# 因为是下三角，对j的搜索只需要从0到i即可
j_start = 0
for j in range(i, -1, -1):
# 列遍历，反向搜索，尝试找最左侧的非零元素
# 注意这里步长是-1
if K[i, j] != 0:
j_start = j
else:
# 一旦遇到0，则再往前都是0，继续循环直到找到非0或到达行首
pass
# 该行的带宽= (i - j_start + 1)
bandwidth = i - j_start + 1
band_width_per_row.append(bandwidth)
# 对于这里的矩阵，band_width_per_row是 [1,2,2,2,2]
# -----------------------------------------------
# 第二步：构建ID数组
# 创建一个全零数组，长度为n
ID = np.zeros(n, dtype=int)
sum_val = 0
for i in range(n):
sum_val += band_width_per_row[i]
# 累加：每行带宽是 [1,2,2,2,2]，本次累加的结果是当前行的对角线位置元素在一维数组中的真实位置（不是索引位置），如第1行是1、第2行是3、第3行是5
ID[i] = sum_val
# 现在 ID[i]-1 就是对角线元素在一维数组的位置(因为索引从0开始)
# 整个一维数组的长度 = sum(band_width_per_row)
total_length = sum(band_width_per_row)
K_1D = np.zeros(total_length, dtype=float)
# -----------------------------------------------
# 第三步：将K的下三角非零元素拷贝到K_1D中
for i in range(n):
# i行的带宽
bw = band_width_per_row[i]
# 对角元素在1D中的位置
diag_pos = ID[i] - 1
# 这一行在1D中对应的起始存储区间是 [diag_pos - (bw-1)] 到 diag_pos
start_pos = diag_pos - (bw - 1)
# 从K中提取下三角区（含对角）元素
# 行i中有效的列范围是 [i - (bw-1), i]
col_start = i - (bw - 1)
col_end = i
# 将这些元素复制到K_1D中
idx = start_pos # 起始存储索引
for col in range(col_start, col_end + 1):
K_1D[idx] = K[i, col] # 对当前行，提取有效列中数据
idx += 1 # 索引+1，往右边下一个非零元素走
# 至此，K_1D和ID数组就能描述原来的带状下三角刚度矩阵。
print("原二维刚度矩阵K:\n", K)
print("带宽信息:", band_width_per_row)
print("ID数组:", ID)
print("一维数组K_1D:", K_1D)
求解矩阵方程
核心方程分析
在结构有限元分析中，核心方程是：
\[\begin{equation}[K]\{\Delta\}= \{F\}\end{equation} \]把已知的外力或分布载荷转化为节点力F之后，可形成上述方程组。
但最初建立时往往是\([K]\{\Delta\}=\{F\}+\{R\}\)。其中R是附加支撑条件（约束）产生的反力。
由于约束条件使得某些位移为0，实际分析时将方程写成核心形式，解出位移\(\Delta\)后回代求R即可。
下一步是施加边界条件。
在刚度矩阵中，对应受约束的自由度的刚度可调成很大的值（如\(10^{50}\)），确保该自由度的最终位移为0。
施加后的方程变成：
\[\begin{equation} [\underline{K}]\{\Delta\}=\{F\} \end{equation} \]要求解的形式变成：
\[\begin{equation} \begin{Bmatrix}\Delta\end{Bmatrix}=\begin{bmatrix}\underline{K}\end{bmatrix}^{-1}\begin{Bmatrix}F\end{Bmatrix} \end{equation} \]求解方法
矩阵求逆
直接计算方程4，但计算过程昂贵、数值上不稳定、且存储需要庞大的内存，因此并不实际使用。
高斯消元
举例：
给定 \([K]\) 和 \(\{F\}\) ，假设 \([K]\) 是 3x3 矩阵，写成增广矩阵 \([K|\{F\}]\)形式：
\[\begin{equation} \left[\begin{array}{ccc|c} k_{11} & k_{12} & k_{13} & f_1 \\ k_{21} & k_{22} & k_{23} & f_2 \\ k_{31} & k_{32} & k_{33} & f_3 \\ \end{array}\right] \end{equation} \]用 \(k_{11}\) 对第二、三行进行消元，使 \(k_{21}\) 和 \(k_{31}\) 变为零；
接着对第二、三行使用第二行主元消元，使 (3,2) 元素变为零，从而最终形成上三角结构。
一旦得到上三角系统 \([U]\{\Delta\} = \{F'\}\)（\(\{F'\}\)是变换后的载荷向量），从末行开始逐个向上回代求解 \(\{\Delta\}\)。
LU分解
代码：
import numpy as np
def lu_decomposition(A):
"""
对给定的方阵 A 进行 LU 分解，分解为下三角矩阵 L 和上三角矩阵 U，使得 A = L * U。
参数：
A (np.ndarray): 要分解的方阵。
返回：
L (np.ndarray): 下三角矩阵，主对角线元素为 1。
U (np.ndarray): 上三角矩阵。
"""
n = A.shape[0] # 获取矩阵的维度（假设 A 为方阵）
L = np.zeros_like(A, dtype=np.float64) # 初始化 L 为与 A 同型的零矩阵，数据类型为 float64
U = np.zeros_like(A, dtype=np.float64) # 初始化 U 为与 A 同型的零矩阵，数据类型为 float64
for i in range(n):
# 计算 U 的第 i 行
for k in range(i, n):
# U[i, k] = A[i, k] - (L[i, 0] * U[0, k] + L[i, 1] * U[1, k] + ... + L[i, i-1] * U[i-1, k])
U[i, k] = A[i, k] - np.dot(L[i, :i], U[:i, k])
# 计算 L 的第 i 列
L[i, i] = 1 # 对角线元素设为 1，L 是单位下三角矩阵
for k in range(i+1, n):
# L[k, i] = (A[k, i] - (L[k, 0] * U[0, i] + L[k, 1] * U[1, i] + ... + L[k, i-1] * U[i-1, i])) / U[i, i]
L[k, i] = (A[k, i] - np.dot(L[k, :i], U[:i, i])) / U[i, i]
return L, U
def forward_substitution(L, f):
"""
使用前向替换法解线性方程组 L * y = f，其中 L 是下三角矩阵。
参数：
L (np.ndarray): 下三角矩阵，通常来自 LU 分解。
f (np.ndarray): 力向量（右侧向量）。
返回：
y (np.ndarray): 中间变量向量，使得 L * y = f。
"""
n = L.shape[0] # 获取矩阵的维度
y = np.zeros_like(f, dtype=np.float64) # 初始化 y 为与 f 同型的零向量
for i in range(n):
# y[i] = f[i] - (L[i, 0] * y[0] + L[i, 1] * y[1] + ... + L[i, i-1] * y[i-1])
y[i] = f[i] - np.dot(L[i, :i], y[:i])
# 由于 L 是单位下三角矩阵，对角线元素为 1，因此无需除以 L[i, i]
return y
def backward_substitution(U, y):
"""
使用后向替换法解线性方程组 U * x = y，其中 U 是上三角矩阵。
参数：
U (np.ndarray): 上三角矩阵，通常来自 LU 分解。
y (np.ndarray): 中间变量向量，使得 U * x = y。
返回：
x (np.ndarray): 解向量，使得 U * x = y。
"""
n = U.shape[0] # 获取矩阵的维度
x = np.zeros_like(y, dtype=np.float64) # 初始化 x 为与 y 同型的零向量
for i in reversed(range(n)):
# x[i] = (y[i] - (U[i, i+1] * x[i+1] + U[i, i+2] * x[i+2] + ... + U[i, n-1] * x[n-1])) / U[i, i]
x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
return x
def solve_displacements(A, f):
"""
通过 LU 分解及前向、后向替换法求解线性方程组 A * x = f。
参数：
A (np.ndarray): 系数矩阵。
f (np.ndarray): 力向量（右侧向量）。
返回：
L (np.ndarray): LU 分解中的下三角矩阵。
U (np.ndarray): LU 分解中的上三角矩阵。
y (np.ndarray): 前向替换得到的中间变量向量。
x (np.ndarray): 解向量，使得 A * x = f。
"""
L, U = lu_decomposition(A) # 进行 LU 分解，得到 L 和 U
y = forward_substitution(L, f) # 通过前向替换求解 L * y = f
x = backward_substitution(U, y) # 通过后向替换求解 U * x = y
return L, U, y, x
if __name__ == "__main__":
# 示例矩阵 A，为 3x3 对称正定矩阵
A = np.array([
[2, 1, 1],
[1, 2, 1],
[1, 1, 2]
], dtype=np.float64)
# 力向量 f
f = np.array([2, 2, 2], dtype=np.float64)
# 求解线性方程组 A * x = f
L, U, y, x = solve_displacements(A, f)
# 输出 LU 分解结果
print("LU分解结果：")
print("L =\n", L)
print("U =\n", U)
# 输出前向替换得到的中间变量 y
print("\n通过前向替换得到 y =\n", y)
# 输出后向替换得到的解向量 x
print("\n通过后向替换得到 位移向量 x =\n", x)
# 验证 LU 分解的正确性，通过 L 和 U 重构原矩阵 A
A_reconstructed = np.dot(L, U)
print("\n重构的 A =\n", A_reconstructed)
# 验证解的正确性，通过 A * x 计算得到的 f_computed 应与原始 f 相同
f_computed = np.dot(A, x)
print("\n验证 Ax =\n", f_computed)
print("原始 f =\n", f)
LDL\(^\text{T}\)分解
代码：
import numpy as np
def ldl_decomposition(A):
"""
对称矩阵 A 的 LDLᵗ 分解。
返回 L 和 D，其中 A = L D Lᵗ
"""
n = A.shape[0]
L = np.identity(n, dtype=np.float64)
D = np.zeros(n, dtype=np.float64)
for i in range(n):
# 计算 D[i]
D[i] = A[i, i] - np.dot(L[i, :i]**2, D[:i])
for j in range(i+1, n):
# 计算 L[j, i]
L[j, i] = (A[j, i] - np.dot(L[j, :i] * L[i, :i], D[:i])) / D[i]
return L, D
def forward_substitution(L, b):
"""
解下三角方程 L z = b
L 为单位下三角矩阵
"""
n = L.shape[0]
z = np.zeros_like(b, dtype=np.float64)
for i in range(n):
z[i] = b[i] - np.dot(L[i, :i], z[:i])
# 因为 L 是单位下三角矩阵，L[i, i] = 1
return z
def backward_substitution(Lt, y):
"""
解上三角方程 Lt x = y
Lt 为 L 的转置，即上三角矩阵
"""
n = Lt.shape[0]
x = np.zeros_like(y, dtype=np.float64)
for i in reversed(range(n)):
x[i] = (y[i] - np.dot(Lt[i, i+1:], x[i+1:])) / Lt[i, i]
return x
def solve_ldlt(A, f):
"""
使用 LDLᵗ 分解求解 A x = f
"""
# LDLᵗ 分解
L, D = ldl_decomposition(A)
# 前向替换，解 L z = f
z = forward_substitution(L, f)
# 对角线缩放，解 D y = z
y = z / D # 因为 D 是对角矩阵
# 后向替换，解 Lᵗ x = y
Lt = L.T
x = backward_substitution(Lt, y)
return L, D, z, y, x
if __name__ == "__main__":
# 示例矩阵
A = np.array([
[2, 1, 1],
[1, 2, 1],
[1, 1, 2]
], dtype=np.float64)
# 力向量
f = np.array([2, 2, 2], dtype=np.float64)
# 求解
L, D, z, y, x = solve_ldlt(A, f)
# 输出结果
np.set_printoptions(precision=8, suppress=True)
print("LDLᵗ分解结果：")
print("L =\n", L)
print("D =\n", D)
print("\n通过前向替换得到 z =\n", z)
print("通过对角线缩放得到 y =\n", y)
print("\n通过后向替换得到 位移向量 x =\n", x)
# 验证
A_reconstructed = L @ np.diag(D) @ L.T
print("\n重构的 A =\n", A_reconstructed)
# 验证 A x = f
f_computed = A @ x
print("\n验证 Ax =\n", f_computed)
print("原始 f =\n", f)
后处理
结构分析的本质，是通过分析结构在荷载作用下的应力分布，评估其性能与安全性。
结构分析的第一步，是确定施加在结构上的外部作用力，即荷载。
结构在荷载作用下，会发生变形，并在内部产生应力，它们是结构的响应。
我们用位移场描述结构中每个点在荷载作用下的位移，并用应变场描述结构中每个点由于变形产生的相对形变。
当我们计算出节点位移后，我们需要按如下步骤进行后处理：
-
从整体节点位移\(\{\Delta\}\)，提取特定单元e在全局坐标系下的节点位移\(\{\Delta\}^e\)；
-
利用\(\{\delta\}^e=[T]^e \{\Delta\}^e\)，把全局坐标系下的节点位移转换到局部坐标系下；
-
对于一维杆单元，计算单元内的应变：
\[\begin{equation}\varepsilon^e = \frac{u_2 - u_1}{L}\end{equation} \] -
计算单元内的应力：
\[\begin{equation}\sigma^e = E \varepsilon^e\end{equation} \] -
求约束反力：
\[\begin{equation} \{R\}=\begin{bmatrix}K\end{bmatrix}\{\Delta\}-\{F\} \end{equation} \]
标签：X.2,销接,刚度,矩阵,print,np,节点,位移,结构 From： https://www.cnblogs.com/yukibvvd/p/18599069/x2-connect-structure-2-zf1lvx