1.程序源码
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define PI 3.141592654
int Gauss_Solve(double *x, double **A, double *b, int n)
{
double mik;
double s;
int i, j, k;
//消元
for(k=0; k<n-1; k++)
{
if(!A[k][k])
return -1;
for(i=k+1; i<n; i++)
{
mik = A[i][k] / A[k][k];
for(j=k; j<n; j++)
{
A[i][j] = A[i][j] - mik * A[k][j];
}
b[i] = b[i] - mik * b[k];
}
}
//printf("消元后的矩阵A: \n");
//for(i=0; i<n; i++)
//printf("%.4lf\n", b[i]);
//回代
x[n-1] = b[n-1] / A[n-1][n-1];
for(k=n-2; k>=0; k--)
{
s = b[k];
for(j=k+1; j<n; j++)
{
s = s - A[k][j] * x[j];
}
x[k] = s / A[k][k];
}
return 0;
}
//计算刚度矩阵
int Bar2D2Node_Stiffness(double **k, double E, double A, double x1, double y1, double x2, double y2, double alpha)
{
double L = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
double x = alpha * PI / 180;
double C = cos(x);
double S = sin(x);
k[0][0] = E * A / L * C * C;
k[0][1] = E * A / L * C * S;
k[0][2] = - E * A / L * C * C;
k[0][3] = - E * A / L * C * S;
k[1][0] = E * A / L * C * S;
k[1][1] = E * A / L * S * S;
k[1][2] = - E * A / L * C * S;
k[1][3] = - E * A / L * S * S;
k[2][0] = - E * A / L * C * C;
k[2][1] = - E * A / L * C * S;
k[2][2] = E * A / L * C * C;
k[2][3] = E * A / L * C * S;
k[3][0] = - E * A / L * C * S;
k[3][1] = - E * A / L * S * S;
k[3][2] = E * A / L * C * S;
k[3][3] = E * A / L * S * S;
return 0;
}
//单元刚度矩阵的组装, 0-u1, 1-v1, 2-u2, 3-v2
int Bar2D2Node_Assembly(double **KK, double **k, int i, int j)
{
KK[2*i+0][2*i+0] += k[0][0];
KK[2*i+0][2*i+1] += k[0][1];
KK[2*i+1][2*i+0] += k[1][0];
KK[2*i+1][2*i+1] += k[1][1];
KK[2*i+0][2*j+0] += k[0][2];
KK[2*i+0][2*j+1] += k[0][3];
KK[2*i+1][2*j+0] += k[1][2];
KK[2*i+1][2*j+1] += k[1][3];
KK[2*j+0][2*i+0] += k[2][0];
KK[2*j+0][2*i+1] += k[2][1];
KK[2*j+1][2*i+0] += k[3][0];
KK[2*j+1][2*i+1] += k[3][1];
KK[2*j+0][2*j+0] += k[2][2];
KK[2*j+0][2*j+1] += k[2][3];
KK[2*j+1][2*j+0] += k[3][2];
KK[2*j+1][2*j+1] += k[3][3];
return 0;
}
//边界条件处理
int Bar2D2Node_Condition(double **SS, double *b, double **KK, int *u, double *r)
{
int i, j;
int p = 0, q = 0;
for(i=0; i<8; i++)
{
if(u[i] == 0)
{
b[p] = r[i];
q = 0;
for(j=0; j<8; j++)
{
if(u[j] == 0)
{
SS[p][q] = KK[i][j];
q++;
}
}
p++;
}
}
return q;
}
//合并位移
int Bar2D2Node_Combine(double *r, double *b, int *u)
{
int i, j = 0;
for(i=0; i<8; i++)
{
if(u[i] == 0)
{
r[i] = b[j];
j++;
}
}
return j;
}
//获取单元位移矢量
int Bar2D2Node_Offset(double *v, double *u, int i, int j, int k, int w)
{
v[0] = u[i];
v[1] = u[j];
v[2] = u[k];
v[3] = u[w];
return 0;
}
//计算应力
int Bar2D2Node_Stress(double *stress, double E, double x1, double y1, double x2, double y2, double alpha, double *u)
{
double L = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
double x = alpha * PI / 180;
double C = cos(x);
double S = sin(x);
*stress = E / L * (C*u[2]+S*u[3]-C*u[0]-S*u[1]);
return 0;
}
//计算力矢量
int Bar2D2Node_Force(double *forces, double **KK, double *u, int n)
{
int i = 0;
int j = 0;
for(i=0; i<n; i++)
{
forces[i] = 0;
for(j=0; j<n; j++)
{
forces[i] += (KK[i][j] * u[j]);
}
}
return 0;
}
//输出单元刚度矩阵
int Bar2D2Node_Printk(double **k)
{
int i, j;
printf("[");
for(i=0; i<4; i++)
{
for(j=0; j<4; j++)
{
if(j != 3)
printf("%.8lf ", k[i][j]);
else
printf("%.8lf", k[i][j]);
}
if(i != 3)
printf("\n");
}
printf("]\n");
return 0;
}
//输出整体刚度矩阵
int Bar2D2Node_PrintK(double **k, int n)
{
int i, j;
printf("[");
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
if(j != n-1)
printf("%.8lf ", k[i][j]);
else
printf("%.8lf", k[i][j]);
}
if(i != n-1)
printf("\n");
}
printf("]\n");
return 0;
}
int main()
{
double rK[8][8];
double rS[8][8];
double *K[8];
double *S[8];
double rk[4][4][4];
double *k[4][4];
double b[8];
double c[8];
double r[8];
double x[8];
double v[4];
double stress[4];
double E = 2.95e11;
double A = 0.0001;
double alpha1 = 0;
double alpha2 = 90;
double alpha3 = atan(0.75)*180/PI;
double x1 = 0;
double y1 = 0;
double x2 = 0.4;
double y2 = 0;
double x3 = 0.4;
double y3 = 0.3;
double x4 = 0;
double y4 = 0.3;
int u[8];
int n = 0;
int i, j;
memset(stress, 0x00, sizeof(stress));
memset(rK, 0x00, sizeof(rK));
memset(rS, 0x00, sizeof(rS));
memset(b, 0x00, sizeof(b));
memset(c, 0x00, sizeof(c));
memset(u, 0x00, sizeof(u));
memset(x, 0x00, sizeof(x));
memset(v, 0x00, sizeof(v));
for(i=0; i<8; i++)
{
K[i] = rK[i];
S[i] = rS[i];
}
for(i=0; i<4; i++)
{
for(j=0; j<4; j++)
{
k[i][j] = rk[i][j];
}
}
u[0] = 1;
u[1] = 1;
u[3] = 1;
u[6] = 1;
u[7] = 1;
r[2] = 20000;
r[4] = 0;
r[5] = -25000;
Bar2D2Node_Stiffness(k[0], E, A, x1, y1, x2, y2, alpha1);
//Bar2D2Node_Printk(k[0]);
Bar2D2Node_Stiffness(k[1], E, A, x2, y2, x3, y3, alpha2);
//Bar2D2Node_Printk(k[1]);
Bar2D2Node_Stiffness(k[2], E, A, x1, y1, x3, y3, alpha3);
//Bar2D2Node_Printk(k[2]);
Bar2D2Node_Stiffness(k[3], E, A, x4, y4, x3, y3, alpha1);
//Bar2D2Node_Printk(k[3]);
Bar2D2Node_Assembly(K, k[0], 0, 1);
Bar2D2Node_Assembly(K, k[1], 1, 2);
Bar2D2Node_Assembly(K, k[2], 0, 2);
Bar2D2Node_Assembly(K, k[3], 3, 2);
n = Bar2D2Node_Condition(S, b, K, u, r);
Gauss_Solve(x, S, b, n);
Bar2D2Node_Combine(c, x, u);
Bar2D2Node_Force(r, K, c, 8);
/*for(i=0; i<8; i++)
{
printf("%.8lf ", c[i]);
//printf("%.8lf ", r[i]);
}
printf("\n");*/
Bar2D2Node_Offset(v, c, 0, 1, 2, 3);
//printf("%.8lf, %.8lf, %.8lf, %.8lf\n", v[0], v[1], v[2], v[3]);
Bar2D2Node_Stress(stress, E, x1, y1, x2, y2, alpha1, v);
Bar2D2Node_Offset(v, c, 2, 3, 4, 5);
Bar2D2Node_Stress(stress+1, E, x2, y2, x3, y3, alpha2, v);
Bar2D2Node_Offset(v, c, 0, 1, 4, 5);
Bar2D2Node_Stress(stress+2, E, x1, y1, x3, y3, alpha3, v);
Bar2D2Node_Offset(v, c, 6, 7, 4, 5);
Bar2D2Node_Stress(stress+3, E, x4, y4, x3, y3, alpha1, v);
printf("[%.8lf, %.8lf, %.8lf, %.8lf]\n", stress[0], stress[1], stress[2], stress[3]);
return 0;
}
2.编译源码
$ gcc -o Bar2D2Node Bar2D2Node.c
3.运行程序
./Bar2D2Node
[200000000.04486638, -218749999.98461720, -52083333.35897125, 41666666.64231063]
网友评论