## 位移法解平面桁架

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

$$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}=\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)

$$\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)

$$\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)

$$\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)

$$\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)

$$\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)

$$\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)

$$\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)

$$\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)

$$\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)

// ================================================================
//       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;
}


单元数 节点数 集中荷载数 约束数量

......

......

......

......