c语言计算平面桁架内力计算程序.docx
《c语言计算平面桁架内力计算程序.docx》由会员分享，可在线阅读，更多相关《c语言计算平面桁架内力计算程序.docx（13页珍藏版）》请在冰点文库上搜索。
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;}
升级会员 