运用c语言 求解有限元单元刚度矩阵

/*函数功能：计算单元刚度矩阵

函数入口参数：实型数组nodeCoo[8][3]，存储单元的8个结点的总体自由度编号 实型数组elementK[24][24]，存储单元刚度矩阵

函数返回值：无

*/

void Getnumber(double nodeCoo[8][3],double elementK[24][24])

{

double Gauss[2]={0.57735026919,-0.57735026919}; //高斯点

double WeiG[2]={1.0,1.0}; //权重

int local[3][8]={{1,1,-1,-1,1,1,-1,-1},{-1,1,1,-1,-1,1,1,-1},{-1,-1,-1,-1,1,1,1,1}}; //结点局部坐标

double diffN[3][8]; //插值函数关于局部坐标偏导数

double diffNInverJ[3][8]={0}; //插值函数关于整体坐标偏导数

double DJ;

double J[3][3];

double *MatFirstadd;

double InverJ[3][3]; //存储J的逆矩阵 double B[6][24]={0}; //存储B矩阵 double BInver[24][6]; //存储B的转置 double *MatrixInverFirstadd;

double BInverD[24][6]={0}; //存储B的转置与D的乘积

int i,j,ksi,eta,zeta,m,n;

double Gaussksi,Gausseta,Gausszeta;

double G=E/(2*(1+v)),lam=E*v/((1+v)*(1-2*v)); //已知常数E，v

double

D[6][6]={{lam+2*G,lam,lam,0,0,0},{lam,lam+2*G,lam,0,0,0},{lam,lam,lam+2*G,0,0,0},{0,0,0,G,0,0},{0,0,0,0,G,0},{0,0,0,0,0,G}}; //D矩阵

double BTDBJ[24][24]={0};

for(ksi=0;ksi<2;ksi++) //以ksi方向高斯点计算积分

{

Gaussksi=Gauss[ksi];

for(eta=0;eta<2;eta++) //以eta方向高斯点计算积分

{

Gausseta=Gauss[eta];

for(zeta=0;zeta<2;zeta++) //以zeta方向高斯点计算积分

{

Gausszeta=Gauss[zeta];

for(i=0;i<3;i++) //以坐标进行循环

{

if (i==0)

{

for(j=0;j<8;j++) //以节点进行循环

{

运用c语言 求解有限元单元刚度矩阵

diffN[0][j]=(1.0/8)*Gaussksi*(1+Gausseta*local[1][j])*(1+Gausszeta*local[2][j]); //求插值函数关于ksi的偏导数

}

}

else if(i==1)

{

for(j=0;j<8;j++) //以节点进行循环

{

diffN[1][j]=(1.0/8)*Gausseta*(1+Gaussksi*local[0][j])*(1+Gausszeta*local[2][j]); //求插值函数关于eta的偏导数

}

}

else

{

for(j=0;j<8;j++) //以节点进行循环

{

diffN[2][j]=(1.0/8)*Gausszeta*(1+Gaussksi*local[0][j])*(1+Gausseta*local[1][j]); //求插值函数关于zeta的偏导数

}

}

}

for(i=0;i<3;i++) //计算雅可比矩阵

{

for(j=0;j<3;j++)

{

J[i][j]=0;

for(m=0;m<8;m++)

{

J[i][j]=J[i][j]+diffN[i][m]*nodeCoo[m][j];

}

}

}

m=3;

DJ=Surplus(J[0],m,m); //计算行列式J的值

MatFirstadd=MatrixOpp(J[0],m,m); //返回矩阵J[3][3]逆矩阵的首地址 for(i=0;i<3;i++)

{

for(j=0;j<3;j++)

{

InverJ[i][j]=*(MatFirstadd+i*3+j); //利用返回的首地址求出逆

运用c语言 求解有限元单元刚度矩阵

矩阵

} } for(i=0;i<3;i++) { for(j=0;j<8;j++) { for(m=0;m<3;m++) {

diffNInverJ[i][j]=diffNInverJ[i][j]+InverJ[i][m]*diffN[m][j]; //计算插值函数关于整体坐标xyz的偏导数

}

}

}

for(j=0;j<8;j++) //计算B矩阵

{

B[0][j*3]=diffNInverJ[0][j];

B[1][1+j*3]=diffNInverJ[1][j];

B[2][2+3*j]=diffNInverJ[2][j];

B[3][j*3]=B[1][1+j*3];

B[3][1+j*3]=B[0][j*3];

B[4][1+j*3]=B[2][2+j*3];

B[4][2+j*3]=B[1][1+j*3];

B[5][j*3]=B[2][2+j*3];

B[5][2+j*3]=B[0][j*3];

}

m=6;

n=24;

MatrixInverFirstadd=MatrixInver(B[0],m,n);

for(i=0;i<24;i++)

{

for(j=0;j<6;j++)

{

求B矩阵的转置

BInver[i][j]=*(MatrixInverFirstadd+i*6+j); //利用返回的首地址 } } for(i=0;i<24;i++) { for(j=0;j<6;j++) { for(m=0;m<6;m++)

运用c语言 求解有限元单元刚度矩阵

{

BInverD[i][j]=BInverD[i][j]+BInver[i][m]*D[m][j]; //计算B的转置与D的乘积

}

}

}

for(i=0;i<24;i++)

{

for(j=0;j<24;j++)

{

for(m=0;m<6;m++)

{

BTDBJ[i][j]=BTDBJ[i][j]+BInverD[i][m]*B[m][j]*DJ; //计算B的转置,D与B，雅可比矩阵的乘积

}

}

}

for(i=0;i<24;i++)

{

for(j=0;j<24;j++)

{

elementK[i][j]=BTDBJ[i][j]+elementK[i][j];

}

}

}

}

}

}