Abaqus对无位移约束的节点进行反力的提取-NFORC-freebody symbol详解

---

前言

在很多情况下要对节点的力进行提取，特别是对边界上的节点进行力的提取。然而如果abaqus没有对节点施加位移约束的情况下，那么常规的RF的提取方式提取到的反力均为0。

本文以对一个shell单元类型的TPMS施加周期性边界条件进行力学均匀化的过程为例讲解。周期性边界条件是通过节点间用equation限制位移关系来实现变形的，而在load模块中没有对任何实体节点施加位移边界条件（仅对参考点施加了位移）。因而可以看到RF的云图全为0值，如下图

NFORC提取节点力

此时应当打开场输出中的NFORC而非RF，如下图

abaqus帮助文档对nforc的解释为 “Forces at the nodes of the element caused by the stress in the element (internal forces in the global coordinate system).” 表明这个力是根据单元的应力插值到节点上得到的。有一些视频或者博客例如视频1和博客2对这个进行了一定的解释。

对于shell单元NFORC有六个分量（123表示xyz三个方向的力，456表示shell的扭转）结果如下：

注意事项1：云图上提取的NFORC并非真实的节点力

如果是为了计算平面上的合力，那么需要对平面上所有节点的力进行求和，然后除以面积，就是平面上的平均应力。
nforc在用abaqus cae对odb进行显示的时候，显示的是节点的nforc的一个平均值。假设节点a同时被3个单元共享，那么云图上显示的大概是节点a在单元1、2、3上各自的分力按照一定规则进行平均后的结果。
然而实际上我想求的是节点上的等效合力，所以是该节点提取出来的nforc乘以被共享的单元个数。这种在cae中不太方便实现，我写了python脚本来提（pythonset.txt是我想要提取的节点的集set编号），例如：

# -*- coding: mbcs -*-
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
import math
import os
import numpy as np
import xyPlot


#odb_path = 'D:/australia2024/TPMS/tpms-matlabcode/HJPminisurf/test/main1.odb'
odb_path = os.path.abspath(__doc__)+'/' + 'main1.odb'
myodb = openOdb(path=odb_path)
lastframe=myodb.steps['Step-1'].frames[-1]
rf1=lastframe.fieldOutputs['NFORC1']
rf2=lastframe.fieldOutputs['NFORC2']
rf3=lastframe.fieldOutputs['NFORC3']
file_path=odb_path = os.path.abspath(__doc__)+'/' + 'pythonset.txt'
file = open(file_path, 'r')
infos = file.readlines()
file.close()
setid = [int(line.strip()) for line in infos]
setRF=[]
for tmpid in setid:
    tmpregion = myodb.rootAssembly.nodeSets['B-'+str(tmpid)]
    setrf1 = rf1.getSubset(region=tmpregion).values
    setrf1 = [item.data for item in setrf1]
    setrf1=sum(setrf1)
    setrf2 = rf2.getSubset(region=tmpregion).values
    setrf2 = [item.data for item in setrf2]
    setrf2=sum(setrf2)
    setrf3 = rf3.getSubset(region=tmpregion).values
    setrf3 = [item.data for item in setrf3]
    setrf3=sum(setrf3)
    setRF.append([tmpid,setrf1,setrf2,setrf3])
# resultfile = open(os.path.abspath(__doc__)+'/'+'result1.txt','w')
# np.savetxt(resultfile, np.array(setRF), fmt='%d %.9f %.9f %.9f', delimiter=' ')
file=os.path.abspath(__doc__)+'/'+'result1.txt'
with open(file, 'w') as f:
    for row in setRF:
        row_str = ' '.join(map(str, row))
        f.write(row_str + '\n')
f.close()
myodb.close()
print('finish')
# print(np.array(setRF))

注意事项2：NFORC的方向和我们设想的相反

在对TPMS进行均匀化处理时，其实是把TPMS单胞等效成一个均质的立方体。我对TPMS沿x方向进行拉伸位移加载。然后我对面法向分别为xyz轴正方向的3个面上的所有节点力分别提取出来，并用matlab求和除以面积得到如下的应力

结果如下：

可以发现sigma11竟然（其他几个力同理）为负数，这显然不可能的。查了一下资料，发现nforc提取力的方向和表面力的方向是相反的。freebody的方向与nforc相反，所以其方向才是我们想要的参考Producing a symbol plot of free body nodal forces中的介绍：

”You can create a symbol plot that shows the summation of free body nodal forces at individual nodes in your model. The symbol field output variable FREEBODY is the sum of the individual NFORCn components at a node, but with the vector’s direction reversed. The FREEBODY output variable also takes into account the displayed elements by default; in contrast, contour plots of NFORC1 consider the whole model.“

我们可以在cae中进行检查free body symbol，如下图。可以看到法向为x正方向的面上的freebody的x方向力与我们想要的一致（从而说明NFORC的方向和我们设想的相反）

修改后的matlab程序如下（result1.txt为python代码提取后的输出文件）：

data=table2array(readtable('.\test\result1.txt'));
%第一列为set的编号，后三列对应的Fx、Fy、Fz
pfaceXplus;%BB1C1C面上的所有集合编号，包含点线面
pfaceYplus;%ABCD面上的所有集合编号，包含点线面
pfaceZplus;%AA1B1B面上的所有集合编号，包含点线面
ylength = 1; % 你需要的值
xlength = 1; % 你需要的值
zlength = 1; % 你需要的值

% 查找索引列表对应的行并提取对应的行
[is_found, idx] = ismember(pfaceXplus, data(:, 1));
forcefaceXplus = data(idx(is_found), :);
forcefaceXplus=forcefaceXplus(:,2:4);
[is_found, idx] = ismember(pfaceYplus, data(:, 1));
forcefaceYplus = data(idx(is_found), :);
forcefaceYplus=forcefaceYplus(:,2:4);
[is_found, idx] = ismember(pfaceZplus, data(:, 1));
forcefaceZplus = data(idx(is_found), :);
forcefaceZplus=forcefaceZplus(:,2:4);
%加负号的原因是nforc的结果和free body的方向相反
sigma11=-sum(forcefaceXplus(:,1))/ylength/zlength;
sigma12=-sum(forcefaceXplus(:,2))/ylength/zlength;
sigma13=-sum(forcefaceXplus(:,3))/ylength/zlength;
sigma21=-sum(forcefaceYplus(:,1))/xlength/zlength;
sigma22=-sum(forcefaceYplus(:,2))/xlength/zlength;
sigma23=-sum(forcefaceYplus(:,3))/xlength/zlength;
sigma31=-sum(forcefaceZplus(:,1))/xlength/ylength;
sigma32=-sum(forcefaceZplus(:,2))/xlength/ylength;
sigma33=-sum(forcefaceZplus(:,3))/xlength/ylength;

额外说明

可以看到abaqus在对TPMS施加周期性边界条件后进行均匀化弹性张量的计算过程中，结果确实有
sigma12=sigma21这种对称性，所以我们选择一个面按照约定的法则进行提取分量即可