c语言计算平面桁架内力计算程序.docx
《c语言计算平面桁架内力计算程序.docx》由会员分享,可在线阅读,更多相关《c语言计算平面桁架内力计算程序.docx(13页珍藏版)》请在冰豆网上搜索。
![c语言计算平面桁架内力计算程序.docx](https://file1.bdocx.com/fileroot1/2023-6/26/33dd1674-de94-46ff-8ed7-6821fccd96bb/33dd1674-de94-46ff-8ed7-6821fccd96bb1.gif)
c语言计算平面桁架内力计算程序
c语言计算平面桁架内力计算程序
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#defineM5
intn,nc,nn,m,e,f;//节点总数,固定节点数,自由度数,杆件数intio,jo;//单根杆对号指示数
intihl[M],ihr[M];//杆件左右节点号
doublea[M];//各杆截面积
doublemm[M];//杆件质量
doubleea[M];//杆件EA的值
doublex[M],y[M];//节点坐标
doubledp[M];//总体系下的节点载荷
doublet[2];//0,1分别为坐标转换矩阵的cos(),sin()
doublec[2][2];//总体系下的单刚
doubleclxy[3];//0,1,2分别为杆长,正弦,余弦
doubleh[M];//杆件轴力
doubler[M][M];//总刚度阵
doublerd;//桁架轴力杆局部系单刚
doubleu[M];//桁架节点位移
doublev[2];//存放节点位移差
doubled[M];//LDLT分解时的D矩阵的对角线元素
doublel[M][M];////LDLT分解时的D矩阵的对角线元素doublefdp[M];//总体系下支座反力
voidiojo(intk)//计算对号指示数io,jo
{
inti,j;
i=ihl[k-1];//k号杆左节点号进入i
j=ihr[k-1];//k号杆节点右号进入i
io=2*(i-nc-1);//uxi前未知位移的个数
jo=2*(j-nc-1);//uyi前未知位移的个数
}
voidch(intk)//计算杆长与方向余弦函数
{
inti,j;
i=ihl[k-1];//k号杆左节点进入i
j=ihr[k-1];//k号杆右节点进入j
clxy[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轴余弦
}
voidstif(intk)//计算k号杆件总体系下的单元刚度阵
{
inti,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;
}
}
}
voiddor()//总体系下的总刚度阵的组集
{
inti,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的单刚
}
}
}
elseif(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的单刚
}
}
}
}
}
voidcholdlt()//总刚度阵的LDLT分解
{
inti,j,k;
doublem,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];
}
}
}
voidtrildlt()//回代求节点位移
{
inti,k;
doublem,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]}
}
}
voiddoh()//计算杆件的轴力
{
inti,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]}}}
voiddowt()竖直向上{intk,ko,i;doubleg;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);二分之一重力dp[jo+1]=dp[jo+1]-(g/2);二分之一重力}elseif(jo>=0)//考虑自重,且规定y轴//g为重力//各杆件质量值输入//对桁架杆件循环//调用函数//重力计算公式//左节点为自由节点//左节点的y轴荷载减少//右节点的y轴荷载减少//若右节点为自由节点,
则仅有右节点做如下处理{ko=io+2*nc;//定义反力指示数ko等于2*(固定节点号-1)dp[jo+1]=dp[jo+1]-(g/2);fdp[ko+1]=fdp[ko+1]-(g/2);//支座反力叠加重力的一半}}}
voiddofanli(){intk,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];件轴力反力fdp[ko+1]=fdp[ko+1]-h[k-1]*clxy[2];件轴力反力}}}
voidverify(){intk,i;doublesigema,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++)//计算反力//对杆件循环//引用函数//记录左节点//左节点为固定节点//为第i个节点x轴向力加杆//为第i个节点y轴向力加杆//强度校核函数//定义应力,拉伸许用应//截面面积输入//杆件拉伸许用应力输入//杆件压缩许用应力输入//对杆件循环
{
sigema=h[k-1]/a[k-1];//计算公式轴力与面积之商if(sigema>sigema0||sigema<-1*sigema1)//对应力,许用应力进行比较(注:
压应力为负值,所以不小于压缩许用应力){printf("第%d根杆件超过许用应力,为危险杆件,请增加横截面积或更换其他材料\n",k);}}}
voidassemble(){intk,ko;doublel;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];加P=△l*EA/l乘其方向余弦dp[io+1]=dp[io+1]+ea[k-1]*l*clxy[2]/clxy[0];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];力fdp[ko+1]=fdp[ko+1]+ea[k-1]*l*clxy[1]/clxy[0];}}
//装配应力计算//定义杆件被拉长l//引用函数//储存装配杆件的应力值//左节点x轴方向附加载荷增//y轴方向同上操作//注:
右节点与左节点附加装配//注:
右节点与左节点附加装配//固定节点ihl反力叠加装配应
voidtem()//计算温度应力{
intk,ko;
doublet0,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];
}
}
}
intmain()
{
inti;
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’。
\n");scanf("%d",&e);
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("\n\n请输出各个节点位移列向量[U]:
\n");
for(i=0;i<nn;i++)
{
printf("u[%d]=%f\n",i+1,u[i]);
}
doh();
printf("输出杆件轴力:
\n");
for(i=0;i<m;i++)
{
printf("h[%d]=%8.4f\n",i+1,h[i]);
}
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.4f\n",i,fdp[i]);
}
}
return0;}