桁架有限元分析实例求解答

Question

（1） 描述什么是有限单元方法。
（2） 有限单元方法与其他数值方法相比，有哪些优势。
（3） 推导平面三节点三角形单元的形状函数
（4） 举例说明其在边界上的分布特点
（5） 请说明，为什么说平面三节点三角形单元为常应力单元!

（1） 对结构进行单元划分，列明节点编号和单元编号。
（2） 编写程序，计算各节点变形、应变和应力。
（3） 利用有限元商业软件对结构进行受力分析。
（4） 对不同方法所得到的结果进行对比分析，说明误差产生的主要原因。

Answer 1

引自免费微信小程序：皆我百晓生

（1）有限单元方法 有限单元法（Finite Element Method, FEM）是一种用于解决各种工程和物理问题的数值计算方法。这种方法将复杂的连续体（如结构、流体等）划分为许多互不重叠的子区域（称为有限单元），然后在每个单元内近似求解微分方程，再通过接口条件将所有单元连接起来，最终得到整个系统的近似解。这种方法尤其适用于处理具有复杂几何形状和边界条件的问题。

（2）有限单元方法的优势

• 适应性强：可以处理任意形状的域和复杂的边界条件。
• 局部性：仅需关注每个单元内的行为，无需全局信息。
• 灵活性：可方便地引入材料非线性、几何非线性和动力学效应。
• 并行计算：适合大规模并行计算，提高计算效率。
• 结果可视化：可以通过图形界面展示位移、应力等物理量的分布情况。

（3）平面三节点三角形单元的形状函数 对于一个平面三节点三角形单元，其三个节点分别标记为(A), (B), 和(C)。形状函数(N_i)定义了单元内任一点的位置与各个节点坐标的关联。对于这个三角形，我们可以定义三个形状函数如下：

[ N_A = \frac{(x - x_C)(y - y_B)}{(x_A - x_C)(y_B - y_C) + (x_B - x_C)(y_A - y_C)} ] [ N_B = \frac{(x - x_A)(y - y_C)}{(x_B - x_A)(y_C - y_A) + (x_C - x_A)(y_B - y_A)} ] [ N_C = \frac{(x - x_B)(y - y_A)}{(x_C - x_B)(y_A - y_C) + (x_A - x_B)(y_C - y_A)} ]

其中，(x)和(y)是任意点的坐标，而(x_i)和(y_i)是各个节点的坐标。

（4）形状函数在边界上的分布特点 在边界上，至少有一个节点位于边界上。例如，在边(AB)上，节点(C)位于内部，那么(N_C)在整个边(AB)上为零，只有(N_A)和(N_B)有非零值，且它们在端点处取值为1，中间部分线性变化。类似地，在边(BC)和(AC)上也有类似的特性。

（5）为什么平面三节点三角形单元为常应力单元 平面三节点三角形单元中的应力场在元素内部近似为常数。这意味着在单元内部，无论位置如何，应力都是相同的。这是因为在理想情况下，三角形单元内的位移场由线性插值确定，因此应变也是线性的，进而导致应力也是常数。然而，实际应用中，由于位移场的近似性质，这种假设并不完全准确，但在很多情况下足够精确。

接下来，我们将针对给定的桁架结构系统进行分析：

（6）单元划分和节点编号 我们可以通过将横梁分为两个单元，纵梁也分为两个单元，形成四个三角形单元。节点编号可以按照从左到右、从上到下的顺序依次为(1, 2, 3, 4)，单元编号可以表示为(U1, U2, U3, U4)。

（7）编写程序计算各节点变形、应变和应力 这里提供一个Python代码示例，使用基本的矩阵操作实现：

import numpy as np

# 杆件参数
L = 1e3  # 杆长
A = 50e-6  # 截面面积
E = 205e9  # 杨氏模量
nu = 0.3  # 泊松比
Ay = 235e6  # 屈服强度

# 单元和节点参数
nodes = np.array([[0, 0], [L, 0], [L, L], [0, L]])  # 节点坐标
elements = [(1, 2, 3), (2, 3, 4)]  # 元素列表，每个元素是一个包含三个节点编号的元组
loads = np.array([[-80e3, -90e3], [-80e3, -90e3]])  # 应力载荷

# 形状函数和雅可比矩阵
def shape_functions(x, y):
    return np.array([
        ((x - nodes[2][0]) * (y - nodes[1][1])) / (
            (nodes[0][0] - nodes[2][0]) * (nodes[1][1] - nodes[2][1]) +
            (nodes[1][0] - nodes[2][0]) * (nodes[0][1] - nodes[2][1])
        ),
        ((x - nodes[0][0]) * (y - nodes[2][1])) / (
            (nodes[1][0] - nodes[0][0]) * (nodes[2][1] - nodes[1][1]) +
            (nodes[2][0] - nodes[0][0]) * (nodes[1][1] - nodes[0][1])
        ),
        ((x - nodes[1][0]) * (y - nodes[0][1])) / (
            (nodes[2][0] - nodes[1][0]) * (nodes[0][1] - nodes[2][1]) +
            (nodes[0][0] - nodes[1][0]) * (nodes[2][1] - nodes[0][1])
        )
    ])

def jacobian(nodes, x, y):
    return np.array([
        [shape_functions(x, y)[i, 0], shape_functions(x, y)[i, 1]] for i in range(3)
    ]).T

# 基于刚度矩阵和质量矩阵的计算
def stiffness_matrix(nodes, elements, E, nu):
    K = np.zeros((len(nodes), len(nodes)))
    for e in elements:
        Ke = E * A / (2 * (1 - nu**2)) * jacobian(nodes, nodes[e[1]][0], nodes[e[1]][1]).T @ jacobian(nodes, nodes[e[1]][0], nodes[e[1]][1])
        for i in range(3):
            for j in range(3):
                K[e[i], e[j]] += Ke[i, j]
    return K

def mass_matrix(nodes, elements, rho):
    M = np.zeros((len(nodes), len(nodes)))
    for e in elements:
        Me = rho * A / 6 * np.array([
            [1, 1, 1],
            [1, 1, 1],
            [1, 1, 1]
        ])
        for i in range(3):
            for j in range(3):
                M[e[i], e[j]] += Me[i, j]
    return M

# 解析求解
K = stiffness_matrix(nodes, elements, E, nu)
M = mass_matrix(nodes, elements, 7850)  # 钢材密度
F = np.sum(loads, axis=0)  # 总载荷向量
u = np.linalg.solve(K, M @ F)  # 解得位移向量

# 计算应变和应力
eps = np.zeros(len(nodes))
sigma = np.zeros(len(nodes))
for e in elements:
    for i in range(3):
        eps[e[i]] = u[e[i+1], 0] - u[e[i], 0] / L
        sigma[e[i]] = E * (1 - nu) * eps[e[i]]

print("节点位移:", u)
print("节点应变:", eps)
print("节点应力:", sigma)

这段代码首先定义了形状函数和雅可比矩阵，然后根据这些函数计算了刚度矩阵和质量矩阵。接着，通过解析求解得到了节点位移，并进一步计算出应变和应力。

（8）利用有限元商业软件进行受力分析 这一步骤通常涉及具体软件的使用，例如ABAQUS、ANSYS或SolidWorks Simulation等。你需要导入模型，设置边界条件和载荷，然后运行分析。这部分需要具体的软件知识和操作技能，这里不再详述。

（9）对比分析及误差原因 对比两种方法（解析求解和商业软件分析）的结果，可能会发现一些差异。主要误差来源包括：

• 数值方法的近似误差：有限元方法基于离散化，会引入一定的近似误差。
• 分割网格的影响：单元大小和形状会影响结果精度。
• 商业软件可能考虑了更复杂的物理现象，比如非线性、接触等，而解析解往往简化了这些问题。
• 输入数据的精度：在商业软件中输入的模型参数可能存在小误差。

为了减小误差，可以尝试调整单元大小、优化网格划分，以及确保输入数据的准确性。此外，当解析