c語(yǔ)言計(jì)算平面桁架內(nèi)力計(jì)算程序

1、#include<stdio.h>#include<math.h>#include<stdlib.h>#define M 5int n,nc,nn,m,e,f;/節(jié)點(diǎn)總數(shù)，固定節(jié)點(diǎn)數(shù)，自由度數(shù)，桿件數(shù)int io,jo;/單根桿對(duì)號(hào)指示數(shù)int ihlM,ihrM;/桿件左右節(jié)點(diǎn)號(hào)double aM;/各桿截面積double mmM;/桿件質(zhì)量double eaM;/桿件EA的值double xM,yM;/節(jié)點(diǎn)坐標(biāo)double dpM;/總體系下的節(jié)點(diǎn)載荷double t2;/0,1分別為坐標(biāo)轉(zhuǎn)換矩陣的cos(),sin()double c22;/總體系下的

2、單剛double clxy3;/0,1,2分別為桿長(zhǎng)，正弦，余弦double hM;/桿件軸力double rMM;/總剛度陣double rd;/桁架軸力桿局部系單剛double uM;/桁架節(jié)點(diǎn)位移double v2;/存放節(jié)點(diǎn)位移差double dM;/LDLT分解時(shí)的D矩陣的對(duì)角線元素double lMM;/LDLT分解時(shí)的D矩陣的對(duì)角線元素double fdpM;/總體系下支座反力void iojo(int k)/計(jì)算對(duì)號(hào)指示數(shù) io,jo int i,j; i=ihlk-1;/k號(hào)桿左節(jié)點(diǎn)號(hào)進(jìn)入i j=ihrk-1;/k號(hào)桿節(jié)點(diǎn)右號(hào)進(jìn)入i io=2*(i-nc-1);/uxi前未知

3、位移的個(gè)數(shù) jo=2*(j-nc-1);/uyi前未知位移的個(gè)數(shù)void ch(int k)/計(jì)算桿長(zhǎng)與方向余弦函數(shù) int i,j; i=ihlk-1;/k號(hào)桿左節(jié)點(diǎn)進(jìn)入i j=ihrk-1;/k號(hào)桿右節(jié)點(diǎn)進(jìn)入j clxy1=xj-1-xi-1;/k號(hào)桿x坐標(biāo)差 clxy2=yj-1-yi-1;/k號(hào)桿y坐標(biāo)差 clxy0=sqrt(clxy1*clxy1+clxy2*clxy2);/k號(hào)桿長(zhǎng) clxy1=clxy1/clxy0;/k號(hào)桿件x軸余弦 clxy2=clxy2/clxy0;/k號(hào)桿件y軸余弦void stif(int k)/計(jì)算k號(hào)桿件總體系下的單元?jiǎng)偠汝?int i,j; ch

4、(k);/調(diào)用ch(),計(jì)算k號(hào)桿件桿長(zhǎng)與余弦 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);/調(diào)用k的對(duì)號(hào)指示函數(shù)，從而確定組集位置 stif(k);/第k號(hào)桿件的總體系下的單剛if(io>=0) for(i=io+1;i&l

5、t;=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的負(fù)單剛 rj-1i-1=ri-1j-1; 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列位

6、置累加k的單剛 else if(jo>=0)/如果io<0,即左節(jié)點(diǎn)為固定節(jié)點(diǎn),jo>=0,右端為可動(dòng)節(jié)點(diǎn),則只在jo+1,jo+2對(duì)角分塊位置累加 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,tMM; for(i=0;i<nn;i+) lii=1.0; d0=r00;/d0=r00 for(i=1;i&

7、lt;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; /計(jì)算lij m=0.0; for(k=0;k<i;k+) m+=tik*lik; di=rii-m; /計(jì)算di lji=lij; void trildlt()/回代求節(jié)點(diǎn)位移 int i,k; double m,yM; y0=dp0; for(i=1;i<nn;i+) m=0.0; for(k=0;k<i;k+) m+=lik*yk; yi=dpi-m; /計(jì)算yi unn-1=ynn-1/dn

8、n-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; /計(jì)算ui void doh()/計(jì)算桿件的軸力 int i,k; for(k=1;k<=m;k+) iojo(k);/調(diào)用第k號(hào)桿件的左右端點(diǎn)的位移指示數(shù) for(i=0;i<2;i+)/計(jì)算每個(gè)節(jié)點(diǎn)2個(gè)自由度循環(huán) if(io<0)/把右節(jié)點(diǎn)的2個(gè)位移存入v0,v1 vi=ujo+i; else/把右節(jié)點(diǎn)的2個(gè)位移存入v0,v1 vi=ujo+i-uio+i; stif(k);/計(jì)算第k號(hào)桿件總體系單剛，存入2 hk

9、-1=0.0;/數(shù)組hk-1清零 for(i=1;i<=2;i+)/對(duì)兩個(gè)位移循環(huán) hk-1=hk-1+ti-1*vi-1*rd;/軸力存入hk-1 void dowt() /考慮自重，且規(guī)定y軸豎直向上int k,ko,i;double g; /g為重力printf("請(qǐng)輸入桿件質(zhì)量:n");for(i=0;i<m;i+) /各桿件質(zhì)量值輸入printf("mm%d=",i+1);scanf("%lf",&mmi);for(k=1;k<=m;k+) /對(duì)桁架桿件循環(huán)iojo(k); /調(diào)用函數(shù)ch(k);g

10、=mmk-1*9.80665; /重力計(jì)算公式if(io>=0) /左節(jié)點(diǎn)為自由節(jié)點(diǎn)dpio+1=dpio+1-(g/2); /左節(jié)點(diǎn)的y軸荷載減少二分之一重力dpjo+1=dpjo+1-(g/2); /右節(jié)點(diǎn)的y軸荷載減少二分之一重力else if(jo>=0) /若右節(jié)點(diǎn)為自由節(jié)點(diǎn)，則僅有右節(jié)點(diǎn)做如下處理ko=io+2*nc; /定義反力指示數(shù)ko等于2*（固定節(jié)點(diǎn)號(hào)-1）dpjo+1=dpjo+1-(g/2);fdpko+1=fdpko+1-(g/2); /支座反力疊加重力的一半void dofanli() /計(jì)算反力int k,ko;for(k=1;k<=m;k+)

11、/對(duì)桿件循環(huán)iojo(k); /引用函數(shù)ch(k);ko=io+2*nc; /記錄左節(jié)點(diǎn)if(io<0) /左節(jié)點(diǎn)為固定節(jié)點(diǎn)fdpko=fdpko-hk-1*clxy1; /為第i個(gè)節(jié)點(diǎn)x軸向力加桿件軸力反力fdpko+1=fdpko+1-hk-1*clxy2; /為第i個(gè)節(jié)點(diǎn)y軸向力加桿件軸力反力void verify() /強(qiáng)度校核函數(shù)int k,i;double sigema,sigema0,sigema1; /定義應(yīng)力，拉伸許用應(yīng)力，壓縮許用應(yīng)力 printf("請(qǐng)輸入各桿截面積:n"); /截面面積輸入for(i=0;i<m;i+)printf(&qu

12、ot;a%d=",i+1);scanf("%lf",&ai);printf("請(qǐng)輸入桿件拉伸許用應(yīng)力:n"); /桿件拉伸許用應(yīng)力輸入scanf("%lf",&sigema0);printf("請(qǐng)輸入桿件壓縮許用應(yīng)力(輸入正數(shù)):n"); /桿件壓縮許用應(yīng)力輸入scanf("%lf",&sigema1);for(k=1;k<=m;k+) /對(duì)桿件循環(huán)sigema=hk-1/ak-1; /計(jì)算公式軸力與面積之商if(sigema>sigema0|sig

13、ema<-1*sigema1) /對(duì)應(yīng)力，許用應(yīng)力進(jìn)行比較（注：壓應(yīng)力為負(fù)值，所以不小于壓縮許用應(yīng)力）printf("第%d根桿件超過許用應(yīng)力，為危險(xiǎn)桿件，請(qǐng)?jiān)黾訖M截面積或更換其他材料n",k);void assemble() /裝配應(yīng)力計(jì)算int k,ko;double l; /定義桿件被拉長(zhǎng)lprintf("請(qǐng)輸入存在裝配應(yīng)力的桿件號(hào):n");scanf("%d",&k);printf("請(qǐng)輸入桿件裝配時(shí)的拉長(zhǎng)長(zhǎng)度:n");scanf("%lf",&l);iojo(k)

14、; /引用函數(shù)ch(k);hk-1=hk-1+eak-1*l/clxy0; /儲(chǔ)存裝配桿件的應(yīng)力值if(io>=0)dpio=dpio+eak-1*l*clxy1/clxy0; /左節(jié)點(diǎn)x軸方向附加載荷增加P=l*EA/l乘其方向余弦dpio+1=dpio+1+eak-1*l*clxy2/clxy0; /y軸方向同上操作dpjo=dpjo-eak-1*l*clxy1/clxy0; /注：右節(jié)點(diǎn)與左節(jié)點(diǎn)附加裝配應(yīng)力相反 dpjo+1=dpjo+1-eak-1*l*clxy2/clxy0;elseko=io+2*nc; dpjo=dpjo-eak-1*l*clxy1/clxy0; /注：右節(jié)

15、點(diǎn)與左節(jié)點(diǎn)附加裝配應(yīng)力相反 dpjo+1=dpjo+1-eak-1*l*clxy2/clxy0;printf("%f,%fn",dpjo,dpjo+1);fdpko=fdpko+eak-1*l*clxy1/clxy0; /固定節(jié)點(diǎn)ihl反力疊加裝配應(yīng)力fdpko+1=fdpko+1+eak-1*l*clxy1/clxy0;void tem() /計(jì)算溫度應(yīng)力int k,ko;double t0,arf,t1,t2,detal; /定義變量，溫差，熱膨脹系數(shù)，初始溫度，最終溫度，溫變引起的長(zhǎng)度變化printf("請(qǐng)輸入初始溫度n"); /變量輸入scanf

16、("%lf",&t1);printf("請(qǐng)輸入最終溫度n");scanf("%lf",&t2);printf("請(qǐng)輸入桿件的熱膨脹系數(shù)n");scanf("%lf",&arf);t0=t2-t1;for(k=1;k<=m;k+)iojo(k); /引用函數(shù)ch(k);detal=-1*arf*t0*clxy0; /等效為裝配應(yīng)力桿件受壓 hk-1=hk-1+eak-1*detal/clxy0;if(io>=0)dpio=dpio+eak-1*detal*clx

17、y1/clxy0; /左節(jié)點(diǎn)x軸方向附加載荷增加P=l*EA/l乘其方向余弦dpio+1=dpio+1+eak-1*detal*clxy2/clxy0; /y軸方向同上操作dpjo=dpjo-eak-1*detal*clxy1/clxy0; /注：右節(jié)點(diǎn)與左節(jié)點(diǎn)附加溫度應(yīng)力相反 dpjo+1=dpjo+1-eak-1*detal*clxy2/clxy0;elseko=io+2*nc; dpjo=dpjo-eak-1*detal*clxy1/clxy0; /注：右節(jié)點(diǎn)與左節(jié)點(diǎn)附加溫度應(yīng)力相反 dpjo+1=dpjo+1-eak-1*detal*clxy2/clxy0;fdpko=fdpko+ea

18、k-1*detal*clxy1/clxy0; /固定節(jié)點(diǎn)ihl反力疊加溫度應(yīng)力fdpko+1=fdpko+1+eak-1*detal*clxy1/clxy0;int main() int i; printf(" *求解平面桁架節(jié)點(diǎn)位移與桿端力程序*nn");printf("n -輸入數(shù)據(jù)時(shí)請(qǐng)用空格或回車符間隔-n"); printf("nn請(qǐng)輸入桁架節(jié)點(diǎn)總數(shù)n,固定節(jié)點(diǎn)數(shù)nc,桿件數(shù)m:n"); scanf("%d,%d,%d",&n,&nc,&m); nn=2*(n-nc); printf

19、("請(qǐng)輸入節(jié)點(diǎn)坐標(biāo):n"); for(i=0;i<n;i+) printf("x%d=",i+1); scanf("%lf",&xi); printf("y%d=",i+1); scanf("%lf",&yi); printf("輸入桿件左右節(jié)點(diǎn)號(hào):n"); for(i=0;i<m;i+) printf("ihl%d=",i+1); scanf("%d",&ihli); printf("ih

20、r%d=",i+1); scanf("%d",&ihri); printf("請(qǐng)輸入桿件的EA值ea:n"); for(i=0;i<m;i+) printf("ea%d=",i+1); scanf("%lf",&eai); printf("請(qǐng)輸入節(jié)點(diǎn)載荷dp:n"); for(i=0;i<nn;i+) printf("dp%d=",i+1); scanf("%lf",&dpi); dor(); choldlt();printf("是否需要考慮自重，需要請(qǐng)輸入1，不需要請(qǐng)輸入0。n");scanf("%d",&e);if(e=1)dowt();e=0;printf("