探索莫尔圆与强度包络线：一份Python实现指南

在岩土工程领域，莫尔圆 (Mohr Circles) 和 强度包络线 (Strength Envelope) 是分析土体应力状态的重要工具。如果你对土力学感兴趣，或者想深入了解应力分析的世界，这篇文章将为你提供一个极具吸引力的Python实现，帮助你自动生成莫尔圆并找到它们的公切线，也就是材料的强度包络线。

在这篇博客中，我们将介绍：

1. 莫尔圆与强度包络线的简要概述。
2. 寻找莫尔圆公切线的数学原理。
3. 一份自动生成莫尔圆并计算强度包络线的Python代码。
4. 结合代码和可视化的实际示例。

什么是莫尔圆和强度包络线？

莫尔圆是一种用于图形化表示应力状态的工具，尤其在岩土工程中，它帮助我们理解不同应力条件下土体的行为。每一个莫尔圆代表了某种应力下的状态，而这些圆的公切线（即强度包络线）代表了材料的极限状态，或者说材料的强度。

强度包络线可以用来确定材料在不同法向应力下的屈服强度，从而在实际工程中为设计提供参考。通过绘制莫尔圆并找到其公切线，我们可以计算材料的凝聚力 (cohesion) 和摩擦角 (friction angle)，这是土力学中两个至关重要的强度参数。

数学背景

对于三个莫尔圆来说，公切线的存在意味着必须满足一定的条件。我们通过以下步骤来判断是否存在公切线：

1. 对每两个圆之间的应力差进行计算。
2. 根据圆心坐标和半径的差值来判断是否能够满足包络线的几何条件。
3. 计算切线的斜率与截距，并验证是否满足第三个圆的条件。

该方法的复杂性在于需要遍历所有可能的符号组合来验证切线的可行性。

Python实现

在这份代码中，我们首先生成了一组可能的圆心与半径组合，并检查这些组合是否满足几何条件。如果满足条件，我们会进一步计算强度包络线的参数：摩擦角和凝聚力，并绘制出相应的莫尔圆与公切线。

核心代码部分

检查条件

def check_conditions(x_list, r_list):
    # 所有可能的符号组合
    signs = list(itertools.product([-1, 1], repeat=len(x_list)))
    for s in signs:
        s1, s2, s3 = s
        x1, x2, x3 = x_list
        r1, r2, r3 = r_list
        # 计算Δx和Δr
        dx12 = x1 - x2
        dx13 = x1 - x3
        dr12 = s1 * r1 - s2 * r2
        dr13 = s1 * r1 - s3 * r3
        # 计算分母，检查是否为正
        denominator1 = dx12 ** 2 - dr12 ** 2
        denominator2 = dx13 ** 2 - dr13 ** 2
        if denominator1 <= 0 or denominator2 <= 0:
            continue
        # 计算m^2
        m_squared1 = dr12 ** 2 / denominator1
        m_squared2 = dr13 ** 2 / denominator2
        # 检查两个m^2是否相等
        if not np.isclose(m_squared1, m_squared2, atol=1e-6):
            continue
        m_squared = m_squared1
        if m_squared <= 0:
            continue
        m = np.sqrt(m_squared)
        # 计算b
        sqrt_term = np.sqrt(m ** 2 + 1)
        b = s1 * r1 * sqrt_term - m * x1
        # 验证第三个方程
        left = m * x3 + b
        right = s3 * r3 * sqrt_term
        if not np.isclose(left, right, atol=1e-6):
            continue
        return True
    return False

寻找公切线

def find_common_tangent(x_list, r_list):
    solutions = []
    signs = list(itertools.product([-1, 1], repeat=len(x_list)))
    for s in signs:
        # 计算相应的斜率与截距
        m_squared1 = ...
        m_squared2 = ...
        if not np.isclose(m_squared1, m_squared2, atol=1e-6):
            continue
        m = np.sqrt(m_squared1)
        b = s1 * r1 * sqrt_term - m * x1
        left = m * x3 + b
        right = s3 * r3 * sqrt_term
        if not np.isclose(left, right, atol=1e-6):
            continue
        phi = np.degrees(np.arctan(m))
        cohesion = b / sqrt_term
        solutions.append((m, b, s, phi, cohesion))
    return solutions

生成组合并绘制图形

def plot_circles_and_tangent(x_list, r_list, solutions):
    fig, ax = plt.subplots(figsize=(10, 6))
    # 绘制莫尔圆
    for idx, (x_c, r) in enumerate(zip(x_list, r_list)):
        theta = np.linspace(0, np.pi, 200)
        x = x_c + r * np.cos(theta)
        y = r * np.sin(theta)
        ax.plot(x, y, color='blue')
        ax.plot(x_c, 0, 'bo')
        ax.text(x_c, -3, f'Center: ({x_c}, 0)\nRadius: {r}', ha='center')
    # 绘制公切线
    if solutions:
        for solution in solutions:
            m, b, s, phi, cohesion = solution
            x_vals = np.linspace(-10, max(x_list) + max(r_list) + 10, 400)
            y_vals = m * x_vals + b
            ax.plot(x_vals, y_vals, color='red', linestyle='--')
    else:
        ax.text(0.5, 0.9, '未找到公切线。', transform=ax.transAxes, fontsize=14, ha='center', color='red')
    plt.show()

运行示例

运行该程序后，你可以输入圆心坐标和半径的组合，程序将计算出相应的强度包络线，并绘制出莫尔圆图及其公切线。如下图所示：

# 自动生成满足条件的组合
valid_combinations = generate_valid_combinations()

if valid_combinations:
    print("以下是满足条件的圆心和半径组合：")
    for idx, (x_list, r_list) in enumerate(valid_combinations):
        print(f"组合 {idx + 1}:")
        print(f"  Centers (σ): {x_list}")
        print(f"  Radii (τ): {r_list}")
    # 用户选择组合
    choice = int(input(f"\n请输入要选择的组合编号 (1 - {len(valid_combinations)}): "))
    selected_x_list, selected_r_list = valid_combinations[choice - 1]
    solutions = find_common_tangent(selected_x_list, selected_r_list)
    plot_circles_and_tangent(selected_x_list, selected_r_list, solutions)
else:
    print("未找到满足条件的组合。")

总结

通过这份Python实现，我们可以快速生成莫尔圆，并自动寻找它们的公切线，从而计算出土体的强度参数。这为我们提供了一个强大且灵活的工具，用于分析材料的应力和强度特性。不论是在教学中，还是在实际的工程应用中，这样的程序都能带来极大的帮助。

如果你对岩土工程感兴趣，或者正在寻找编程与土力学的结合点，这份代码将是你不容错过的参考工具。希望这篇文章能为你带来启发，帮助你更好地理解莫尔圆与强度包络线的奥秘！

欢迎在评论区分享你的想法和改进建议！