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

见贤思齐

Abaqus对无位移约束的节点进行反力的提取-NFORC-freebodysymbol详解前言NFORC提取节点力注意事项1：云图上提取的NFORC并非真实的节点力注意事项2：NFORC的方向和我们设想的相反额外说明前言在很多情况下要对节点的力进行提取，特别是对边界上的节点进行力的提取。然而如果abaqus没有对节点施加位移约束的情况下，那么常规的RF的提取方式提取到的反力均为0。本文以对一个shell单元类型的TPMS施加周期性边界条件进行力学均匀化的过程为例讲解。周期性边界条件是通过节点间用equation限制位移关系来实现变形的，而在load模块中没有对任何实体节点施加位移边界条件（仅对参考点施加了位移）。因而可以看到RF的云图全为0值，如下图NFORC提取节点力此时应当打开场输出中的NFORC而非RF，如下图abaqus帮助文档对nforc的解释为“Forcesatthenodesoftheelementcausedbythestressintheelement(internalforcesintheglobalcoordinatesystem).”表明这个力是根据单元的应力插值到节点上得到的。有一些视频或者博客例如视频1和博客2对这个进行了一定的解释。对于shell单元NFORC有六个分量（123表示xyz三个方向的力，456表示shell的扭转）结果如下：注意事项1：云图上提取的NFORC并非真实的节点力如果是为了计算平面上的合力，那么需要对平面上所有节点的力进行求和，然后除以面积，就是平面上的平均应力。nforc在用abaquscae对odb进行显示的时候，显示的是节点的nforc的一个平均值。假设节点a同时被3个单元共享，那么云图上显示的大概是节点a在单元1、2、3上各自的分力按照一定规则进行平均后的结果。然而实际上我想求的是节点上的等效合力，所以是该节点提取出来的nforc乘以被共享的单元个数。这种在cae中不太方便实现，我写了python脚本来提（pythonset.txt是我想要提取的节点的集set编号），例如：#-*-coding:mbcs-*-frompartimport*frommaterialimport*fromsectionimport*fromassemblyimport*fromstepimport*frominteractionimport*fromloadimport*frommeshimport*fromoptimizationimport*fromjobimport*fromsketchimport*fromvisualizationimport*fromconnectorBehaviorimport*importmathimportosimportnumpyasnpimportxyPlot#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())forlineininfos]setRF=[]fortmpidinsetid:tmpregion=myodb.rootAssembly.nodeSets['B-'+str(tmpid)]setrf1=rf1.getSubset(region=tmpregion).valuessetrf1=[item.dataforiteminsetrf1]setrf1=sum(setrf1)setrf2=rf2.getSubset(region=tmpregion).valuessetrf2=[item.dataforiteminsetrf2]setrf2=sum(setrf2)setrf3=rf3.getSubset(region=tmpregion).valuessetrf3=[item.dataforiteminsetrf3]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'withopen(file,'w')asf:forrowinsetRF:row_str=''.join(map(str,row))f.write(row_str+'\n')f.close()myodb.close()print('finish')#print(np.array(setRF))1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556注意事项2：NFORC的方向和我们设想的相反在对TPMS进行均匀化处理时，其实是把TPMS单胞等效成一个均质的立方体。我对TPMS沿x方向进行拉伸位移加载。然后我对面法向分别为xyz轴正方向的3个面上的所有节点力分别提取出来，并用matlab求和除以面积得到如下的应力结果如下：可以发现sigma11竟然（其他几个力同理）为负数，这显然不可能的。查了一下资料，发现nforc提取力的方向和表面力的方向是相反的。freebody的方向与nforc相反，所以其方向才是我们想要的参考Producingasymbolplotoffreebodynodalforces中的介绍：”Youcancreateasymbolplotthatshowsthesummationoffreebodynodalforcesatindividualnodesinyourmodel.ThesymbolfieldoutputvariableFREEBODYisthesumoftheindividualNFORCncomponentsatanode,butwiththevector’sdirectionreversed.TheFREEBODYoutputvariablealsotakesintoaccountthedisplayedelementsbydefault;incontrast,contourplotsofNFORC1considerthewholemodel.“我们可以在cae中进行检查freebodysymbol，如下图。可以看到法向为x正方向的面上的freebody的x方向力与我们想要的一致（从而说明NFORC的方向和我们设想的相反）修改后的matlab程序如下（result1.txt为python代码提取后的输出文件）：data=table2array(readtable('.\test\result1.txt'));%第一列为set的编号，后三列对应的Fx、Fy、FzpfaceXplus;%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的结果和freebody的方向相反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;1234567891011121314151617181920212223242526272829额外说明可以看到abaqus在对TPMS施加周期性边界条件后进行均匀化弹性张量的计算过程中，结果确实有sigma12=sigma21这种对称性，所以我们选择一个面按照约定的法则进行提取分量即可