基于MATLAB编程实现的平面三角形单元通用有限元分析程序

基于MATLAB编程实现的平面三角形单元通用有限元分析程序

这是本人学习有限元课程的一次大作业，由个人独立完成，开源以供学弟学妹以及感兴趣的同学参考学习，本文仅供学习使用，禁止任何形式的转载

原大作业题目

解决思路

使用MATLAB软件，按照有限元的一般编程方法编程，为了提高泛用性，程序从inp文件读取模型信息，这一逻辑与abaqus类似，此外，本人还编程实现了自动网格划分，可以针对给定的单元数量划分不同密度的网格。此外还实现了全过程的可视化，能够直观地展示计算结果。最后与有限元仿真结果进行了对比。

inp文件格式说明

inp示例文件如下。

* Config
dimension, 2

* Material
210e3, 0.2

* Nodes
5
1, 0.000000, -0.000000
2, 1000.000000, -0.000000
3, 0.000000, -1000.000000
4, 1000.000000, -1000.000000
5, 2000.000000, -1000.000000

* Elements
tri2d_stress, 3, 10.000000
1, 1, 3, 4
2, 1, 4, 2
3, 2, 4, 5

* Forces
1
1, 10000.000000, 0.000000

* Boundaries
3
5, all
4, all
3, all

*后跟定义的关键字，其后再跟数据的内容。
值得说明的是这里的规则都是自己定义的，不是死板的，因为博主平时用ABAQUS多，因此这里的规则设定也是参考了ABAQUS中的INP文件的思路，这种从文件读写类的程序通用性都比较好，而且我在代码中考虑了单元类型维度等等因素，尽管目前代码只能处理平面三角形单元的有限元问题，但是在代码的框架中考虑了未来的更新拓展，以后要在这个框架内实现其他单元也是十分方便的。

INP文件中，Config定义一些全局的常数，Material定义材料相关的参数，Nodes定义节点信息，Elements定义单元信息，Forces定义外力信息，Boundaries定义边界条件信息，all关键字表示全方向固定。

源代码

为了编写的高效以及后期维护方便，使用了多文件编程方法。同时建议工科学生（非计算机的）不要光闷着头编程，要及时学习一些高效的计算机编程技术与框架，比如面向对象编程，在我看来这是非常好的编程框架，但是我接触的大多数工科生还在用函数式编程，甚至函数式编程都不会的，但是却写出了及其复杂的程序，然后在项目交接的时候，维护他代码的人就会非常痛苦。 所以既然我们使用了编程工具，还是希望大家把他打磨到精湛的程度，能够写出不输专业人员水准的代码。

前置类Ref.m

定义传递使用句柄的类，用来实例化 整体刚度矩阵 和 整体载荷列阵，这样这两个变量在函数间可以按句柄传递，不用开辟新的内存空间，可以节约内存提高运行效率。这步是因为在MATLAB里貌似函数参数要按引用传递（非计算机可能不知道我在说什么，这里可以简单理解为省内存），这种用句柄的子类是比较方便的做法，其他语言比如C语言，实现这个可能就会用更方便的指针工具。

%% Copyright: YunHao Liu
% Define Ref
classdef Ref < handle
    properties
        value;
    end

    methods
        function obj = Ref(value)
            obj.value = value;
        end
    end
end

① main.m

主程序入口，定义一些使用到的参数，调用其他函数

%% Copyright: YunHao Liu
%% Variable Define
plot_color = false;
num_seeds = 1;
path = 'inp.txt';
%% Generate the Inp file
GenerateInpFromNumSeeds(num_seeds)

%% Read Input DATA
[material, nodes, elements, forces, boundaries, main_config] = ReadData(path);
% parameters from config
dimension = int64(str2double(main_config(1, 1)));
element_type = main_config(1, 2);
thickness = str2double(main_config(1, 3));

%% Initial Display
triang = triangulation(elements, nodes);  % mesh on 2d object
figure;
triplot(triang, '-xb');  % 2d tri plot
axis equal;
hold on

%% Get element_stiffness and Assemble it to Global K
switch element_type
    case "tri2d_stress"
        % Define Ref K
        K = Ref(zeros(length(nodes)*2, length(nodes)*2));
        % Define Ref F
        F = Ref(zeros(length(nodes)*2, 1));
        % Get D
        D = Ref(GetD("plane_stress", material));    
        % Calculate the K, F
        FigureOutKF(K, F, nodes, elements, D, forces, thickness)
        % Boundary condition process
        PenaltyBounary(K, F, boundaries)
        % Solve
        U = SolveKF(K, F);
        % Assign U
        nodes_end = AssignU(nodes, 3000 * U);
        % Calculate the stress
        sigma = FigureOutStress(U, D, nodes, elements);
        mises = ToMises(sigma);

end

%% Finial Display
% Plot patch according to mises
if plot_color
    % With color
    colormap('turbo')
    clim([min(mises), max(mises)])
    size_elements = size(elements);
    len_elements = size_elements(1);
    for e_dof = 1:len_elements
        % Get Dof of each Nodes
        n1 = elements(e_dof, 1);
        n2 = elements(e_dof, 2);
        n3 = elements(e_dof, 3);
        % Get the x, y on each element node
        e_nodes = zeros(3, 2);
        e_nodes(1, :) = nodes_end(n1, :);
        e_nodes(2, :) = nodes_end(n2, :);
        e_nodes(3, :) = nodes_end(n3, :);
        % Patch
        patch(e_nodes(:,1), e_nodes(:,2), ...
            [mises(e_dof, 1); mises(e_dof, 1); mises(e_dof, 1)], ...
            'FaceAlpha', 0.9, 'EdgeColor', 'none')
    end
else
    % No color
    triang = triangulation(elements, nodes_end);  % mesh on 2d object
    triplot(triang, '--or');  % 2d tri plot
end
axis equal;

%% Calculate Von Mises Stress
function mises = ToMises(sigma)
    size_sigma = size(sigma);
    len = size_sigma(1);
    mises = NaN * zeros(len, 1);
    for i = 1:len
        sig_m = (sigma(i, 1) + sigma(i, 2))/3;
        mises(i, 1) = sqrt(1.5*((sigma(i, 1)-sig_m)^2 + (sigma(i, 2)-sig_m)^2 + sig_m^2 + 2*sigma(i, 3)^2));
    end
end

② ReadData.m

从inp文件中读取信息，并转化为变量。

%% Copyright: YunHao Liu
% Read the information from the input file
function [Materials, Nodes, Elements, Forces, Boundaries, Configs] = ReadData(path)
    % Open the input file
    f = fopen(path);
    
    % Define some needed variables
    input_var = "";
    element_type = "";
    thickness = 0.0;
    num_nodes_per_element = 0;
    dimension = int8(2);
    in_data = false;
    Materials = zeros(2);
    num_nodes = 0;
    num_elements = 0;
    num_boudaries = 0;

    %% Read the line in the file
    line = fgetl(f);
    while ischar(line)
        %% Get the key variable name
        if ~isempty(line) && line(1)=='*'
            input_var = lower(line(2:end));
            input_var = strrep(input_var, ' ', '');
            in_data = false;
            line = fgetl(f);  % Next line and continue
            continue
        end

        %% Read the data of the variable
        if ~isempty(input_var) && ~isempty(line)
            %% Read the Config Parameters
            if input_var == "config"
                % no indata check, means no data, only need parameters in
                % this variable
                line = remove_split_convert(line, ' ', ',', false);
                % config parameter define as the key, value
                % for example: dimension, 2
                switch string(line(1, 1))
                    case 'dimension'
                        dimension = int8(str2double(string(line(1, 2))));
                end
            end

            %% Read Material Parameters
            if input_var == "material"  
                % no indata check, means no data, only need parameters in
                % this variable
                Materials = remove_split_convert(line, ' ', ',', true);
            end
            
            %% Read Nodes information
            if input_var =="nodes"
                % indata check
                if in_data == false  % if false, read pre parameters first
                    % Pre parameters define the size of the following data
                    % or some other needed information
                    num_nodes = remove_split_convert(line, ' ', ',', true);
                    num_nodes = num_nodes(1, 1);
                    Nodes = zeros(num_nodes, dimension);
                    in_data = true;
                else  % if true, read data
                    for i = 1:num_nodes
                        line = remove_split_convert(line, ' ', ',', true);
                        node_dof = int64(line(1, 1));
                        % each of dimension
                        for dim = 1:dimension
                            Nodes(node_dof, dim) = line(1, dim + 1);
                        end
                        % next line
                        if i ~= num_nodes
                            line = fgetl(f);
                        end
                    end
                end
            end

            %% Read Element information
            if input_var == "elements"
                % indata check
                if in_data == false
                    line = remove_split_convert(line, ' ', ',', false);
                    element_type = string(line(1,1));
                    num_elements = int64(str2double(line(1, 2)));
                    switch element_type
                        case "tri2d_stress"
                            num_nodes_per_element = 3;
                            thickness = str2double(line(1, 3));
                            Elements = zeros(num_elements, 3);
                    end
                    in_data = true;
                else 
                    for i = 1:num_elements
                        line = remove_split_convert(line, ' ', ',', true);
                        line = int64(line);
                        element_dof = line(1, 1);
                        for i_node = 1:num_nodes_per_element
                            Elements(element_dof, i_node) = ...
                                line(1, i_node + 1);
                        end
                        % next line
                        if i ~= num_elements
                            line = fgetl(f);
                        end
                    end
                end
            end

            %% Read Forces information
            if input_var == "forces"
                % indata check
                if in_data == false
                    line = remove_split_convert(line, ' ', ',', true);
                    num_forces = line(1, 1);
                    Forces = zeros(num_forces, 3);
                    in_data = true;
                else
                    for i = 1:num_forces
                        line = remove_split_convert(line, ' ', ',', true);
                        line_length = length(line);
                        % dof of force
                        Forces(i, 1) = int64(line(1, 1));
                        for i_param = 2:line_length
                            Forces(i, i_param) = line(1, i_param);
                        end
                    end
                end
            end

            %% Read Boundaries information
            if input_var == "boundaries"
                % indata check
                if in_data == false
                    line = remove_split_convert(line, ' ', ',', true);
                    num_boudaries = line(1, 1);
                    Boundaries = NaN * ones(num_boudaries, 3);
                    in_data = true;
                else
                    for i = 1:num_boudaries
                        line = remove_split_convert(line, ' ', ',', false);
                        line_length = length(line);
                        % dof of boundary
                        Boundaries(i, 1) = int64(str2double(line(1, 1)));
                        % for example:
                        % 1, x, y=10
                        for i_param = 2:line_length
                            cond_words = string(line(1, i_param));
                            switch cond_words
                                case 'x'
                                    Boundaries(i, 2) = 0;
                                case 'y'
                                    Boundaries(i, 3) = 0;
                                case 'all'
                                    Boundaries(i, 2) = 0;
                                    Boundaries(i, 3) = 0;
                                otherwise
                                    cond_words = strsplit(cond_words, '=');
                                    cond_key = string(cond_words(1, 1));
                                    cond_value = str2double(cond_words(1, 2));
                                    switch cond_key
                                        case 'x'
                                            Boundaries(i, 2) = cond_value;
                                        case 'y'
                                            Boundaries(i, 3) = cond_value;
                                    end
                            end
                        end
                        % next line
                        if i ~= num_boudaries
                            line = fgetl(f);
                        end
                    end
                end
            end
        end
        
        % Get new line
        line = fgetl(f);
    end
    
    fclose(f);  % Close file

    %% Output to Configs
    % Configs = [dimension, element type, thickness]
    Configs = [
        int2str(dimension),...
        element_type,...
        string(thickness)
    ];
    %% disp for Debug
    % disp(Configs)
    % disp(Materials)
    % disp(Nodes)
    % disp(Elements)
end

%% Function to handle a str
% remove the space, split, and convert data type if needed
function data_out = remove_split_convert(string_in, remove, split, convert)
    % Remove the space
    string_out = strrep(string_in, remove, '');

    % Split with char given
    string_out = strsplit(string_out, split);
    
    % Convert to double type
    if convert == true
    data_out = zeros(1, length(string_out));
        for i = 1:length(string_out)
            data_out(1, i) = str2double(string_out(1, i));
        end
    else
        data_out = string_out;
    end
end

③ FigureOutKF.m

求解出整体刚度矩阵和整体载荷列阵，该函数是一个封装函数，整合了ElementStiffness.m，Assemble2dGlobalStiffness.m和Assemble2dGlobalForce.m

%% Copyright: YunHao Liu
% FigureOut the K and F, both them are Ref class, global F and global K
function FigureOutKF(K, F, nodes, elements, D, forces, thickness)
    size_elements = size(elements);
    % Calculate the GlobalK
    for i_elem = 1:size_elements(1)
        e_nodes = zeros(3, 2);
        e_nodes(1, :) = nodes(elements(i_elem, 1), :);
        e_nodes(2, :) = nodes(elements(i_elem, 2), :);
        e_nodes(3, :) = nodes(elements(i_elem, 3), :);
        k_e = ElementStiffness("tri2d_stress", e_nodes, D, thickness);
        % Assemble to Global K
        Assemble2dGlobalStiffness(k_e, K, elements(i_elem, :))
    end
    % Assemble Global F
    Assemble2dGlobalForce(forces, F)
end

④ ElementStiffness.m

求解单元刚度矩阵

%% Copyright: YunHao Liu
% Calculate the element stiffness K (K_a actually)
function K = ElementStiffness(element_type, nodes, D, param)
    switch element_type
        case "tri2d_stress"
            % Get parameters
            thickness = param(1, 1);
            % Define Nodes
            x1 = nodes(1, 1);
            y1 = nodes(1, 2);
            x2 = nodes(2, 1);
            y2 = nodes(2, 2);
            x3 = nodes(3, 1);
            y3 = nodes(3, 2);
            % Define Jacobi
            Jacobi = [
                x1 - x3, y1 - y3;
                x2 - x3, y2 - y3;
            ];
            det_Jacobi = det(Jacobi);
            % Define B
            B = [
                y2-y3, 0, y3-y1, 0, y1-y2, 0;
                0, x3-x2, 0, x1-x3, 0, x2-x1;
                x3-x2, y2-y3, x1-x3, y3-y1, x2-x1, y1-y2;
            ]/det_Jacobi;
            % Calculate K
            K = 0.5 * thickness * abs(det_Jacobi) * B' * D.value * B;
    end
end

⑤ Assemble2dGlobalStiffness.m

将单刚组装到总刚

%% Copyright: YunHao Liu
% Assemble 2d global stiffness, k_e is the element stiffness, K is a Ref
% class, K is the global stiffness, location is the node dof of the element
function Assemble2dGlobalStiffness(k_e, K, location)
    num_nodes = length(location);
    for i = 1:num_nodes
        for j = 1:num_nodes
            dim1 = location(i);
            dim2 = location(j);
            K.value(2*dim1-1:2*dim1, 2*dim2-1:2*dim2) = ...
                K.value(2*dim1-1:2*dim1, 2*dim2-1:2*dim2) ...
                + k_e(2*i-1:2*i, 2*j-1:2*j);
        end
    end
end

⑥ Assemble2dGlobalForce.m

将节点载荷组装到总体载荷列阵

%% Copyright: YunHao Liu
% Assemble global force, F is a Ref class, f_e is the node dof of the
% element
function Assemble2dGlobalForce(f_e, F)
    dof = f_e(1, 1);
    F.value(2*dof-1) = F.value(2*dof-1) + f_e(1, 2);
    F.value(2*dof) = F.value(2*dof) + f_e(1, 3);
end

⑦ PenaltyBounary.m

利用罚函数法处理边界条件

%% Copyright: YunHao Liu
% Penalty method to process the Bounary conditions
function PenaltyBounary(K, F, boundaries)
    size_boundaries = size(boundaries);
    num_boundaries = size_boundaries(1);
    penalty = 1e5 * max(K.value, [], 'all');
    for i = 1:num_boundaries
        dof = boundaries(i, 1);
        x_value = boundaries(i, 2);
        y_value = boundaries(i, 3);
        if ~isnan(x_value)
            K.value(2*dof-1, 2*dof-1) = ...
                K.value(2*dof-1, 1) + penalty;
            F.value(2*dof-1, 1) = ...
                F.value(2*dof-1, 1) + penalty * x_value;
        end
        if ~isnan(y_value)
            K.value(2*dof, 2*dof) = ...
                K.value(2*dof, 2*dof) + penalty;
            F.value(2*dof, 1) = ...
                F.value(2*dof, 1) + penalty * y_value;
        end
    end
end

⑧ SolveKF.m

求解结构的位移

%% Copyright: YunHao Liu
function U = SolveKF(K, F)
    U = K.value\F.value;
end

⑨ GetD.m

根据材料常数 求解问题的D矩阵

%% Copyright: YunHao Liu
% Get D matrix from the constant value and specific type
function D = GetD(type, material)
    E = material(1, 1);
    Mu = material(1, 2);
    switch type
        case "plane_stress"
            D = (E/(1.0-Mu^2))*[
                1.0,  Mu, 0.0;
                Mu , 1.0, 0.0;
                0.0, 0.0, 0.5*(1.0-Mu);
            ];
    end
end

⑩ AssignU.m

将求解出的位移附加到原结构上（用于可视化变形后的结构）

%% Copyright: YunHao Liu
% Assign the U (displacement) to the nodes information, the output
% nodes_end assigned U is the nodes after deformation.
function nodes_end = AssignU(nodes, U)
    nodes_end = nodes;
    size_nodes = size(nodes);
    for dof = 1:size_nodes(1)
        % Assign x
        nodes_end(dof, 1) = nodes_end(dof, 1) + U(2*dof-1);
        % Assign y
        nodes_end(dof, 2) = nodes_end(dof, 2) + U(2*dof);
    end
end

⑪ FigureOutStress.m

求解整个结构的应变信息

function sigma = FigureOutStress(Q, D, nodes, elements)
    size_elements = size(elements);
    n_elements = size_elements(1);
    sigma = NaN * zeros(n_elements, 3);
    for e_dof = 1:n_elements
        % Get Dof of each Nodes
        n1 = elements(e_dof, 1);
        n2 = elements(e_dof, 2);
        n3 = elements(e_dof, 3);
        % Get the x, y on each element node
        e_nodes = zeros(3, 2);
        e_nodes(1, :) = nodes(n1, :);
        e_nodes(2, :) = nodes(n2, :);
        e_nodes(3, :) = nodes(n3, :);
        % Get the displacement on each node and direction
        q_e = zeros(6, 1);
        q_e(1, 1) = Q(2*n1-1, 1);
        q_e(2, 1) = Q(2*n1, 1);
        q_e(3, 1) = Q(2*n2-1, 1);
        q_e(4, 1) = Q(2*n2, 1);
        q_e(5, 1) = Q(2*n3-1, 1);
        q_e(6, 1) = Q(2*n3, 1);
        sigma(e_dof, :) = ElementStress("tri2d_stress", e_nodes, D, q_e)';
        disp('sigma_e')
        disp(sigma(e_dof, :))
    end

end

⑫ ElementStress.m

求解单元的应变

%% Copyright: YunHao Liu
% Calculate the stress
function sigma = ElementStress(element_type, nodes, D, q_e)
    switch element_type
        case "tri2d_stress"
            % Define Nodes
            x1 = nodes(1, 1);
            y1 = nodes(1, 2);
            x2 = nodes(2, 1);
            y2 = nodes(2, 2);
            x3 = nodes(3, 1);
            y3 = nodes(3, 2);
            % Define Jacobi
            Jacobi = [
                x1 - x3, y1 - y3;
                x2 - x3, y2 - y3;
            ];
            det_Jacobi = det(Jacobi);
            % Define B
            B = [
                y2-y3, 0, y3-y1, 0, y1-y2, 0;
                0, x3-x2, 0, x1-x3, 0, x2-x1;
                x3-x2, y2-y3, x1-x3, y3-y1, x2-x1, y1-y2;
            ]/abs(det_Jacobi);
            % Calculate Stress
            sigma = D.value * B * q_e;
    end

end

⑬ GenerateInpFromNumSeeds.m

根据网格密度信息自动划分网格，并输出信息到inp文件。

%% Copyright: YunHao Liu
% Generate the Inp file with the number of seeds on specific side
function GenerateInpFromNumSeeds(num_of_seeds)
    %% Variable pre-define    
    % The Geometry is Fixed with this Function
    path = 'inp.txt';
    side_len = 1000.0/num_of_seeds;
    thickness = 10.0;
    n_nodes = (3*num_of_seeds+2)*(num_of_seeds+1)/2;
    n_elements = 3*num_of_seeds^2;
    % nodes = NaN * zeros(n_nodes, 2);
    % elements = NaN * zeros(n_elements, 3);
    
    f = fopen("inp.txt", 'w');

    %% Write Config Information
    fprintf(f, "* Config\ndimension, 2\n");
    fprintf(f, "\n");

    %% Write Material Information
    fprintf(f, "* Material\n210e3, 0.2\n");
    fprintf(f, "\n");
    
    %% Write Nodes Information while Loop in Geometry
    fprintf(f, "* Nodes\n%d\n", n_nodes);
    n_dof = 1;
    for row = 1:num_of_seeds+1
        for column = 1 : num_of_seeds + row
            x = side_len * (column-1);
            y = -side_len * (row-1);
            fprintf(f, "%d, %f, %f\n", n_dof, x, y);
            n_dof = n_dof + 1;
        end
    end
    fprintf(f, "\n");

    %% Write Elements Information while Loop in Geometry
    fprintf(f, "* Elements\ntri2d_stress, %d, %f\n", n_elements, thickness);
    n_dof = 1;
    e_dof = 1;
    for row = 1:num_of_seeds
        for column = 1 : num_of_seeds + row
            % Lower Tri
            fprintf(f,'%d, %d, %d, %d\n', e_dof,...
                 n_dof, n_dof+num_of_seeds+row, n_dof+num_of_seeds+row+1);
            e_dof = e_dof + 1;
            % Upper Tri
            if column ~= num_of_seeds + row
                fprintf(f,'%d, %d, %d, %d\n', ...
                    e_dof, n_dof, n_dof+num_of_seeds+row+1, n_dof+1);
                e_dof = e_dof + 1;
            end

            n_dof = n_dof + 1;
        end
    end
    fprintf(f,'\n');

    %% Write Force Information
    fprintf(f, "* Forces\n%d\n%d, %f, %f\n", 1, 1, 1e4, 0.0);
    fprintf(f, "\n");

    %% Write Boundary Conditions
    fprintf(f, "* Boundaries\n%d\n", 3);
    fprintf(f,"%d, all\n", n_nodes);
    fprintf(f,"%d, all\n", n_nodes-num_of_seeds);
    fprintf(f,"%d, all\n", n_nodes-2*num_of_seeds);
    fprintf(f,"\n");
end

求解结果展示

先按照题意，划分为3个单元求解。

变形前后的结构简图如下所示（位移已放大3000倍），蓝色实线为原结构，红色虚线为变形后的结构。

现每个短边上布置6个节点，重新求解，结构变形简图如下所示。应力已用颜色大小表示，每个单元的应力大小用Von Mises应力度量。

顺便给出怎么用 应力分量 去求解 Von Mises 应力的公式，在塑性力学假设下它只与应力偏量有关，这部分的知识可以在塑性力学中找到。其中${\sigma }_{ij}$为应力张量，$S_{ij}$为应力张量的偏张量部分
$$\begin{array}{rl}{S}_{ij}& ={\sigma }_{ij}-\frac{1}{3}{\sigma }_{tt}\\ {\sigma }_i& =\sqrt{3J_2}=\sqrt{\frac{3}{2}{S}_{ij}{S}_{ij}}\\ & =\sqrt{\frac{3}{2}}\cdot \sqrt{S_x^2+S_y^2+S_z^2-2(S_{xy}^2+S_{yz}^2+S_{zx}^2)}\end{array}$$

20节点的情况如下图所示。

ABAQUS有限元模拟与结果对比

ABAQUS有限元只模拟题目中示意的网格划分情况，求解的结果云图如下所示

原始数据（节点位移、单元Mises应力）：

mise_matlab = [
    1.76399629932817,
    1.13143429427077,
    1.04806787555210
]

mise_abaqus = [
    1.764,
    1.13143,
    1.04806
]

u_matlab = [
    0.0113836526324610,
    0.00183434582386225,
    0.00688774908394240,
    -0.000467582400313694,
    1.75726182710415e-08,
    2.28134362652684e-08,
    1.47466737515054e-08,
    -1.29738113060470e-08,
    3.34030828479390e-10,
    -9.83962495922142e-09
]

u_abaqus = [
    0.0113836,
    0.00183433,
    0.00688771,
    -0.000467573,
    5.38159e-33,
    6.98662e-33,
    4.51613e-33,
    -3.97325e-33,
    1.02282e-34,
    -3.01338e-33
]

对比5个节点的10个位移值，绘制分析对比图如下所示。

对比3个单元的Mises应力的求解结果，绘制分析对比图如下图所示。

对比绝对误差（以abaqus模拟结果作为真值）。

位移的绝对误差图如下。

Mises应力的绝对误差图如下。

可以看到结果吻合十分良好。