MATLAB 向量化实现：三维实体单元表面压强如何等效为节点荷载？

前两期的推文介绍了三维单元表面压强荷载如何等效为节点荷载的理论部分和Abaqus的一些相关Python操作。

本期推文将着重介绍如何将前期分享的理论串起来，凝结于程序中，使得读者更加具体的明白理论在程序中的具体实现过程。

数值实现

单个单元

以C3D20单元为例，代码参考结构分析的有限元法与 MATLAB 程序设计^[1]

function edf = EquivalentDistPressure(node, pressure, gNode)
%   计算分布压力的等效节点力
%   输入参数:
%      node ----------  结点号
%      pressure ------  跟结点号对应的压力值
%      gNode ---------  所有节点的坐标矩阵
%   返回值:
%      edf -----------  等效节点力向量

    x = gNode(node, 1);
    y = gNode(node, 2);
    z = gNode(node, 3);

    edf = zeros(8*3, 1);
    ngp = 3;
    elemType = 'Q8';
    samp = gauss(ngp);
    
    for ig = 1:ngp
        wi = samp(ig,2);
        for jg = 1:ngp
            wj = samp(jg,2);
            [N, NDerivatives] = shapeFunctionQuad(samp, ig, jg, elemType); 
            N_xi = NDerivatives(:,1)';
            N_eta = NDerivatives(:,2)';
            
            x_xi = N_xi * x;y_xi = N_xi * y;z_xi = N_xi * z;
            x_eta = N_eta * x;y_eta = N_eta * y;z_eta = N_eta * z;
            
            px = y_xi * z_eta - y_eta * z_xi;
            py = z_xi * x_eta - z_eta * x_xi;
            pz = x_xi * y_eta - x_eta * y_xi;
            
            sig = pressure * N;  
            
            I = eye(3);
            edf = edf + wi * wj * kron(N,I) * sig * [px; py; pz];
        end
    end
end

根据前期推文的方法索引到这些节点的顺序，写入函数 EquivalentDistPressure 的形参中：

gNode = [
     5.0,  5.0,  7.5;  % 节点 1
     5.0,  0.0,  7.5;  % 节点 2
     5.0,  5.0,  0.0;  % 节点 3
     5.0,  0.0,  0.0;  % 节点 4
     0.0,  5.0,  7.5;  % 节点 5
     0.0,  0.0,  7.5;  % 节点 6
     0.0,  5.0,  0.0;  % 节点 7
     0.0,  0.0,  0.0;  % 节点 8
     0.0,  5.0,  3.75; % 节点 9
     0.0,  2.5,  0.0;  % 节点 10
     0.0,  0.0,  3.75; % 节点 11
     0.0,  2.5,  7.5;  % 节点 12
     5.0,  2.5,  7.5;  % 节点 13
     5.0,  0.0,  3.75; % 节点 14
     5.0,  2.5,  0.0;  % 节点 15
     5.0,  5.0,  3.75; % 节点 16
     2.5,  0.0,  7.5;  % 节点 17
     2.5,  5.0,  7.5;  % 节点 18
     2.5,  0.0,  0.0;  % 节点 19
     2.5,  5.0,  0.0   % 节点 20
];

node = [1, 3, 7, 5, 16, 20, 9, 18];
pressure =  -ones(1, 8);  % 假设施加均匀单位压力
edf = EquivalentDistPressure(node, pressure, gNode);

for i = 1:length(node)
    node_number = node(i);
    force_x = edf((node_number-1)*3 + 1);
    force_y = edf((node_number-1)*3 + 2);
    force_z = edf((node_number-1)*3 + 3);
    
    fprintf('Node %2d: Fx = %8.4f, Fy = %8.4f, Fz = %8.4f\n', node_number, force_x, force_y, force_z);
end

运行结果：

Node  1: Fx =  -0.0000, Fy =   3.1250, Fz =  -0.0000
Node  3: Fx =  -0.0000, Fy =   3.1250, Fz =   0.0000
Node  7: Fx =   0.0000, Fy =   3.1250, Fz =   0.0000
Node  5: Fx =   0.0000, Fy =   3.1250, Fz =  -0.0000
Node 16: Fx =  -0.0000, Fy = -12.5000, Fz =   0.0000
Node 20: Fx =  -0.0000, Fy = -12.5000, Fz =   0.0000
Node  9: Fx =  -0.0000, Fy = -12.5000, Fz =   0.0000
Node 18: Fx =  -0.0000, Fy = -12.5000, Fz =   0.0000

这时可以验证数值的正确性，将等效的节点累加起来，是否等于压强值* 表面面积。

还可以使用 EMM 插件提取 Abaqus 中单元的荷载向量，得到：

这时发现数值上都可以对上，但是对应的节点编号对应不上，这是为什么呢？

细心意向，这里的节点编号是单元内部的局部节点编号，映射在全局节点的位置，还需要看 INP 文件的单元节点连接关系：

*Element, type=C3D20
1,  5,  6,  8,  7,  1,  2,  4,  3, 12, 11, 10,  9, 13, 14, 15,
   16, 18, 17, 19, 20

将之对应起来后，便可与函数 EquivalentDistPressure 计算得到的等效节点荷载对上。

多个单元

相同的道理，可以将之拓展至多个单元情况：原来的 Python 脚本只需要更改对应的模型名字、实例名字、surface 名字即可，核心语句不需要变动，即可输出：

Element 1 face nodes: [1, 4, 13, 10, 35, 39, 28, 37]
Element 3 face nodes: [4, 7, 16, 13, 53, 55, 48, 39]
Element 5 face nodes: [10, 13, 22, 19, 28, 68, 61, 66]
Element 7 face nodes: [13, 16, 25, 22, 48, 78, 74, 68]

将上面罗列的节点顺序，写入函数 EquivalentDistPressure 即可。

曲面边界

模型的曲面边界上承受压强荷载，通过本节的函数 EquivalentDistPressure 也同样可以计算出等效的节点荷载。

我们索引出包含 surfaces 的单元后可以编制以下脚本：

% 初始化全局等效节点力向量
num_nodes = size(gNode, 1);
global_edf = zeros(num_nodes * 3, 1);  % 每个节点有3个分量

% 对每个单元面进行循环
for i = 1:size(element_faces, 1)
    node = element_faces(i, :);  % 当前单元面的节点编号
    
    % 计算当前单元面的等效节点力
    local_edf = EquivalentDistPressure(node, pressure, gNode);
    
    % 将局部等效力累加到全局力向量中
    % 生成全局和局部索引
    node_indices = (node - 1) * DOF;  % 每个节点对应的起始索引
    
    global_indices = [node_indices + 1; node_indices + 2; node_indices + 3];
    local_indices = (0:length(node)-1) * DOF;  % 局部起始索引
    local_indices = [local_indices + 1; local_indices + 2; local_indices + 3];
    
    % 向量化操作
    global_edf(global_indices(:)) = global_edf(global_indices(:)) + local_edf(local_indices(:));

end

对每个单元进行循环，将节点荷载累加至全局节点荷载中，最终得到的全局节点荷载与 EMM 插件提取到的全局节点荷载向量是一致的。可以将两个数组变量进行一个差异化分析，程序中的全局等效节点荷载变量为 global_edf EMM 插件得到的全局荷载变量为 globalLoad：

tolerance = 1e-12;
isEqual = all(abs(global_edf - globalLoad) < tolerance);
disp(isEqual);  % 如果在容差范围内，返回 true；否则返回 false

最终输出为 1，表明两个数组内的元素相差均在 1e-12 范围内，验证了本节的方法得到的表面节点荷载与 Abaqus 处理表面压强的方法得到额表面节点荷载精度一致。

本次推文的单元类型为 C3D20，读者也可用相同的道理拓展至 C3D8 单元或者其他的三维单元类型。

注：EMM插件为 团队在技 术邻公布的Abaqus插件，可以一键输出刚度矩阵和质量矩阵以及载荷向量至Matlab中，使得用户更能直观地感受矩阵的形式，便于和自研程序做对比。下载地址：https://www.jishulink.com/post/341364