有限元模型中非节点外力的处理用C++实现

有限元建模时，在网格划分过程中会尽量在外力和外部约束作用处设置节点。因为有限元求解问题求解的刚度方程F=KD中的位移向量D和力向量F指的都是节点位移和节点力。但是，有时候模型中的外力作用位置处确实不好添加节点，此时就需要采用把有限元模型中的非节点力转化为节点力的技术，即“等效节点力”法。消去模型中的非节点力，之后才能求解刚度方程F=KD。

（1）等效节点力的原理

现有一个分布外荷载作用在杆件单元的非节点上，如下图（a）,我们想办法找出和图（a）中的力对有限元模型的作用效果相同的节点力，如下图（b）。根据图（a）中力的参数求得图（b）中节点力的参数。这样，就可以用图（b）中的节点力代替图（a）中的非节点力，称图（b）中的节点力为等效节点力。

1）等功法（形函数法）

我们可以用等功法求等效节点力。原理为：对于任意节点位移，分布荷载w(x)通过位移场v(x)所做的功，等于节点荷载fiy（剪力）和mi（弯矩）通过节点位移vi（位移）和Φi（转角）所做的功。

分布荷载所做的功为：

式中的v(x)为杆件单元的形函数。

离散节点荷载所做的功为：

令Wdiscrete=Wdistributed，即可得到fiy和mi的值，即为等效节点力。

2）查表法

实际上等效节点力问题是结构力学中经常求解的一类问题。已经有一部分现成的图表可查。如果找不到图表，可以依照等功法计算。对于刚性连接的杆单元还可以通过结构的固端反力推得等效节点力，而常见结构受力的固端反力也有表可查。常见的等效节点力如下表所示：

（2）将输入有限元程序的非节点力按照“等效节点力”方法转化为等效节点力的C++实现

该程序将输入的当前非节点力的外荷载计算，得到等效节点力。由于“等效节点力”是以“单元”为集合计算的，即计算“一次等效节点力”就是计算了“一个单元”的等效节点力。故对于平面杆系结构，计算完成后要用一个长度为6（1个单元有两个节点共6个自由度）的一维数组存储等效节点力。

由于等效节点力是在单元局部坐标下计算的，在累加到总体节点力F向量之前需要将其转换到全局坐标下。

1）前期准备：

a)存储了非节点力的结构体数组pLoad。pLoad结构体中的属性为：当前载荷的载荷类型号、载荷作用方向、载荷值/位移值、载荷作用的单元号、载荷作用的节点号、载荷作用位置或分布长度、等等所有计算等效节点力时可能用到的参数。

b)存储等效节点荷载的一维数组，pFixedEndF[6]

2）程序流程图：

3）C++程序为：

void FixedEndForceCalcu(Element *pElem,Material *pMate,Section *pSect,Load *pLoad,double *pFixedEndF,int i) //将平面结构作为整体时，所受的外力以“等价节点力”的形式转移到单元的节点上
//此处pLoad是从其他函数的（pLoad+i）传进来的，指向载荷数组中的当前荷载！
//pElem指向单元数组；pMate指向材料数组；pSect指向截面数组
//pFixedEndF指向一维数组，是外界传入的
{
	int iMateType,iSectType,iLoadType,iLoadedElem;
	double dE,dAlpha;
	double dA,dIz,dH;
	double dt0,dt1,dBuf;
	double da,dQ,dc,dg,dl,db,ds;

	double **pT= TowArrayDoubAlloc(6,6); //申请了6行6列的矩阵的内存空间
	double **pTT= TowArrayDoubAlloc(6,6);  //申请了6行6列的矩阵的内存空间
	double * pTemp=new double[6];  //申请了有6个元素的一维数组

	double dXi=0;
	double dXj=0;
	double dYi=0; //存储梁端剪力，全局以y方向为梁的断面方向，x方向为梁的延伸方向
	double dYj=0;
	double dMi=0;
	double dMj=0;

	iLoadedElem=pLoad->iLoadedElem;//将当前荷载的载荷作用的单元号放入iLoadedElem中
	iLoadType=pLoad->iType; //将当前荷载的载荷类型号放入iLoadType中

//计算-------------------------------------------------------------------------
//直接作用在杆上的非节点载荷，包括温度载荷，先计算固端力---------------------------------
	da=pLoad->dPosition; //取出当前荷载或力矩的载荷作用位置或分布长度放入da中
	dl=(pElem+iLoadedElem)->dLength; //将当前荷载的作用单元的长度放入dl中
	dQ=pLoad->dValue; //将当前荷载的荷载值或力矩值/位移值放入dQ中
	dc=da/dl; //求荷载作用的长度占杆单元总长度的百分比
	dg=dc*dc;
	db=dl-da; //无载荷作用长度占杆单元总长度的百分比
	switch(iLoadType) //判断当前荷载的荷载类型
	//以下计算的都是“剪力”和“弯矩”
	{
	case LATERAL_FORCE: //当前荷载为梁上横向集中力
	//按有限元方法中，梁单元作用有横向荷载时的等价节点力法将集中外力分配到梁两端节点上
		ds=db/dl; //即，ds=(dl-da)/dl,即，（杆长-荷载在杆上的作用位置或分布长度）/杆长
		dYi=-dQ*ds*ds*(1.0+2.0*dc); //梁一端的等价节点力
		dYj=-dQ*dg*(1.0+2.0*ds); //梁另一端节点上的等价节点力
		dMi=-dQ*ds*ds*da; //梁一端的等价节点弯矩
		dMj=dQ*db*dg; //梁另一端的等价节点弯矩
		break; //退出此次循环
	case LATERAL_UNIFORM_PRESSURE: //当前荷载为梁上均布荷载
		ds=dQ*da*0.5;
		dYi=-ds*(2.0-2.0*dg+dc*dg);  //求节点i的等效剪力
		dYj=-ds*dg*(2.0-dc);
		ds=ds*da/6.0;
		dMi=-ds*(6.0-8.0*dc+3.0*dg); //求节点i的等效弯矩
		dMj=ds*dc*(4.0-3.0*dc);
		break;
	case MOMENT_ON_A_POINT: //当前荷载为梁上集中力矩
		ds=db/dl;
		dYi=-6.0*dQ*dc*ds/dl;
		dYj=-dYi;
		dMi=dQ*ds*(2.0-3.0*ds);
		dMj=dQ*dc*(2.0-3.0*dc);
		break;
	case LATERAL_LINEARLY_PRESSURE:  //当前荷载为梁上三角形分布荷载
		ds=dQ*da*0.25;
		dYi=-ds*(2.0-3.0*dg+1.6*dg*dc);
		dYj=-ds*dg*(3.0-1.6*dc);
		ds*=da;
		dMi=ds*(2.0-3.0*dc+1.2*dg)/1.5;
		dMj=-ds*dc*(1.0-0.5*dc);
		break;
	case AXIAL_PRESSURE:  //当前荷载为梁上轴向分布荷载
		dXi=-dQ*da*(1.0-0.5*dc); //梁的轴力
		dXj=-0.5*dQ*dc*da;
		break;
	case AXIAL_FORCE:  //当前荷载为梁上轴向集中力
		dXi=-dQ*db/dl; //材料力学里的梁的剪力！
		dXj=-dQ*dc;
		break;
	case MOMENT_ON_BEAM: //当前荷载为梁上分布力矩
		ds=db/dl;
		dYi=-dQ*dg*(3.0*ds+dc);
		dYj=-dYi;
		ds*=db/dl;
		dMi=-dQ*ds*da;
		dMj=dQ*dg*db;
		break;
	case TEMPERATURE: //当前荷载为梁上温度荷载
		iMateType=(pElem+iLoadedElem)->iMaterial; //将当前荷载作用的单元的材料在材料数组中的下标放入iMateType中
		dE=(pMate+iMateType)->dE; //将当前单元的材料的弹性模量放入dE中
		dAlpha=(pMate+iMateType)->dAlpha; //将当前材料的线膨胀系数放入dAlpha中
		iSectType=(pElem+iLoadedElem)->iSection; //将当前荷载作用的单元的截面类型数组下标放入iSectType中
		dA=(pSect+iSectType)->dA; //取出当前截面类型的截面积
		dIz=(pSect+iSectType)->dIz; //取出当前截面的横截面惯性矩
		dH=(pSect+iSectType)->dH; //取出当前截面的横截面高
		dt0=pLoad->dT0; //取出当前荷载作用的杆上表面温度变化值
		dt1=pLoad->dT1; //取出当前荷载作用的杆下表面温度变化值
		dBuf=0.5*dAlpha*(dt0+dt1)*dE*dA; //计算梁端剪力
		dXi=dBuf;
		dXj=-dBuf;
		dBuf=dAlpha*(dt1-dt0)*dE*dIz/dH; //计算梁端弯矩
		dMi=dBuf;
		dMj=-dBuf;
		break;
	}

//经反号（说明上面计算等效节点力用的是材料力学的正方向！）得到局部坐标系下的等效节点力----------------------------------------
	VectorZeroize(6,pTemp); //将有6个元素的一维数组pTemp置零
	pTemp[0]=0.0;//轴向等效节点力为0
	pTemp[1]=-dYi; //横向等效节点力
	pTemp[2]=-dMi;//等效节点弯矩
	pTemp[3]=0.0;
	pTemp[4]=-dYj;
	pTemp[5]=-dMj;
//然后将其转换到整体坐标系下--------------------------------------------------
	MatrixZeroize(6,6,pT); //将6行6列的pT矩阵置零
	//构建局部坐标和全局坐标的变换矩阵（平面刚架考虑轴向影响的变换矩阵是6阶）
	pT[2][2]=pT[5][5]=1.0; //
	pT[0][0]=pT[1][1]=pT[3][3]=pT[4][4]=(pElem+iLoadedElem)->dCos; //将当前荷载作用的单元的局部x方向和全局x方向的夹角的余弦取出
	pT[0][1]=pT[3][4]=(pElem+iLoadedElem)->dSin; //将当前荷载作用的单元的局部x方向与全局x方向的夹角的正弦值取出
	pT[1][0]=pT[4][3]=-(pElem+iLoadedElem)->dSin; //取正弦值的相反数

	MatrixTrans(6,6,pT,pTT); //将pT矩阵转置，结果放入pTT矩阵中
	MatrixVectorMultiply(6,6,pTT,pTemp,pFixedEndF); //矩阵pTT*列向量pTemp，得到结果列向量放入pFixedEndF中

//释放存储空间，防止内存泄漏----------------------------------------------------------------
	TowArrayFree(6,pT);
	TowArrayFree(6,pTT);
	delete []pTemp; //释放一维数组，数组长度不必指出
}