?? ganjian.cpp
字號:
#include<stdio.h>
#include<math.h>
#define ne 3 //單元數
#define nj 4 //節點數
#define nz 6 //支撐數
#define npj 0 //節點載荷數
#define npf 1 //非節點載荷數
#define nj3 12 //節點位移總數
#define dd 6 //半帶寬
#define e0 2.1E8 //彈性模量
#define a0 0.008 //截面積
#define i0 1.22E-4 //單元慣性距
#define pi 3.141592654
int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/
double gc[ne+1]={0.0,1.0,2.0,1.0};
double gj[ne+1]={0.0,90.0,0.0,90.0};
double mj[ne+1]={0.0,a0,a0,a0};
double gx[ne+1]={0.0,i0,i0,i0};
int zc[nz+1]={0,1,2,3,10,11,12};
double pj[npj+1][3]={{0.0,0.0,0.0}};
double pf[npf+1][5]={{0,0,0,0,0},{0,-20,1.0,2.0,2.0}};
double kz[nj3+1][dd+1],p[nj3+1];
double pe[7],f[7],f0[7],t[7][7];
double ke[7][7],kd[7][7];
//**kz[][]—整體剛度矩陣
//**ke[][]—整體坐標下的單元剛度矩陣
//**kd[][]—局部坐標下的單位剛度矩陣
//**t[][]—坐標變換矩陣
//**這是函數聲明
void jdugd(int);
void zb(int);
void gdnl(int);
void dugd(int);
//**主程序開始
void main()
{
int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;
double cl,wy[7];
int im,in,jn;
//***********************************************
//<功能：形成矩陣P>
//***********************************************
if(npj>0)
{
for(i=1;i<=npj;i++)
{ //把節點載荷送入P
j=pj[i][2];
p[j]=pj[i][1];
}
}
if(npf>0)
{
for(i=1;i<=npf;i++)
{ //求固端反力F0
hz=i;
gdnl(hz);
e=(int)pf[hz][3];
zb(e); //求單元號碼
for(j=1;j<=6;j++) //求坐標變換矩陣T
{
pe[j]=0.0;
for(k=1;k<=6;k++) //求等效節點載荷
{
pe[j]=pe[j]-t[k][j]*f0[k];
}
}
al=jm[e][1];
bl=jm[e][2];
p[3*al-2]=p[3*al-2]+pe[1]; //將等效節點載荷送到P中
p[3*al-1]=p[3*al-1]+pe[2];
p[3*al]=p[3*al]+pe[3];
p[3*bl-2]=p[3*bl-2]+pe[4];
p[3*bl-1]=p[3*bl-1]+pe[5];
p[3*bl]=p[3*bl]+pe[6];
}
}
//*************************************************
//<功能：生成整體剛度矩陣kz[][]>
for(e=1;e<=ne;e++) //按單元循環
{
dugd(e); //求整體坐標系中的單元剛度矩陣ke
for(i=1;i<=2;i++) //對行碼循環
{
for(ii=1;ii<=3;ii++)
{
h=3*(i-1)+ii; //元素在ke中的行碼
dh=3*(jm[e][i]-1)+ii; //該元素在KZ中的行碼
for(j=1;j<=2;j++)
{
for(jj=1;jj<=3;jj++) //對列碼循環
{
l=3*(j-1)+jj; //元素在ke中的列碼
zl=3*(jm[e][j]-1)+jj; //該元素在KZ中的行碼
dl=zl-dh+1; //該元素在KZ*中的行碼
if(dl>0)
kz[dh][dl]=kz[dh][dl]+ke[h][l]; //剛度集成
}
}
}
}
}
//**引入邊界條件**
for(i=1;i<=nz;i++) //按支撐循環
{
z=zc[i]; //支撐對應的位移數
kz[z][l]=1.0; //第一列置1
for(j=2;j<=dd;j++)
{
kz[z][j]=0.0; //行置0
}
if((z!=1))
{
if(z>dd)
j0=dd;
else if(z<=dd)
j0=z; //列（45度斜線）置0
for(j=2;j<=j0;j++)
kz[z-j+1][j]=0.0;
}
p[z]=0.0; //P置0
}
for(k=1;k<=nj3-1;k++)
{
if(nj3>k+dd-1) //求最大行碼
im=k+dd-1;
else if(nj3<=k+dd-1)
im=nj3;
in=k+1;
for(i=in;i<=im;i++)
{
l=i-k+1;
cl=kz[k][l]/kz[k][1]; //修改KZ
jn=dd-l+1;
for(j=1;j<=jn;j++)
{
m=j+i-k;
kz[i][j]=kz[i][j]-cl*kz[k][m];
}
p[i]=p[i]-cl*p[k]; //修改P
}
}
p[nj3]=p[nj3]/kz[nj3][1]; //求最后一個位移分量
for(i=nj3-1;i>=1;i--)
{
if(dd>nj3-i+1)
j0=nj3-i+1;
else j0=dd; //求最大列碼j0
for(j=2;j<=j0;j++)
{
h=j+i-1;
p[i]=p[i]-kz[i][j]*p[h];
}
p[i]=p[i]/kz[i][1]; //求其他位移分量
}
printf("\n");
printf("_____________________________________________________________\n");
printf("NJ U V CETA \n"); //輸出位移
for(i=1;i<=nj;i++)
{
}
printf("_____________________________________________________________\n");
//*根據E的值輸出相應E單元的N,Q,M(A,B)的結果**
printf("E N Q M \n");
//*計算軸力N,剪力Q,彎矩M*
for(e=1;e<=ne;e++) //按單元循環
{
jdugd(e); //求局部單元剛度矩陣kd
zb(e); //求坐標變換矩陣T
for(i=1;i<=2;i++)
{
for(ii=1;ii<=3;ii++)
{
h=3*(i-1)+ii;
dh=3*(jm[e][i]-1)+ii; //給出整體坐標下單元節點位移
wy[h]=p[dh];
}
}
for(i=1;i<=6;i++)
{
f[i]=0.0;
for(j=1;j<=6;j++)
{
for(k=1;k<=6;k++) //求由節點位移引起的單元節點力
{
f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];
}
}
}
if(npf>0)
{
for(i=1;i<=npf;i++) //按非節點載荷數循環
if(pf[i][3]==e) //找到荷載所在的單元
{
hz=i;
gdnl(hz); //求固端反力
for(j=1;j<=6;j++) //將固端反力累加
{
f[j]=f[j]+f0[j];
}
}
}
printf("%-3d(A) %9.5f %9.5f %9.5f\n",e,f[1],f[2],f[3]); //輸出單元A（i）端內力
printf(" (B) %9.5f %9.5f %9.5f\n",f[4],f[5],f[6]); //輸出單元B（i）端內力
}
return;
}
//**主程序結束**
//************************************************************
//<功能：將非節點載荷下的桿端力計算出來存入f0[]>
//************************************************************
void gdnl(int hz)
{
int ind,e;
double g,c,l0,d;
g=pf[hz][1]; //載荷值
c=pf[hz][2]; //載荷位置
e=(int)pf[hz][3]; //作用單元
ind=(int)pf[hz][4]; //載荷類型
l0=gc[e]; //桿長
d=l0-c;
if(ind==1)
{
f0[1]=0.0;
f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2; //均布載荷的固端反力
f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12;
f0[4]=0.0;
f0[5]=-g*c-f0[2];
f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);
}
else
{
if(ind==2) //橫向集中力的固端反力
{
f0[1]=0.0;
f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);
f0[3]=-(g*c*d*d)/(l0*l0);
f0[4]=0.0;
f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);
f0[6]=(g*c*c*d)/(l0*l0);
}
else
{
f0[1]=-(g*d/l0); //縱向集中力的固端反力
f0[2]=0.0;
f0[3]=0.0;
f0[4]=-g*c/l0;
f0[5]=0.0;
f0[6]=0.0;
}
}
}
//************************************************************
//<功能：構成坐標變換矩陣>
//************************************************************
void zb(int e)
{
double ceta,co,si;
int i,j;
ceta=(gj[e]*pi)/180; //角度變弧度
co=cos(ceta);
si=sin(ceta);
t[1][1]=co; //計算T右上角元素
t[1][2]=si;
t[2][1]=-si;
t[2][2]=co;
t[3][3]=1.0;
for(i=1;i<=3;i++)
{
for(j=1;j<=3;j++) //計算T的左下角元素
{
t[i+3][j+3]=t[i][j];
}
}
}
//*****************************************************
//<計算局部坐標下單元剛度矩陣kd[][]>
//*****************************************************
void jdugd(int e)
{
double A0,l0,j0;
int i;
int j;
A0=mj[e]; //面積
l0=gc[e]; //桿長
j0=gx[e]; //慣性鉅
for(i=0;i<=6;i++)
for(j=0;j<=6;j++) //kd清0
kd[i][j]=0.0;
kd[1][1]=e0*A0/l0;
kd[2][2]=12*e0*j0/pow(l0,3);
kd[3][2]=6*e0*j0/pow(l0,2);
kd[3][3]=4*e0*j0/l0;
kd[4][1]=-kd[1][1];
kd[4][4]=kd[1][1];
kd[5][2]=-kd[2][2]; //計算kd左下角各元素
kd[5][3]=-kd[3][2];
kd[5][5]=kd[2][2];
kd[6][2]=kd[3][2];
kd[6][3]=2*e0*j0/l0;
kd[6][5]=-kd[3][2];
kd[6][6]=kd[3][3];
for(i=1;i<=6;i++)
for(j=1;j<=i;j++) //將kd左下角元素按對稱原則送到右下角
kd[j][i]=kd[i][j];
}
//**********************************************************
//<計算整體坐標下單元剛度矩陣ke[][]>
//**********************************************************
void dugd(int e)
{
int i,k,j,m;
jdugd(e); //計算局部單元剛度矩陣kd
zb(e); //計算坐標變換矩陣T
for(i=1;i<=6;i++)
{
for(j=1;j<=6;j++)
{
ke[i][j]=0.0;
for(k=1;k<=6;k++)
{
for(m=1;m<=6;m++)
{
ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j]; //計算剛度坐標內單元剛度矩陣ke
}
}
}
}
}
//**程序結束**
?? 快捷鍵說明
復制代碼
Ctrl + C
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -