c语言计算平面桁架内力计算程序

#include<stdio.h>#include<math.h>#include<stdlib.h>#define M 5int n,nc,nn,m,e,f;//节点总数，固定节点数，自由度数，杆件数int io,jo;//单根杆对号指示数int ihl[M],ihr[M];//杆件左右节点号double a[M];//各杆截面积double mm[M];//杆件质量double ea[M];//杆件EA的值double x[M],y[M];//节点坐标double dp[M];//总体系下的节点载荷double t[2];//0,1分别为坐标转换矩阵的cos(),sin()double c[2][2];//总体系下的单刚double clxy[3];//0,1,2分别为杆长，正弦，余弦double h[M];//杆件轴力double r[M][M];//总刚度阵double rd;//桁架轴力杆局部系单刚double u[M];//桁架节点位移double v[2];//存放节点位移差double d[M];//LDLT分解时的D矩阵的对角线元素double l[M][M];////LDLT分解时的D矩阵的对角线元素double fdp[M];//总体系下支座反力void iojo(int k)//计算对号指示数io,jo{int i,j;i=ihl[k-1];//k号杆左节点号进入ij=ihr[k-1];//k号杆节点右号进入iio=2*(i-nc-1);//uxi前未知位移的个数jo=2*(j-nc-1);//uyi前未知位移的个数}void ch(int k)//计算杆长与方向余弦函数{int i,j;i=ihl[k-1];//k号杆左节点进入ij=ihr[k-1];//k号杆右节点进入jclxy[1]=x[j-1]-x[i-1];//k号杆x坐标差clxy[2]=y[j-1]-y[i-1];//k号杆y坐标差clxy[0]=sqrt(clxy[1]*clxy[1]+clxy[2]*clxy[2]);//k号杆长clxy[1]=clxy[1]/clxy[0];//k号杆件x轴余弦clxy[2]=clxy[2]/clxy[0];//k号杆件y轴余弦void stif(int k)//计算k号杆件总体系下的单元刚度阵{int i,j;ch(k);//调用ch(),计算k号杆件杆长与余弦t[0]=clxy[1];t[1]=clxy[2];rd=ea[k-1]/clxy[0];for(i=0;i<2;i++){for(j=0;j<2;j++){c[i][j]=t[i]*t[j]*rd;}}}void dor()//总体系下的总刚度阵的组集{int i,j,k;for(i=0;i<nn;i++){for(j=0;j<nn;j++){r[i][j]=0.0;//总刚度阵清0}}for(k=1;k<=m;k++){iojo(k);//调用k的对号指示函数，从而确定组集位置stif(k);//第k号杆件的总体系下的单刚if(io>=0){for(i=io+1;i<=io+2;i++){for(j=io+1;j<=io+2;j++){r[i-1][j-1]+=c[i-io-1][j-io-1];//在r中io+1，io+2行以及io+1，io+2列位置累加k的单刚}for(j=jo+1;j<=jo+2;j++){r[i-1][j-1]-=c[i-io-1][j-jo-1];//在r中io+1，io+2行以及jo+1，jo+2列位置累加k的负单刚r[j-1][i-1]=r[i-1][j-1];}}for(i=jo+1;i<=jo+2;i++){for(j=jo+1;j<=jo+2;j++){r[i-1][j-1]+=c[i-jo-1][j-jo-1];//在r中jo+1，jo+2行以及jo+1，jo+2列位置累加k的单刚}}}else if(jo>=0)//如果io<0,即左节点为固定节点,jo>=0,右端为可动节点,则只在jo+1,jo+2对角分块位置累加{for(i=jo+1;i<=jo+2;i++){for(j=jo+1;j<=jo+2;j++){r[i-1][j-1]+=c[i-jo-1][j-jo-1];//在r中jo+1，jo+2行以及jo+1，jo+2列位置累加k的单刚}}}}}void choldlt()//总刚度阵的LDLT分解{int i,j,k;double m,t[M][M];for(i=0;i<nn;i++){l[i][i]=1.0;}d[0]=r[0][0];//d[0]=r[0][0]for(i=1;i<nn;i++){for(j=0;j<i;j++){m=0.0;for(k=0;k<j;k++){m+=t[i][k]*l[j][k];}t[i][j]=r[i][j]-m;l[i][j]=t[i][j]/d[j]; //计算l[i][j]m=0.0;for(k=0;k<i;k++){m+=t[i][k]*l[i][k];d[i]=r[i][i]-m; //计算d[i]}l[j][i]=l[i][j];}}}void trildlt()//回代求节点位移{int i,k;double m,y[M];y[0]=dp[0];for(i=1;i<nn;i++){m=0.0;for(k=0;k<i;k++){m+=l[i][k]*y[k];y[i]=dp[i]-m; //计算y[i] }}u[nn-1]=y[nn-1]/d[nn-1];for(i=nn-1;i>=0;i--){m=0.0;for(k=i+1;k<nn;k++){m+=l[k][i]*u[k];u[i]=y[i]/d[i]-m; //计算u[i]}}}void doh()//计算杆件的轴力{int i,k;for(k=1;k<=m;k++){iojo(k);//调用第k号杆件的左右端点的位移指示数for(i=0;i<2;i++)//计算每个节点2个自由度循环{if(io<0)//把右节点的2个位移存入v[0],v[1]{v[i]=u[jo+i];}else//把右节点的2个位移存入v[0],v[1]{v[i]=u[jo+i]-u[io+i];}}stif(k);//计算第k号杆件总体系单刚，存入[2]h[k-1]=0.0;//数组h[k-1]清零for(i=1;i<=2;i++)//对两个位移循环{h[k-1]=h[k-1]+t[i-1]*v[i-1]*rd;//轴力存入h[k-1]}}}void dowt() //考虑自重，且规定y轴竖直向上{int k,ko,i;double g; //g为重力printf("请输入杆件质量:\n");for(i=0;i<m;i++) //各杆件质量值输入{printf("mm[%d]=",i+1);scanf("%lf",&mm[i]);}for(k=1;k<=m;k++) //对桁架杆件循环{iojo(k); //调用函数ch(k);g=mm[k-1]*9.80665; //重力计算公式if(io>=0) //左节点为自由节点{dp[io+1]=dp[io+1]-(g/2); //左节点的y轴荷载减少二分之一重力dp[jo+1]=dp[jo+1]-(g/2); //右节点的y轴荷载减少二分之一重力}else if(jo>=0) //若右节点为自由节点，则仅有右节点做如下处理{ko=io+2*nc; //定义反力指示数ko 等于2*（固定节点号-1）dp[jo+1]=dp[jo+1]-(g/2);fdp[ko+1]=fdp[ko+1]-(g/2); //支座反力叠加重力的一半}}}void dofanli() //计算反力{int k,ko;for(k=1;k<=m;k++) //对杆件循环{iojo(k); //引用函数ch(k);ko=io+2*nc; //记录左节点if(io<0) //左节点为固定节点{fdp[ko]=fdp[ko]-h[k-1]*clxy[1]; //为第i个节点x轴向力加杆件轴力反力fdp[ko+1]=fdp[ko+1]-h[k-1]*clxy[2]; //为第i个节点y轴向力加杆件轴力反力}}}void verify() //强度校核函数{int k,i;double sigema,sigema0,sigema1; //定义应力，拉伸许用应力，压缩许用应力printf("请输入各杆截面积:\n"); //截面面积输入for(i=0;i<m;i++){printf("a[%d]=",i+1);scanf("%lf",&a[i]);}printf("请输入杆件拉伸许用应力:\n"); //杆件拉伸许用应力输入scanf("%lf",&sigema0);printf("请输入杆件压缩许用应力(输入正数):\n"); //杆件压缩许用应力输入scanf("%lf",&sigema1);for(k=1;k<=m;k++) //对杆件循环{sigema=h[k-1]/a[k-1]; //计算公式轴力与面积之商if(sigema>sigema0||sigema<-1*sigema1) //对应力，许用应力进行比较（注：压应力为负值，所以不小于压缩许用应力）{printf("第%d根杆件超过许用应力，为危险杆件，请增加横截面积或更换其他材料\n",k);}}}void assemble() //装配应力计算{int k,ko;double l; //定义杆件被拉长l printf("请输入存在装配应力的杆件号:\n");scanf("%d",&k);printf("请输入杆件装配时的拉长长度:\n");scanf("%lf",&l);iojo(k); //引用函数ch(k);h[k-1]=h[k-1]+ea[k-1]*l/clxy[0]; //储存装配杆件的应力值if(io>=0){dp[io]=dp[io]+ea[k-1]*l*clxy[1]/clxy[0]; //左节点x轴方向附加载荷增加P=△l*EA/l乘其方向余弦dp[io+1]=dp[io+1]+ea[k-1]*l*clxy[2]/clxy[0]; //y轴方向同上操作dp[jo]=dp[jo]-ea[k-1]*l*clxy[1]/clxy[0]; //注：右节点与左节点附加装配应力相反dp[jo+1]=dp[jo+1]-ea[k-1]*l*clxy[2]/clxy[0];}else{ko=io+2*nc;dp[jo]=dp[jo]-ea[k-1]*l*clxy[1]/clxy[0]; //注：右节点与左节点附加装配应力相反dp[jo+1]=dp[jo+1]-ea[k-1]*l*clxy[2]/clxy[0];printf("%f,%f\n",dp[jo],dp[jo+1]);fdp[ko]=fdp[ko]+ea[k-1]*l*clxy[1]/clxy[0]; //固定节点ihl反力叠加装配应力fdp[ko+1]=fdp[ko+1]+ea[k-1]*l*clxy[1]/clxy[0];}}void tem() //计算温度应力{int k,ko;double t0,arf,t1,t2,detal; //定义变量，温差，热膨胀系数，初始温度，最终温度，温变引起的长度变化printf("请输入初始温度\n"); //变量输入scanf("%lf",&t1);printf("请输入最终温度\n");scanf("%lf",&t2);printf("请输入杆件的热膨胀系数\n");scanf("%lf",&arf);t0=t2-t1;for(k=1;k<=m;k++){iojo(k); //引用函数ch(k);detal=-1*arf*t0*clxy[0]; //等效为装配应力杆件受压h[k-1]=h[k-1]+ea[k-1]*detal/clxy[0];if(io>=0){dp[io]=dp[io]+ea[k-1]*detal*clxy[1]/clxy[0];//左节点x轴方向附加载荷增加P=△l*EA/l乘其方向余弦dp[io+1]=dp[io+1]+ea[k-1]*detal*clxy[2]/clxy[0]; //y轴方向同上操作dp[jo]=dp[jo]-ea[k-1]*detal*clxy[1]/clxy[0]; //注：右节点与左节点附加温度应力相反dp[jo+1]=dp[jo+1]-ea[k-1]*detal*clxy[2]/clxy[0];}else{ko=io+2*nc;dp[jo]=dp[jo]-ea[k-1]*detal*clxy[1]/clxy[0]; //注：右节点与左节点附加温度应力相反dp[jo+1]=dp[jo+1]-ea[k-1]*detal*clxy[2]/clxy[0];fdp[ko]=fdp[ko]+ea[k-1]*detal*clxy[1]/clxy[0]; //固定节点ihl反力叠加温度应力fdp[ko+1]=fdp[ko+1]+ea[k-1]*detal*clxy[1]/clxy[0];}}}int main(){int i;printf(" *****************求解平面桁架节点位移与杆端力程序*****************\n\n");printf("\n ------------------输入数据时请用空格或回车符间隔------------------\n");printf("\n\n请输入桁架节点总数n,固定节点数nc,杆件数m:\n");scanf("%d,%d,%d",&n,&nc,&m);nn=2*(n-nc);printf("请输入节点坐标:\n");for(i=0;i<n;i++){printf("x[%d]=",i+1);scanf("%lf",&x[i]);printf("y[%d]=",i+1);scanf("%lf",&y[i]);}printf("输入杆件左右节点号:\n");for(i=0;i<m;i++){printf("ihl[%d]=",i+1);scanf("%d",&ihl[i]);printf("ihr[%d]=",i+1);scanf("%d",&ihr[i]);}printf("请输入杆件的EA值[ea]:\n");for(i=0;i<m;i++){printf("ea[%d]=",i+1);scanf("%lf",&ea[i]);}printf("请输入节点载荷[dp]:\n");for(i=0;i<nn;i++){printf("dp[%d]=",i+1);scanf("%lf",&dp[i]);}dor();choldlt();printf("是否需要考虑自重，需要请输入‘1’，不需要请输入‘0’。