1、c语言计算平面桁架内力计算程序c语言计算平面桁架内力计算程序 #include<stdio.h> #include<math.h> #include<stdlib.h> #define M 5 int n,nc,nn,m,e,f;/节点总数，固定节点数，自由度数，杆件数 int io,jo;/单根杆对号指示数 int ihlM,ihrM;/杆件左右节点号 double aM;/各杆截面积 double mmM;/杆件质量 double eaM;/杆件EA的值 double xM,yM;/节点坐标 double dpM;/总体系下的节点载荷 double t2
2、;/0,1分别为坐标转换矩阵的cos(),sin() double c22;/总体系下的单刚 double clxy3;/0,1,2分别为杆长，正弦，余弦 double hM;/杆件轴力 double rMM;/总刚度阵 double rd;/桁架轴力杆局部系单刚 double uM;/桁架节点位移 double v2;/存放节点位移差 double dM;/LDLT分解时的D矩阵的对角线元素 double lMM;/LDLT分解时的D矩阵的对角线元素 double fdpM;/总体系下支座反力 void iojo(int k)/计算对号指示数 io,jo int i,j; i=ihlk-1;
3、/k号杆左节点号进入i j=ihrk-1;/k号杆节点右号进入i io=2*(i-nc-1);/uxi前未知位移的个数 jo=2*(j-nc-1);/uyi前未知位移的个数 void ch(int k)/计算杆长与方向余弦函数 int i,j; i=ihlk-1;/k号杆左节点进入i j=ihrk-1;/k号杆右节点进入j clxy1=xj-1-xi-1;/k号杆x坐标差 clxy2=yj-1-yi-1;/k号杆y坐标差 clxy0=sqrt(clxy1*clxy1+clxy2*clxy2);/k号杆长 clxy1=clxy1/clxy0;/k号杆件x轴余弦 clxy2=clxy2/clxy0
4、;/k号杆件y轴余弦 void stif(int k)/计算k号杆件总体系下的单元刚度阵 int i,j; ch(k);/调用ch(),计算k号杆件杆长与余弦 t0=clxy1; t1=clxy2; rd=eak-1/clxy0; for(i=0;i<2;i+) for(j=0;j<2;j+) cij=ti*tj*rd; void dor()/总体系下的总刚度阵的组集 int i,j,k; for(i=0;i<nn;i+) for(j=0;j<nn;j+) rij=0.0;/总刚度阵清0 for(k=1;k<=m;k+) iojo(k);/调用k的对号指示函数，从
5、而确定组集位置 stif(k);/第k号杆件的总体系下的单刚 if(io>=0) for(i=io+1;i<=io+2;i+) for(j=io+1;j<=io+2;j+) ri-1j-1+=ci-io-1j-io-1;/在r中io+1，io+2行以及io+1，io+2列位置累加k的单刚 for(j=jo+1;j<=jo+2;j+) ri-1j-1-=ci-io-1j-jo-1;/在r中io+1，io+2行以及jo+1，jo+2列位 置累加k的负单刚 rj-1i-1=ri-1j-1; for(i=jo+1;i<=jo+2;i+) for(j=jo+1;j<=
6、jo+2;j+) ri-1j-1+=ci-jo-1j-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+) ri-1j-1+=ci-jo-1j-jo-1;/在r中jo+1，jo+2行以及jo+1，jo+2列位置累加k的单刚 void choldlt()/总刚度阵的LDLT分解 int i,j,k; double m,
7、tMM; for(i=0;i<nn;i+) lii=1.0; d0=r00;/d0=r00 for(i=1;i<nn;i+) for(j=0;j<i;j+) m=0.0; for(k=0;k<j;k+) m+=tik*ljk; tij=rij-m; lij=tij/dj; /计算lij m=0.0; for(k=0;k<i;k+) m+=tik*lik; di=rii-m; /计算di lji=lij; void trildlt()/回代求节点位移 int i,k; double m,yM; y0=dp0; for(i=1;i<nn;i+) m=0.0; f
8、or(k=0;k<i;k+) m+=lik*yk; yi=dpi-m; /计算yi unn-1=ynn-1/dnn-1; for(i=nn-1;i>=0;i-) m=0.0; for(k=i+1;k<nn;k+) m+=lki*uk; ui=yi/di-m; /计算ui 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个位移存入v0,v1 vi=ujo+i; else/把右
9、节点的2个位移存入v0,v1 vi=ujo+i-uio+i; stif(k);/计算第k号杆件总体系单刚，存入2 hk-1=0.0;/数组hk-1清零 for(i=1;i<=2;i+)/对两个位移循环 hk-1=hk-1+ti-1*vi-1*rd;/轴力存入hk-1 void dowt() 竖直向上 int k,ko,i; double g; printf(请输入杆件质量:n); for(i=0;i<m;i+) printf(mm%d=,i+1); scanf(%lf,&mmi); for(k=1;k<=m;k+) iojo(k); ch(k); g=mmk-1*9.
10、80665; if(io>=0) dpio+1=dpio+1-(g/2); 二分之一重力 dpjo+1=dpjo+1-(g/2); 二分之一重力 else if(jo>=0) /考虑自重，且规定y轴 /g为重力 /各杆件质量值输入 /对桁架杆件循环 /调用函数 /重力计算公式 /左节点为自由节点 /左节点的y轴荷载减少/右节点的y轴荷载减少 /若右节点为自由节点， 则仅有右节点做如下处理 ko=io+2*nc; /定义反力指示数ko等于2*（固定节点号-1） dpjo+1=dpjo+1-(g/2); fdpko+1=fdpko+1-(g/2); /支座反力叠加重力的一半 void
11、dofanli() int k,ko; for(k=1;k<=m;k+) iojo(k); ch(k); ko=io+2*nc; if(io<0) fdpko=fdpko-hk-1*clxy1; 件轴力反力 fdpko+1=fdpko+1-hk-1*clxy2; 件轴力反力 void verify() int k,i; double sigema,sigema0,sigema1; 力，压缩许用应力 printf(请输入各杆截面积:n); for(i=0;i<m;i+) printf(a%d=,i+1); scanf(%lf,&ai); printf(请输入杆件拉伸许用
12、应力:n); scanf(%lf,&sigema0); printf(请输入杆件压缩许用应力(输入正数):n); scanf(%lf,&sigema1); for(k=1;k<=m;k+) /计算反力 /对杆件循环 /引用函数 /记录左节点 /左节点为固定节点 /为第i个节点x轴向力加杆/为第i个节点y轴向力加杆 /强度校核函数 /定义应力，拉伸许用应 /截面面积输入 /杆件拉伸许用应力输入 /杆件压缩许用应力输入 /对杆件循环 sigema=hk-1/ak-1; /计算公式轴力与面积之商 if(sigema>sigema0|sigema<-1*sigema1
13、) /对应力，许用应力进行比较（注：压应力为负值，所以不小于压缩许用应力） printf(第%d根杆件超过许用应力，为危险杆件，请增加横截面积或更换其他材料n,k); void assemble() int k,ko; double l; printf(请输入存在装配应力的杆件号:n); scanf(%d,&k); printf(请输入杆件装配时的拉长长度:n); scanf(%lf,&l); iojo(k); ch(k); hk-1=hk-1+eak-1*l/clxy0; if(io>=0) dpio=dpio+eak-1*l*clxy1/clxy0; 加P=l*EA/
14、l乘其方向余弦 dpio+1=dpio+1+eak-1*l*clxy2/clxy0; dpjo=dpjo-eak-1*l*clxy1/clxy0; 应力相反 dpjo+1=dpjo+1-eak-1*l*clxy2/clxy0; else ko=io+2*nc; dpjo=dpjo-eak-1*l*clxy1/clxy0; 应力相反 dpjo+1=dpjo+1-eak-1*l*clxy2/clxy0; printf(%f,%fn,dpjo,dpjo+1); fdpko=fdpko+eak-1*l*clxy1/clxy0; 力 fdpko+1=fdpko+1+eak-1*l*clxy1/clxy0
15、; /装配应力计算 /定义杆件被拉长l /引用函数 /储存装配杆件的应力值 /左节点x轴方向附加载荷增/y轴方向同上操作 /注：右节点与左节点附加装配/注：右节点与左节点附加装配/固定节点ihl反力叠加装配应 void tem() /计算温度应力 int k,ko; double t0,arf,t1,t2,detal; /定义变量，温差，热膨胀系数，初始温度，最终温度，温变引起的长度变化 printf(请输入初始温度n); /变量输入 scanf(%lf,&t1); printf(请输入最终温度n); scanf(%lf,&t2); printf(请输入杆件的热膨胀系数n);
16、scanf(%lf,&arf); t0=t2-t1; for(k=1;k<=m;k+) iojo(k); /引用函数 ch(k); detal=-1*arf*t0*clxy0; /等效为装配应力杆件受压 hk-1=hk-1+eak-1*detal/clxy0; if(io>=0) dpio=dpio+eak-1*detal*clxy1/clxy0; /左节点x轴方向附加载荷增加P=l*EA/l乘其方向余弦 dpio+1=dpio+1+eak-1*detal*clxy2/clxy0; /y轴方向同上操作 dpjo=dpjo-eak-1*detal*clxy1/clxy0; /
17、注：右节点与左节点附加温度应力相反 dpjo+1=dpjo+1-eak-1*detal*clxy2/clxy0; else ko=io+2*nc; dpjo=dpjo-eak-1*detal*clxy1/clxy0; /注：右节点与左节点附加温度应力相反 dpjo+1=dpjo+1-eak-1*detal*clxy2/clxy0; fdpko=fdpko+eak-1*detal*clxy1/clxy0; /固定节点ihl反力叠加温度应力 fdpko+1=fdpko+1+eak-1*detal*clxy1/clxy0; int main() int i; printf( *求解平面桁架节点位移与
18、杆端力程序 *nn); printf(n -输入数据时请用空格或回车符间隔-n); printf(nn请输入桁架节点总数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,&xi); printf(y%d=,i+1); scanf(%lf,&yi); printf(输入杆件左右节点号:n); for(i=0;i<m;i+) printf(ihl%d=,i
19、+1); scanf(%d,&ihli); printf(ihr%d=,i+1); scanf(%d,&ihri); printf(请输入杆件的EA值ea:n); for(i=0;i<m;i+) printf(ea%d=,i+1); scanf(%lf,&eai); printf(请输入节点载荷dp:n); for(i=0;i<nn;i+) printf(dp%d=,i+1); scanf(%lf,&dpi); dor(); choldlt(); printf(是否需要考虑自重，需要请输入1，不需要请输入0。n); scanf(%d,&e);
20、 if(e=1) dowt(); e=0; printf(是否需要考虑装配应力，需要请输入1，不需要请输入0。n); scanf(%d,&e); if(e=1) assemble(); e=0; printf(是否需要考虑温度应力，需要请输入1，不需要请输入0。n); scanf(%d,&e); if(e=1) tem(); e=0; /插入结束 trildlt(); printf(nn请输出各个节点位移列向量U:n); for(i=0;i<nn;i+) printf(u%d=%fn,i+1,ui); doh(); printf(输出杆件轴力:n); for(i=0;i<m;i+) printf(h%d=%8.4fn,i+1,hi); printf(是否需要考虑强度校核，需要请输入1，不需要请输入0。n); scanf(%d,&e); if(e=1) verify(); e=0; printf(是否需要计算支座反力，需要请输入1，不需要请输入0。n); scanf(%d,&f); if(f=1) dofanli(); printf(输出节点反力:n); for(i=0;i<=2*nc-1;i+) printf(fdp%d=%8.4fn,i,fdpi); return 0;
