位移法解平面桁架

最近稍稍闲下来一点,准备系统的学习一下有限元。
一般认为,矩阵位移法是有限元的雏形,并且矩阵位移法与杆系有限元在操作步骤上几乎完全一样。
这里以下图所示桁架为例简单写一下推导过程,详细的教材上都有。

1

单元编号 单元节点编号 节点局部编号
1 1,3 i,j
2 1,2 i,j
3 3,4 i,j
4 2,4 i,j
5 3,2 i,j

一.单元刚度方程组(矩阵)
记杆件\(m\)\(x\)轴夹角为\(\alpha\),则

\(f_{i}^{m}=\frac{{E^m}{A^m}}{L_{0}^{m}}{\left[{\left(u_i-u_j\right)}{cos\alpha}+{\left(v_i-v_j\right)sin\alpha} \right]}\)

因为\(f_{ix}^{m}=f_{i}^{m}cos\alpha\),\(f_{iy}^{m}=f_{i}^{m}sin\alpha\),所以

\(f_{ix}^{m}=\frac{{E^m}{A^m}}{L_{0}^{m}}{\left[{\left(u_i-u_j\right)}{cos^2\alpha}+{\left(v_i-v_j\right){sin\alpha}{cos\alpha}}\right]}\) (1)

\(f_{iy}^{m}=\frac{{E^m}{A^m}}{L_{0}^{m}}{\left[{\left(u_i-u_j\right)}{cos\alpha}{sin\alpha}+{\left(v_i-v_j\right){sin^2\alpha}}\right]}\) (2)

\(f_{jx}^{m}=-\frac{{E^m}{A^m}}{L_{0}^{m}}{\left[{\left(u_i-u_j\right)}{cos^2\alpha}+{\left(v_i-v_j\right){sin\alpha}{cos\alpha}}\right]}\) (3)

\(f_{jy}^{m}=-\frac{{E^m}{A^m}}{L_{0}^{m}}{\left[{\left(u_i-u_j\right)}{cos\alpha}{sin\alpha}+{\left(v_i-v_j\right){sin^2\alpha}}\right]}\) (4)

将公式(1)~(4)写成矩阵形式,得

\(\frac{{E^m}{A^m}}{{L}_{0}^{m}}\begin{bmatrix}{cos}^2{\alpha} &{sin\alpha}{cos\alpha} &-{cos}^{2}{\alpha} &-{sin\alpha}{cos\alpha}\\ {cos\alpha}{sin\alpha} &{sin}^{2}{\alpha} &-{cos\alpha}{sin\alpha} &-{sin}^2{\alpha}\\ -{cos}^2{\alpha} &-{sin\alpha}{cos\alpha} &-{cos}^{2}{\alpha} &{sin\alpha}{cos\alpha}\\ -{cos\alpha}{sin\alpha} &-{sin}^{2}{\alpha} &{cos\alpha}{sin\alpha} &{sin}^2{\alpha}\\ \end{bmatrix}\begin{Bmatrix}{u}_{i}\\{v}_{i}\\{u}_{j}\\{v}_{j} \end{Bmatrix}=\begin{Bmatrix}{f}_{ix}^{m}\\{f}_{jy}^{m}\\{f}_{jx}^{m}\\{f}_{jy}^{m}\end{Bmatrix}\) (5)

将公式(5)展开为符号形式,得

\(\begin{bmatrix}k_{11}^{m} &k_{12}^{m} &k_{13}^{m} &k_{14}^{m}\\k_{21}^{m} &k_{22}^{m} &k_{23}^{m} &k_{24}^{m}\\k_{31}^{m} &k_{32}^{m} &k_{33}^{m} &k_{34}^{m}\\k_{41}^{m}&k_{42}^{m}&k_{43}^{m}&k_{44}^{m}\\\end{bmatrix}\begin{Bmatrix}u_i\\v_i\\u_j\\v_j\end{Bmatrix}=\begin{Bmatrix}f_{ix}^{m}\\f_{iy}^{m}\\f_{jx}^{m}\\f_{jy}^{m}\end{Bmatrix}\)(6)

将表格中单元与节点对应关系填入公式 (6),得

单元1,m=1,i=1,j=3
\(\begin{bmatrix}k_{11}^{1}&k_{12}^{1}&k_{13}^{1}&k_{14}^{1}\\k_{21}^{1}&k_{22}^{1}&k_{23}^{1}&k_{24}^{1}\\k_{31}^{1}&k_{32}^{1}&k_{33}^{1}&k_{34}^{1}\\k_{41}^{1}&k_{42}^{1}&k_{43}^{1}&k_{44}^{1}\\\end{bmatrix}\begin{Bmatrix}u_1\\v_1\\u_3\\v_3\end{Bmatrix}=\begin{Bmatrix}f_{1x}^{1}\\f_{1y}^{1}\\f_{3x}^{1}\\f_{3y}^{1}\end{Bmatrix}\)(7)

单元2,m=1,i=1,j=2
\(\begin{bmatrix}k_{11}^{2}&k_{12}^{2}&k_{13}^{2}&k_{14}^{2}\\k_{21}^{2}&k_{22}^{2}&k_{23}^{2}&k_{24}^{2}\\k_{31}^{2}&k_{32}^{2}&k_{33}^{2}&k_{34}^{2}\\k_{41}^{2}&k_{42}^{2}&k_{43}^{2}&k_{44}^{2}\\\end{bmatrix}\begin{Bmatrix}u_1\\v_1\\u_2\\v_2\end{Bmatrix}=\begin{Bmatrix}f_{1x}^{2}\\f_{1y}^{2}\\f_{2x}^{2}\\f_{2y}^{2}\end{Bmatrix}\)(8)

单元3,m=1,i=3,j=4
\(\begin{bmatrix}k_{11}^{3}&k_{12}^{3}&k_{13}^{3}&k_{14}^{3}\\k_{21}^{3}&k_{22}^{3}&k_{23}^{3}&k_{24}^{3}\\k_{31}^{3}&k_{32}^{3}&k_{33}^{3}&k_{34}^{3}\\k_{41}^{3}&k_{42}^{3}&k_{43}^{3}&k_{44}^{3}\\\end{bmatrix}\begin{Bmatrix}u_3\\v_3\\u_4\\v_4\end{Bmatrix}=\begin{Bmatrix}f_{3x}^{3}\\f_{3y}^{3}\\f_{4x}^{3}\\f_{4y}^{3}\end{Bmatrix}\)(9)

单元4,m=1,i=2,j=4
\(\begin{bmatrix}k_{11}^{4}&k_{12}^{4}&k_{13}^{4}&k_{14}^{4}\\k_{21}^{4}&k_{22}^{4}&k_{23}^{4}&k_{24}^{4}\\k_{31}^{4}&k_{32}^{4}&k_{33}^{4}&k_{34}^{4}\\k_{41}^{4}&k_{42}^{4}&k_{43}^{4}&k_{44}^{4}\\\end{bmatrix}\begin{Bmatrix}u_2\\v_2\\u_4\\v_4\end{Bmatrix}=\begin{Bmatrix}f_{2x}^{4}\\f_{2y}^{4}\\f_{4x}^{4}\\f_{4y}^{4}\end{Bmatrix}\)(10)

单元5,m=1,i=3,j=2
\(\begin{bmatrix}k_{11}^{5}&k_{12}^{5}&k_{13}^{5}&k_{14}^{5}\\k_{21}^{5}&k_{22}^{5}&k_{23}^{5}&k_{24}^{5}\\k_{31}^{5}&k_{32}^{5}&k_{33}^{5}&k_{34}^{5}\\k_{41}^{5}&k_{42}^{5}&k_{43}^{5}&k_{44}^{5}\\\end{bmatrix}\begin{Bmatrix}u_3\\v_3\\u_2\\v_2\end{Bmatrix}=\begin{Bmatrix}f_{3x}^{5}\\f_{3y}^{5}\\f_{2x}^{5}\\f_{2y}^{5}\end{Bmatrix}\)(11)

二.总刚度方程组(矩阵)

将(7)~(11)矩阵相加合并,可得
\(\begin{bmatrix}k_{11} &k_{12} &k_{13} &k_{14} &k_{15} &k_{16} &k_{17} &k_{18}\\k_{21} &k_{22} &k_{23} &k_{24} &k_{25} &k_{26} &k_{27} &k_{28}\\k_{31} &k_{32} &k_{33} &k_{34} &k_{35} &k_{36} &k_{37} &k_{38}\\k_{41} &k_{42} &k_{43} &k_{44} &k_{45} &k_{46} &k_{47} &k_{48}\\k_{51} &k_{52} &k_{53} &k_{54} &k_{55} &k_{56} &k_{57} &k_{58}\\k_{61} &k_{62} &k_{63} &k_{64} &k_{65} &k_{66} &k_{67} &k_{68}\\k_{71} &k_{72} &k_{73} &k_{74} &k_{75} &k_{76} &k_{77} &k_{78}\\k_{81} &k_{82} &k_{83} &k_{84} &k_{85} &k_{86} &k_{87} &k_{88}\\\end{bmatrix}\begin{Bmatrix}u_{1}\\v_{1}\\u_{2}\\v_{2}\\u_{3}\\v_{3}\\u_{4}\\v_{4}\end{Bmatrix}=\begin{Bmatrix}F_{1x}\\F_{1y}\\F_{2x}\\F_{2y}\\F_{3x}\\F_{3y}\\F_{4x}\\F_{4y}\end{Bmatrix}\)(12)

三.求解

引入约束条件及外力
\(\begin{bmatrix}k_{11} &k_{12} &k_{13} &k_{14} &k_{15} &k_{16} &k_{17} &k_{18}\\k_{21} &k_{22} &k_{23} &k_{24} &k_{25} &k_{26} &k_{27} &k_{28}\\k_{31} &k_{32} &k_{33} &k_{34} &k_{35} &k_{36} &k_{37} &k_{38}\\k_{41} &k_{42} &k_{43} &k_{44} &k_{45} &k_{46} &k_{47} &k_{48}\\k_{51} &k_{52} &k_{53} &k_{54} &k_{55} &k_{56} &k_{57} &k_{58}\\k_{61} &k_{62} &k_{63} &k_{64} &k_{65} &k_{66} &k_{67} &k_{68}\\k_{71} &k_{72} &k_{73} &k_{74} &k_{75} &k_{76} &k_{77} &k_{78}\\k_{81} &k_{82} &k_{83} &k_{84} &k_{85} &k_{86} &k_{87} &k_{88}\\\end{bmatrix}\begin{Bmatrix}u_{1}\\v_{1}\\u_{2}\\v_{2}\\0\\0\\u_{4}\\0\end{Bmatrix}=\begin{Bmatrix}0\\F\\0\\0\\R_{1}\\R_{2}\\0\\R_{3}\end{Bmatrix}\)(13)
其中\(F\)已知,\(R_1,R_2,R_3\)未知。

在边界条件引入后,方程组所包含的未知数与方程个数相等,方程组有解。
对于方程组(13)的求解,可先对如下子矩阵求解,从而得到位移向量。
\(\begin{bmatrix}k_{11} &k_{12} &k_{13} &k_{14} &k_{17}\\k_{21} &k_{22} &k_{23} &k_{24} &k_{27}\\k_{31} &k_{32} &k_{33} &k_{34} &k_{37}\\k_{41} &k_{42} &k_{43} &k_{44} &k_{47}\\k_{71} &k_{72} &k_{73} &k_{74} &k_{77}\end{bmatrix}\begin{Bmatrix}u_{1}\\v_{1}\\u_{2}\\v_{2}\\u_{4}\end{Bmatrix}=\begin{Bmatrix}0\\F\\0\\0\\0\end{Bmatrix}\)(14)

三.Eigen矩阵库
由于编写程序需要大量进行矩阵运算,这里采用Eigen库来简化代码编制过程。Eigen是一个轻量级的矩阵库,其他的矩阵和向量操作都比较完善,而且速度不错。Eigen库没有任何依赖,从官网下载最新版本后,解压,将其中的Eigen文件夹复制到编译器的include目录下就可以使用,门槛比较低。

四.CODE

// =========================================================================
//       Filename:  Truss.cpp
//    Description:  位移法解平面桁架
//        Version:  1.0
//        Created:  2014-5-22 13:43:01
//        Company:  CCCC First Highway Consultants Co.,LTD
//       Revision:  none
//       Compiler:  g++
//         Author:  whchen (http://blog.whchen.net), 
// =========================================================================
 
#include<stdio.h>
#include<iostream>
#include<string.h>
#include<string>
#include<math.h>
#include<sstream>
#include<fstream>
#include<iomanip>
#include<Eigen/Dense>
 
using namespace Eigen;
using namespace std;
 
typedef struct node
{
  double x,y;
} node;
 
typedef struct frame
{
  //iNum,jNum,E,A,L0,M 分别为节点号i,j,弹性模量,截面积,长度,杆单元刚度矩阵
  int iNum,jNum;
  double E,A,L0;
  Matrix4d M;
} frame;
 
typedef struct restrain
{
  int num;
  int x,y;
} restrain;
 
int main()
{
  string sTmp;
  //IO
  ifstream fileInput("input.txt");
  ofstream fileOutput("output.txt");
 
  //变量定义
  int NE,NJ,NP,NS;//单元数、节点数、节点荷载数,约束节点数
 
  //输入单元节点荷载总数
  getline(fileInput,sTmp);
  istringstream line(sTmp);
  line>>NE>>NJ>>NP>>NS;
 
  node joint[NJ];
  //输入节点信息
  for(int i=0;i<NJ;++i)
	{
	  int k;
	  getline(fileInput,sTmp);
	  istringstream line(sTmp);
	  line>>k;
	  k-=1;
	  line>>joint[k].x>>joint[k].y;
	}
 
  frame elem[NE];
  //输入单元信息
  for(int i=0;i<NE;++i)
	{
	  int k,ii,jj;
	  getline(fileInput,sTmp);
	  istringstream line(sTmp);
	  line>>k;
	  k-=1;
	  line>>ii>>jj>>elem[k].E>>elem[k].A;
	  ii-=1;
	  jj-=1;
	  elem[k].iNum=ii;
	  elem[k].jNum=jj;
	  elem[k].L0=sqrt(pow(joint[jj].x-joint[ii].x,2)+pow(joint[jj].y-joint[ii].y,2));
 
	  double cosA,sinA;
	  cosA=(joint[jj].x-joint[ii].x)/elem[k].L0;
	  sinA=(joint[jj].y-joint[ii].y)/elem[k].L0;
 
	  //计算单元刚度矩阵
	  double aMat[4][4]={
		   cosA*cosA,  sinA*cosA, -cosA*cosA, -sinA*cosA,
		   cosA*sinA,  sinA*sinA, -cosA*sinA, -sinA*sinA,
		  -cosA*cosA, -sinA*cosA,  cosA*cosA,  sinA*cosA,
		  -cosA*sinA, -sinA*sinA,  cosA*sinA,  sinA*sinA
	  };
 
	  for(int i=0;i<4;++i)
		for(int j=0;j<4;++j)
		  elem[k].M(i,j)=aMat[i][j];
 
	  elem[k].M*=(elem[k].E*elem[k].A/elem[k].L0);
	}
 
  //叠加获得总刚度矩阵
  MatrixXd KM=MatrixXd::Zero(2*NJ,2*NJ);
  for(int i=0;i<NE;++i)
	{
	  int ii=elem[i].iNum;
	  int jj=elem[i].jNum;
 
	  KM(2*ii  ,2*ii  )+=elem[i].M(0,0);
	  KM(2*ii  ,2*ii+1)+=elem[i].M(0,1);
	  KM(2*ii+1,2*ii  )+=elem[i].M(1,0);
	  KM(2*ii+1,2*ii+1)+=elem[i].M(1,1);
 
	  KM(2*ii  ,2*jj  )+=elem[i].M(0,2);
	  KM(2*ii  ,2*jj+1)+=elem[i].M(0,3);
	  KM(2*ii+1,2*jj  )+=elem[i].M(1,2);
	  KM(2*ii+1,2*jj+1)+=elem[i].M(1,3);
 
	  KM(2*jj  ,2*ii  )+=elem[i].M(2,0);
	  KM(2*jj  ,2*ii+1)+=elem[i].M(2,1);
	  KM(2*jj+1,2*ii  )+=elem[i].M(3,0);
	  KM(2*jj+1,2*ii+1)+=elem[i].M(3,1);
 
	  KM(2*jj  ,2*jj  )+=elem[i].M(2,2);
	  KM(2*jj  ,2*jj+1)+=elem[i].M(2,3);
	  KM(2*jj+1,2*jj  )+=elem[i].M(3,2);
	  KM(2*jj+1,2*jj+1)+=elem[i].M(3,3);
	}
 
  //定义总外力向量,位移向量
  VectorXd VF=VectorXd::Zero(2*NJ);
  VectorXd VV=VectorXd::Zero(2*NJ);
  bool T[2*NJ];
  memset(T,true,sizeof(T));
 
  //输入荷载
  for(int i=0;i<NP;++i)
	{
	  getline(fileInput,sTmp);
	  istringstream line(sTmp);
	  int num;
	  double Px,Py;
	  line>>num>>Px>>Py;
	  num-=1;
	  VF(2*num)=Px;
	  VF(2*num+1)=Py;
	}
 
  restrain rst[NS];
  //输入约束信息
  for(int i=0;i<NS;++i)
	{
	  getline(fileInput,sTmp);
	  istringstream line(sTmp);
	  line>>rst[i].num>>rst[i].x>>rst[i].y;
	  rst[i].num-=1;
	  VF(2*rst[i].num  )=rst[i].x;
	  if(rst[i].x) T[2*rst[i].num  ]=false;
	  VF(2*rst[i].num+1)=rst[i].y;
	  if(rst[i].y) T[2*rst[i].num+1]=false;
	}
 
  //开始写求解部分
  //提取子矩阵,先求位移
  int cnt=0;
  for(int i=0;i<2*NJ;++i)
	if(T[i]) ++cnt;
  MatrixXd subM=MatrixXd::Zero(cnt,cnt);
  int ii,jj;
  ii=jj=0;
  for(int i=0;i<2*NJ;++i)
	for(int j=0;j<2*NJ;++j)
	  if(T[i]&&T[j])
		{
		  subM(ii,jj)=KM(i,j);
		  ++ii;
		  if(ii>=cnt)
			{
			  ii=0;
			  ++jj;
			}
		}
 
  //提取子向量
  VectorXd subVF=VectorXd::Zero(cnt);
  ii=0;
  for(int i=0;i<2*NJ;++i)
	if(T[i]) subVF(ii++)=VF(i);
 
  //求解位移
  VectorXd subVV=subM.colPivHouseholderQr().solve(subVF);
 
  //合并结果进总位移向量
  ii=0;
  for(int i=0;i<2*NJ;++i)
	if(T[i])
	  VV(i)=subVV(ii++);
 
  VF=KM*VV;
 
  //输出节点位移
  fileOutput<<setw(25)<<"节点位移"<<endl;
  fileOutput<<"  Joint"<<"                 X"<<"                 Y"<<endl;
  for(int i=0;i<2*NJ;i+=2)
	{
	  fileOutput<<setw(5)<<i/2+1<<"    ";
	  fileOutput<<fixed<<setprecision(12)<<setw(15)<<VV(i)
		    <<"    "<<setw(15)<<VV(i+1)<<endl;
	}
 
  //输出单元轴力
  fileOutput<<endl<<endl;
  fileOutput<<setw(25)<<"单元轴力"<<endl;
  fileOutput<<"Element"<<"                 F"<<endl;
  for(int i=0;i<NE;++i)
	{
	  int ii,jj;
	  double E,A,L0;
	  ii=elem[i].iNum;
	  jj=elem[i].jNum;
	  E=elem[i].E;
	  A=elem[i].A;
	  L0=elem[i].L0;
	  double xi,xj,yi,yj,ui,uj,vi,vj;
	  xi=joint[ii].x;
	  yi=joint[ii].y;
	  xj=joint[jj].x;
	  yj=joint[jj].y;
	  ui=VV(2*ii);
	  vi=VV(2*ii+1);
	  uj=VV(2*jj);
	  vj=VV(2*jj+1);
	  double fi=E*A/L0*((xj-xi)*(uj-ui)+(yj-yi)*(vj-vi))/L0;
	  fileOutput<<setw(5)<<i+1<<"    ";
	  fileOutput<<fixed<<setprecision(12)<<setw(15)<<fi<<endl;
	}
 
  //输出支座反力
  fileOutput<<endl<<endl;
  fileOutput<<setw(25)<<"支承反力"<<endl;
  fileOutput<<"  Joint"<<"                FX"<<"                FY"<<endl;
  for(int i=0;i<NS;++i)
	{
	  int num=rst[i].num;
	  fileOutput<<setw(5)<<num+1<<"    ";
	  fileOutput<<fixed<<setprecision(12)<<setw(15)<<VF(2*num)
		    <<"    "<<setw(15)<<VF(2*num+1)<<endl;
	}
  return 0;
}

五.输入数据格式

单元数 节点数 集中荷载数 约束数量
节点1 x坐标 y坐标
节点2 x坐标 y坐标
节点3 x坐标 y坐标
节点4 x坐标 y坐标
节点5 x坐标 y坐标
......
单元1 i节点 j节点 弹性模量 截面积
单元2 i节点 j节点 弹性模量 截面积
单元3 i节点 j节点 弹性模量 截面积
单元4 i节点 j节点 弹性模量 截面积
单元5 i节点 j节点 弹性模量 截面积
......
荷载1 x分量 y分量
荷载2 x分量 y分量
荷载3 x分量 y分量
......
约束1 X是否固定 Y是否固定(1或0)
约束2 X是否固定 Y是否固定(1或0)
约束3 X是否固定 Y是否固定(1或0)
......

六.使用sap2000校核
模型
1
对比结果
1
2
3
4

本文中使用到的全部文件Truss.zip

» 本博客采用署名 2.5 中国大陆许可协议进行许可,本文版权归作者所有,欢迎转载,但必须在明显位置给出原文连接。
anyShare分享到:

Leave a Comment

NOTE - You can use these HTML tags and attributes:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>