几何非线性 FEM 程序｜一维杆 TL格式

本文是对 TL 格式弱形式的补充，对于单自由度一维杆进行程序实现。

• 非线性迭代算法采用牛顿-拉夫逊增量迭代算法
• 仿照 Abaqus 的习惯，对于分析步的计算，分别采用固定步长和自动步长的增量迭代计算

数值案例

考虑初始长度    、初始截面积    的两节点一维杆。节点 1 固定，节点 2 施加轴向集中力    ，尝试求解节点 2 的位移    。

固定增量步

将外部载荷分为相等的等份进行逐级加载

import math

import numpy as np

YOUNG_MODULUS = 200.0
INITIAL_LENGTH = 1.0
INITIAL_AREA = 1.0
TARGET_FORCE = 100.0


def deformation_gradient(u):
    """Return the deformation gradient F of the uniaxial bar."""
    displacement_gradient = u / INITIAL_LENGTH

    F = np.array([
        [1.0 + displacement_gradient]
    ])

    return F


def green_lagrange_strain(F):
    """Return the Green-Lagrange strain E = 0.5 * (F.T F - I)."""
    identity = np.eye(1)

    E = 0.5 * (F.T @ F - identity)

    return E


def second_p_k_stress(E):
    """Return the second Piola-Kirchhoff stress S = Young's modulus * E."""
    S = YOUNG_MODULUS * E

    return S


def internal_force(F, S):
    """Return the internal force obtained from the TL internal virtual work."""
    F11 = F[0, 0]
    S11 = S[0, 0]

    force = S11 * F11 * INITIAL_AREA

    return force


def tangent_stiffness(F, S):
    """Return the consistent tangent stiffness Kt = k_mat + k_geo."""
    F11 = F[0, 0]
    S11 = S[0, 0]

    k_mat = (
        YOUNG_MODULUS
        * INITIAL_AREA
        / INITIAL_LENGTH
        * F11**2
    )
    k_geo = INITIAL_AREA / INITIAL_LENGTH * S11
    tangent = k_mat + k_geo

    return tangent


def calculate_state(u):
    """Calculate all physical quantities from the current displacement."""
    F = deformation_gradient(u)

    E = green_lagrange_strain(F)
    S = second_p_k_stress(E)
    force = internal_force(F, S)

    return F, E, S, force


def convergence_residual(residual, target_force):
    """Return the normalized squared residual used in the original example."""
    return residual**2 / (1.0 + target_force**2)


def newton_raphson(
    initial_u,
    target_force,
    tol=1.0e-5,
    max_iter=20,
):
    """Solve one load increment with the Newton-Raphson method."""
    trial_u = float(initial_u)
    iteration = 0

    F, E, S, force = calculate_state(trial_u)
    residual = target_force - force
    conv_residual = convergence_residual(residual, target_force)

    print(
        f"{'iter':>5} {'u2':>10} {'E11':>10} "
        f"{'S11':>12} {'conv_R':>12}"
    )
    print(
        f"{iteration:5d} {trial_u:10.5f} {E[0, 0]:10.5f} "
        f"{S[0, 0]:12.3e} {conv_residual:12.3e}"
    )

    while conv_residual > tol:
        if iteration >= max_iter:
            raise RuntimeError("Newton-Raphson iteration did not converge.")

        tangent = tangent_stiffness(F, S)
        if math.isclose(tangent, 0.0, abs_tol=1.0e-14):
            raise ZeroDivisionError("Zero tangent stiffness.")

        delta_u = residual / tangent
        trial_u += delta_u

        iteration += 1
        F, E, S, force = calculate_state(trial_u)
        residual = target_force - force
        conv_residual = convergence_residual(residual, target_force)

        print(
            f"{iteration:5d} {trial_u:10.5f} {E[0, 0]:10.5f} "
            f"{S[0, 0]:12.3e} {conv_residual:12.3e}"
        )

    return trial_u


def incremental_force_fixed(
    tol=1.0e-5,
    max_iter=20,
    num_increments=10
):
    """Apply the total force in equal increments and solve every increment."""

    force_history = [0.0]
    displacement_history = [0.0]
    accepted_u = 0.0

    print(
        f"\nInitial state: force = {0.0:.3f}, u2 = {accepted_u:.5f}, "
        f"E11 = {0.0:.5f}, S11 = {0.0:.3e}"
    )

    for j in range(1, num_increments + 1):
        current_force = j / num_increments * TARGET_FORCE
        print(
            f"\nIncrement {j}: current force = {current_force:.5f}"
        )

        accepted_u = newton_raphson(
            accepted_u,
            current_force,
            tol,
            max_iter,
        )

        force_history.append(current_force)
        displacement_history.append(accepted_u)

    return force_history, displacement_history, accepted_u


def main():
    incremental_force_fixed()


if __name__ == "__main__":
    main()

结果罗列如下：

自动增量步

根据迭代情况自动调整增量步。

参数：

• 最多增量步：max_num_increment
• 初始增量步大小：initial_increment
• 最小增量步大小：mini_increment
• 最大增量步大小：max_increment

调整规则：

1. 一个 increment在 16 次迭代以内收敛：进入下一个 increment
2. 连续两个 increment都在 5 次迭代以内收敛：下一个 increment增加 50%，但不得超过 max_increment
3. 当前 increment迭代 16 次仍不收敛：将当前 increment折减至原来的 25%，重新计算
4. 出现以下情况时停止计算并报错：
• 折减次数超过 5 次
• 折减后的 increment小于 mini_increment
• num_increment超过 max_num_increment

def incremental_force_automatic(
    tol=1.0e-5,
    max_num_increments=20,
    initial_increment=0.1,
    mini_increment=1.0e-4,
    max_increment=1.0,
    max_iter=16,
    max_cutbacks=5,
    load_tol=1.0e-12,
):
    """Apply TARGET_FORCE using automatic normalized load increments."""
    if mini_increment <= 0 or max_increment <= 0 or initial_increment <= 0:
        raise ValueError("Increment sizes must be positive.")
    if mini_increment > max_increment:
        raise ValueError("mini_increment cannot be larger than max_increment.")
    if initial_increment > max_increment:
        raise ValueError("initial_increment cannot be larger than max_increment.")
    if max_num_increments <= 0:
        raise ValueError("max_num_increments must be positive.")

    force_history = [0.0]
    displacement_history = [0.0]
    accepted_u = 0.0

    print(
        f"\nInitial state: force = {0.0:.3f}, u2 = {accepted_u:.5f}, "
        f"E11 = {0.0:.5f}, S11 = {0.0:.3e}"
    )

    current_load_factor = 0.0
    next_increment_size = initial_increment
    accepted_increment = 0
    fast_convergence_count = 0

    while current_load_factor < 1.0 - load_tol:
        if accepted_increment >= max_num_increments:
            raise RuntimeError(
                f"Reached max_num_increments={max_num_increments} "
                f"before TARGET_FORCE."
            )

        cutback_count = 0
        current_increment_size = min(
            next_increment_size,
            1.0 - current_load_factor,
        )

        while True:
            target_load_factor = (
                current_load_factor + current_increment_size
            )
            current_force = target_load_factor * TARGET_FORCE

            print(
                f"\nIncrement {accepted_increment + 1}: attempt {cutback_count + 1}, "
                f"increment size = {current_increment_size:.5f}, "
                f"load factor = {target_load_factor:.5f}, "
                f"current force = {current_force:.5f}"
            )

            trial_u, iteration, converged = newton_raphson(
                accepted_u,
                current_force,
                tol,
                max_iter,
            )

            if converged:
                accepted_u = trial_u
                current_load_factor = target_load_factor
                accepted_increment += 1

                print(
                    f"Converged force = {current_force:.5f}, "
                    f"iteration = {iteration}"
                )
                force_history.append(current_force)
                displacement_history.append(accepted_u)

                if cutback_count == 0 and iteration <= 5:
                    fast_convergence_count += 1
                else:
                    fast_convergence_count = 0

                if fast_convergence_count >= 2:
                    next_increment_size = min(
                        1.5 * current_increment_size,
                        max_increment,
                    )
                    fast_convergence_count = 0
                else:
                    next_increment_size = min(current_increment_size, max_increment)

                break

            cutback_count += 1
            fast_convergence_count = 0

            if cutback_count > max_cutbacks:
                raise RuntimeError(
                    f"Increment {accepted_increment + 1} failed after "
                    f"{max_cutbacks} cutbacks."
                )

            current_increment_size *= 0.25

            if current_increment_size < mini_increment:
                raise RuntimeError(
                    f"Increment size {current_increment_size:.3e} is smaller than "
                    f"mini_increment {mini_increment:.3e}."
                )

    return force_history, displacement_history, accepted_u