一维弹塑性分析程序：理论分析与实现

  帖子给出了一维弹塑性分析程序，基于弹塑性力学理论和有限元方法，采用隐式算法（Newton-Raphson 迭代法）求解非线性方程组。

1. 弹塑性力学基本理论

1.1 应力 - 应变关系

• 弹性阶段：应力与应变满足胡克定律      其中，       是应力，       是弹性模量，       是弹性应变。
• 弹塑性阶段：总应变分解为弹性应变和塑性应变      其中，       是总应变，       是塑性应变。

1.2 屈服准则

  采用 von Mises 屈服准则（适用于金属材料）：    其中，     是屈服应力。当  时，材料进入塑性状态。

1.3 硬化模型

  采用线性硬化模型，屈服应力随塑性应变增加：    其中，     是初始屈服应力，     是硬化模量。

2. 隐式算法（Newton-Raphson 迭代法）

2.1 非线性方程组求解

  弹塑性问题的平衡方程是非线性的，需迭代求解：    其中，     是内力向量，     是外力向量，     是位移向量。

2.2 迭代过程

1. 假设当前位移       ，计算残差力：      
2. 计算切线刚度矩阵       。
3. 求解位移增量方程：      
4. 更新位移：      
5. 重复上述步骤直到残差        小于容差。

3. 有限元离散化

3.1 单元划分

  将一维杆划分为      个单元，每个单元有 2 个节点。单元内的位移采用线性插值：    其中，     和      是形函数，     和      是节点位移。

3.2 应变计算

  单元应变与节点位移的关系：    其中，     是应变矩阵，     是单元节点位移向量。

4. 应力更新算法

4.1 弹性预测 - 塑性修正

1. 弹性预测：假设当前增量步为弹性，计算弹性应力试值：      
2. 屈服检查：计算屈服函数       。
3. 塑性修正：若 ，通过塑性乘子        修正应力：                  

4.2 切线模量

  塑性状态下的切线模量：    用于构建刚度矩阵，保证迭代收敛。

5. 边界条件处理

5.1 位移边界条件

  程序中采用固定 - 自由边界：

• 左端节点（节点 1）固定：      
• 右端节点（节点      ）施加位移：      通过修改刚度矩阵和力向量实现：

1. 将对应行和列清零，对角线元素设为 1。
2. 力向量对应位置设为指定位移值。

6. 收敛性控制

  采用双重收敛判据：

1. 残差力范数：      
2. 位移增量范数：      若超过最大迭代次数仍未收敛，则认为计算失败。

7.程序

  下面是matlab版本的程序。

%% 一维弹塑性分析程序 - 隐式算法求解
clear all; close all; clc;

%% 材料参数
E = 200e3;        % 弹性模量 (MPa)
sigma_y0 = 250;   % 初始屈服应力 (MPa)
H = 10e3;         % 硬化模量 (MPa)
nu = 0.3;         % 泊松比

%% 几何和网格参数
L = 1.0;          % 杆长度 (m)
num_elements = 10; % 单元数量
num_nodes = num_elements + 1; % 节点数量
dx = L / num_elements;      % 单元长度

%% 载荷和求解控制参数
max_disp = 0.02 * L;  % 最大位移 (m)
num_steps = 50;       % 载荷步数
tolerance = 1e-6;     % 收敛容差
max_iter = 20;        % 最大迭代次数

%% 初始化
% 节点位移、速度和加速度
u = zeros(num_nodes, 1);
v = zeros(num_nodes, 1);
a = zeros(num_nodes, 1);

% 单元应力、应变和塑性应变
sigma = zeros(num_elements, 1);
epsilon = zeros(num_elements, 1);
epsilon_p = zeros(num_elements, 1);

% 存储结果
results.u = zeros(num_nodes, num_steps+1);
results.sigma = zeros(num_elements, num_steps+1);
results.epsilon = zeros(num_elements, num_steps+1);
results.epsilon_p = zeros(num_elements, num_steps+1);
results.time = zeros(1, num_steps+1);

% 记录初始状态
results.u(:, 1) = u;
results.sigma(:, 1) = sigma;
results.epsilon(:, 1) = epsilon;
results.epsilon_p(:, 1) = epsilon_p;
results.time(1) = 0;

%% 主循环 - 隐式增量求解
for step = 1:num_steps
    % 当前载荷步的目标位移
    target_disp = max_disp * step / num_steps;
    
    % 位移边界条件 - 右端施加位移，左端固定
    bc_nodes = [1, num_nodes];  % 边界节点
    bc_values = [0, target_disp]; % 边界值
    
    % 牛顿迭代求解
    converged = false;
    iter = 0;
    u_prev = u;
    
    while ~converged && iter < max_iter
        iter = iter + 1;
        
        % 计算单元应变
        for e = 1:num_elements
            % 基于节点位移计算单元应变 (线性插值)
            epsilon(e) = (u(e+1) - u(e)) / dx;
        end
        
        % 计算残差力向量和切线刚度矩阵
        F_int = zeros(num_nodes, 1);  % 内力向量
        K_tangent = zeros(num_nodes, num_nodes); % 切线刚度矩阵
        
        for e = 1:num_elements
            % 高斯积分点 (单点积分)
            xi = 0;  % 自然坐标
            N = [0.5*(1-xi), 0.5*(1+xi)];  % 形函数
            dN_dx = [-1/dx, 1/dx];  % 形函数导数
            
            % 更新应力状态 (返回应力和切线模量)
            [sigma(e), epsilon_p(e), D_tangent] = update_stress(epsilon(e), epsilon_p(e), sigma(e), E, sigma_y0, H);
            
            % 单元刚度矩阵和内力向量
            K_e = D_tangent * dx * [1, -1; -1, 1];  % 单元刚度
            F_e = sigma(e) * dx * [1; -1];  % 单元内力
            
            % 组装全局刚度矩阵和内力向量
            indices = [e, e+1];
            F_int(indices) = F_int(indices) + F_e;
            K_tangent(indices, indices) = K_tangent(indices, indices) + K_e;
        end
        
        % 应用位移边界条件
        [K_tangent, F_int] = apply_bc(K_tangent, F_int, bc_nodes, bc_values);
        
        % 计算位移增量
        du = -K_tangent \ F_int;
        
        % 更新位移
        u = u_prev + du;
        
        % 检查收敛性
        residual_norm = norm(F_int);
        disp_norm = norm(du);
        
        if residual_norm < tolerance && disp_norm < tolerance
            converged = true;
            fprintf('载荷步 %d, 迭代 %d: 收敛成功 (残差 = %e, 位移增量 = %e)\n', step, iter, residual_norm, disp_norm);
        else
            fprintf('载荷步 %d, 迭代 %d: 残差 = %e, 位移增量 = %e\n', step, iter, residual_norm, disp_norm);
        end
        
        % 更新上一步位移
        u_prev = u;
    end
    
    % 记录当前载荷步结果
    results.u(:, step+1) = u;
    results.sigma(:, step+1) = sigma;
    results.epsilon(:, step+1) = epsilon;
    results.epsilon_p(:, step+1) = epsilon_p;
    results.time(step+1) = step / num_steps;
    
    % 检查是否收敛
    if ~converged
        warning('载荷步 %d 未收敛！', step);
    end
end

%% 可视化结果
visualize_results(results, E, sigma_y0, H);

%% 应力更新函数
function [sigma_new, epsilon_p_new, D_tangent] = update_stress(epsilon, epsilon_p_old, sigma_old, E, sigma_y0, H)
    % 弹性应变 = 总应变 - 塑性应变
    epsilon_e = epsilon - epsilon_p_old;
    
    % 计算弹性应力预测
    sigma_trial = sigma_old + E * epsilon_e;
    
    % 检查是否屈服
    sigma_y = sigma_y0 + H * epsilon_p_old;  % 考虑硬化的屈服应力
    f = abs(sigma_trial) - sigma_y;          % 屈服函数
    
    if f <= 0
        % 弹性状态
        sigma_new = sigma_trial;
        epsilon_p_new = epsilon_p_old;
        D_tangent = E;  % 弹性切线模量
    else
        % 塑性状态 - 需要迭代求解
        delta_gamma = f / (E + H);  % 塑性乘子
        sigma_new = sigma_trial - sign(sigma_trial) * E * delta_gamma;
        epsilon_p_new = epsilon_p_old + sign(sigma_trial) * delta_gamma;
        D_tangent = E * H / (E + H);  % 塑性切线模量
    end
end

%% 施加边界条件函数
function [K_modified, F_modified] = apply_bc(K, F, bc_nodes, bc_values)
    % 修改刚度矩阵和力向量以应用位移边界条件
    K_modified = K;
    F_modified = F;
    
    for i = 1:length(bc_nodes)
        node = bc_nodes(i);
        value = bc_values(i);
        
        % 修改刚度矩阵
        K_modified(node, :) = 0;
        K_modified(:, node) = 0;
        K_modified(node, node) = 1;
        
        % 修改力向量
        F_modified(node) = value;
    end
end

%% 结果可视化函数
function visualize_results(results, E, sigma_y0, H)
    % 提取结果
    num_steps = length(results.time) - 1;
    num_elements = size(results.sigma, 1);
    
    % 计算平均应力-应变曲线
    avg_strain = zeros(1, num_steps+1);
    avg_stress = zeros(1, num_steps+1);
    avg_epsilon_p = zeros(1, num_steps+1);
    
    for i = 1:num_steps+1
        avg_strain(i) = mean(results.epsilon(:, i));
        avg_stress(i) = mean(results.sigma(:, i));
        avg_epsilon_p(i) = mean(results.epsilon_p(:, i));
    end
    
    % 创建图形
    figure('Position', [100, 100, 1000, 800]);
    
    % 绘制应力-应变曲线
    subplot(2, 2, 1);
    plot(avg_strain, avg_stress, 'b-', 'LineWidth', 2);
    hold on;
    % 绘制弹性极限
    elastic_limit = sigma_y0 / E;
    plot([0, elastic_limit], [0, sigma_y0], 'r--', 'LineWidth', 1.5);
    grid on;
    xlabel('应变');
    ylabel('应力 (MPa)');
    title('应力-应变曲线');
    
    % 绘制塑性应变发展
    subplot(2, 2, 2);
    plot(results.time, avg_epsilon_p, 'r-', 'LineWidth', 2);
    grid on;
    xlabel('加载步');
    ylabel('塑性应变');
    title('塑性应变发展');
    
    % 绘制最终应力分布
    subplot(2, 2, 3);
    x_pos = linspace(0, 1, num_elements);
    plot(x_pos, results.sigma(:, end), 'g-', 'LineWidth', 2);
    grid on;
    xlabel('位置');
    ylabel('应力 (MPa)');
    title('最终应力分布');
    
    % 绘制应力-塑性应变曲线
    subplot(2, 2, 4);
    plot(avg_epsilon_p, avg_stress, 'm-', 'LineWidth', 2);
    grid on;
    xlabel('塑性应变');
    ylabel('应力 (MPa)');
    title('应力-塑性应变关系');
    
    % 显示材料参数
    sgtitle(sprintf('材料参数: E = %.1f GPa, σ_y0 = %.1f MPa, H = %.1f GPa', E/1e3, sigma_y0, H/1e3));
end    

总结

  该程序通过有限元方法将一维杆离散为多个单元，采用隐式 Newton-Raphson 迭代法求解非线性方程组。在每个载荷步内，使用弹性预测 - 塑性修正算法更新应力状态，并考虑材料硬化效应。通过合理处理边界条件和收敛性控制，确保计算结果的准确性和稳定性。

来源：有限元先生