📄 平面应力分析.cpp
字号:
#include<iostream.h>
#include<math.h>
#include<iomanip.h>
const int sunyongle=0;
class CST{
public:
double RY[60],RX[60],B[3],C[3],DEL;
double SK[12000],EK[21];
double D1,D2,D3,EO,PO,t;
double X[500],Y[500],U[1000];
int NRX[60],NRY[60],JU[50],JV[50],JE[700][3],JD[1000],JLL[6];
int NN,NE,NF,LK,KU,KV,KRX,KRY;
int K1,K2,K3,KK;
double BM[3][6],CM[3][6];
public:
//数据初始化
CST(){
t=1;
}
//数据输入程序
void dataInput1(){
//输入弹性模量EO,泊松比PO,厚度t
//输入结点数NN,单元总数NE
//输入结点坐标X[N],Y[N]
//输入结点编号JE[N][3]
//输入约束：
//输入X约束结点数KU
//输入X方向位移约束的结点号JU[KU]
//输入Y约束结点数KV
//输入Y方向位移约束的结点号JV[KV]
//输入载荷：
//输入X载荷结点数KRX
//输入X方向载荷结点编号NRX[KRX]
//输入X方向结点载荷大小NX[KRX]
//输入Y方向载荷结点数KRY
//输入Y方向载荷结点编号NRY[KRY]
//输入Y方向结点载荷大小NY[KRY]
}
void dataInput2(){
cout<<"请输入弹性模量EO,泊松比PO,厚度t:";
cin>>EO>>PO>>t;
//输入结点和单元数
cout<<"请输入结点数NN,单元总数NE:";
cin>>NN>>NE;
//输入结点坐标
cout<<"请输入"<<NN<<"个结点坐标X:";
for(int i=0;i<NN;i++)
cin>>X[i];
cout<<"请输入"<<NN<<"个结点坐标Y:";
for(int j=0;j<NN;j++)
cin>>Y[j];
//输入结点编号
cout<<"请输入"<<NN<<"×3个结点编号JE:"<<endl;
for(i=0;i<NE;i++)
{
cout<<"第"<<i+1<<"个单元结点编号:";
for(j=0;j<3;j++){
cin>>JE[i][j];
JE[i][j]--;
}
}
//输入约束
cout<<"请输入X约束结点数KU,Y约束结点数KV:";
cin>>KU>>KV;
cout<<"请输入X方向位移约束的结点号JU:";
for(i=0;i<KU;i++){
cin>>JU[i];
JU[i]--;
}
cout<<"请输入Y方向位移约束的结点号JV:";
for(i=0;i<KV;i++){
cin>>JV[i];
JV[i]--;
}
//输入载荷
cout<<"请输入X载荷结点数KRX,Y载荷结点数KRY:";
cin>>KRX>>KRY;
cout<<"请输入X方向载荷结点编号NRX:";
for(i=0;i<KRX;i++){
cin>>NRX[i];
NRX[i]--;
}
cout<<"请输入X方向结点载荷大小NX:";
for(i=0;i<KRX;i++)
cin>>RX[i];
cout<<"请输入Y方向载荷结点编号NRY:";
for(i=0;i<KRY;i++){
cin>>NRY[i];
NRY[i]--;
}
cout<<"请输入Y方向结点载荷大小NY:";
for(i=0;i<KRY;i++)
cin>>RY[i];
}
//数据测试1
void dataTest1(){
EO=6.93;PO=0.3;
NN=4;NE=2;t=2;
//输入结点坐标
X[0]=0;X[1]=0;X[2]=2;X[3]=2;
Y[0]=1;Y[1]=0;Y[2]=0;Y[3]=1;
//输入结点编号
JE[0][0]=0;JE[0][1]=1;JE[0][2]=3;
JE[1][0]=1;JE[1][1]=2;JE[1][2]=3;
//输入约束
KU=KV=2;
JU[0]=0;JU[1]=1;
JV[0]=0;JV[1]=1;
//输入载荷
KRX=2;KRY=0;
NRX[0]=2;NRX[1]=3;
RX[0]=1;RX[1]=1;
}
//数据测试2
void dataTest2(){
EO=1;PO=0;
NN=6;NE=4;//t=1;
//输入结点坐标
X[0]=0;X[1]=0;X[2]=1;X[3]=0;X[4]=1;X[5]=2;
Y[0]=2;Y[1]=1;Y[2]=1;Y[3]=0;Y[4]=0;Y[5]=0;
//输入结点编号
JE[0][0]=2;JE[0][1]=0;JE[0][2]=1;
JE[1][0]=4;JE[1][1]=1;JE[1][2]=3;
JE[2][0]=1;JE[2][1]=4;JE[2][2]=2;
JE[3][0]=5;JE[3][1]=2;JE[3][2]=4;
//输入约束
KU=KV=3;
JU[0]=0;JU[1]=1;JU[2]=3;
JV[0]=3;JV[1]=4;JV[2]=5;
//输入载荷
KRX=2;KRY=2;
NRX[0]=2;NRX[1]=5;
RX[0]=-1;RX[1]=-0.5;
NRY[0]=0;NRY[1]=2;
RY[0]=-1.5;RY[1]=-1;
}
//数据测试3
void dataTest3(){
EO=1;PO=0;
NN=3;NE=1;t=1;
//输入结点坐标
X[0]=3;X[1]=0;X[2]=0;
Y[0]=0;Y[1]=0;Y[2]=4;
//输入结点编号
JE[0][0]=0;JE[0][1]=2;JE[0][2]=1;
//输入约束
KU=2;KV=1;
JU[0]=1;JU[1]=2;
JV[0]=1;
//输入载荷
KRX=0;KRY=1;
NRY[0]=0;
RY[0]=-1;
}
//核对数据信息
bool dataCheck(){
int s=0;char key;bool command=false;
cout<<"请确认信息："<<endl;
cout<<"弹性模量EO＝"<<EO<<"\t"<<"泊松比PO＝"<<PO<<"\t"<<"厚度t＝"
<<t<<"\t"<<"结点数NN＝"<<NN<<"\t"<<"单元总数NE＝"<<NE<<endl;
//输出结点坐标
cout<<NN<<"个结点坐标:";
for(int i=0;i<NN;i++){
cout<<i+1<<"("<<X[i]<<","<<Y[i]<<")"<<"\t";
}
cout<<endl;
//输出结点编号
for(i=0;i<NE;i++)
{
cout<<"第"<<i+1<<"个单元结点编号:";
for(int j=0;j<3;j++){
s=JE[i][j];
s++;
cout<<s<<"\t";
}
cout<<endl;
}
//输出约束
cout<<"X约束结点数KU＝"<<KU<<":"<<"\t";
for(i=0;i<KU;i++){
s=JU[i];
s++;
cout<<s<<"\t";
}
cout<<endl;
cout<<"Y约束结点数KV＝"<<KV<<":"<<"\t";
for(i=0;i<KV;i++){
s=JV[i];
s++;
cout<<s<<"\t";
}
cout<<endl;
//输出载荷
//X:
cout<<"X载荷结点数KRX＝"<<KRX<<":"<<"\t";
for(i=0;i<KRX;i++){
s=NRX[i];
s++;
cout<<s<<"\t";
}
cout<<endl;
cout<<"X方向结点载荷大小NX"<<"\t";
for(i=0;i<KRX;i++)
cout<<RX[i]<<"\t";
cout<<endl;
//Y:
cout<<"Y载荷结点数KRY＝"<<KRY<<":"<<"\t";
for(i=0;i<KRY;i++){
s=NRY[i];
s++;
cout<<s<<"\t";
}
cout<<endl;
cout<<"Y方向结点载荷大小NY"<<"\t";
for(i=0;i<KRY;i++)
cout<<RY[i]<<"\t";
cout<<endl;
//请求用户确认
s=1;
while(s){
cout<<"是否要进行计算？Y/N "<<"\t";
cin>>key;
if(key=='y'||key=='Y'){
command=true;
s=0;
}
else {
if(key=='n'||key=='N'){
command=false;
s=0;
}
else
cout<<"请输入正确的选择!"<<endl;
}
}
return command;
}
//数据处理程序
void dataTackle(){
NF=2*NN;//结点位移总数
SKDD();
//平面弹性矩阵[D]
D1=EO/(1.0-PO*PO);
D2=D1*PO;
D3=D1*(1.0-PO)/2;
//总刚度矩阵初始化
for(int i=0;i<LK;i++)
SK[i]=0;
//储存载荷至矩阵[U]
for(i=0;i<NF;i++)
U[i]=0;
for(i=0;i<KRX;i++)
{
K1=2*NRX[i];
U[K1]=RX[i];
}
for(i=0;i<KRY;i++)
{
K1=2*NRY[i]+1;
U[K1]=RY[i];
}
for(i=0;i<NE;i++)
{
SHAPE(i);//计算第i个单元形状参数
FEK(i);//计算单元刚度矩阵值EK
SKKE();
}
KDisplay(); //输出总刚度矩阵K
//处理单元约束
for(i=0;i<KU;i++)
FIXED(2*JU[i]);
for(i=0;i<KV;i++)
FIXED(2*JV[i]+1);
FDisplay(); //输出处理约束后的载荷
SOLVE();//计算结点位移
UDisplay();//输出结点位移
//计算单元应力
for(i=0;i<NE;i++){
SHAPE(i);
STRESS();
}
}
//SKDD():计算总刚度矩阵主对角元一维存储序号
void SKDD(){
JD[0]=0;
JD[1]=2;
int JO,IO,KO,M1,M2,M3;
for(int N=1;N<NN;N++)
{
//计算半带宽KK
KK=0;
for(int i=0;i<NE;i++){
IO=JE[i][0];
JO=JE[i][1];
KO=JE[i][2];
if(IO==N||JO==N||KO==N){
M1=N-IO;
M2=N-JO;
M3=N-KO;
if(M1>KK) KK=M1;
if(M2>KK) KK=M2;
if(M3>KK) KK=M3;
}
}
KK=2*KK;
JD[2*N]=JD[2*N-1]+KK+1;
JD[2*N+1]=JD[2*N]+KK+2;
LK=JD[NF-1]+1;
}
cout<<"弹性矩阵SK的长度为:LK="<<LK<<endl;
}
//SHAPE(int):计算单元形状参数
void SHAPE(int N){
double XO[6],YO[6];
for(int I=0;I<3;I++){
K1=JE[N][I];
K2=I+3;
XO[I]=X[K1];
XO[K2]=XO[I];
YO[I]=Y[K1];
YO[K2]=YO[I];
K2=2*I;
K3=K2+1;
JLL[K2]=K1*2;
JLL[K3]=K1*2+1;
}
for(I=0;I<3;I++){
K1=I+1;
K2=I+2;
B[I]=YO[K1]-YO[K2];
C[I]=XO[K2]-XO[K1];
}
DEL=(B[1]*C[2]-B[2]*C[1])/2;
}
//FEK(int):计算单元刚度矩阵值EK
void FEK(int i){
double Z;
for(int I=0;I<3;I++){
for(int J=0;J<6;J++){
BM[I][J]=0;
CM[I][J]=0;
}
}
for(I=0;I<3;I++){
K2=I+I+1;
K1=K2-1;
BM[0][K1]=B[I];
BM[2][K2]=B[I];
BM[1][K2]=C[I];
BM[2][K1]=C[I];
CM[0][K1]=D1*B[I];
CM[0][K2]=D2*C[I];
CM[1][K1]=D2*B[I];
CM[1][K2]=D1*C[I];
CM[2][K1]=D3*C[I];
CM[2][K2]=D3*B[I];
}
for(I=0;I<6;I++){
K1=I*(I+1)/2;
for(int II=0;II<=I;II++){
Z=0;
for(int JJ=0;JJ<3;JJ++){
Z=Z+BM[JJ][I]*CM[JJ][II];
EK[K1+II]=Z/DEL/4;
}
EK[K1+II]*=t;
}
}
//输出单元刚度矩阵
EKDisplay(i);
}
//SKKE():单元刚度矩阵向总刚度矩阵传送
void SKKE(){
int K;
for(int I=0;I<6;I++){
K1=JLL[I];
KK=I*(I+1)/2;
for(int J=0;J<=I;J++){
K2=JLL[J];
if(K2<K1) K=JD[K1]-K1+K2;
else
K=JD[K2]-K2+K1;
SK[K]=SK[K]+EK[KK+J];
}
}
}
//处理约束条件
void FIXED(int K){
int L,M,IA,IB;
L=JD[K];
if(K>0){
M=JD[K-1];
IA=M+1;
IB=L-1;
for(int I=IA;I<=IB;I++)
SK[I]=0;
}
IA=K+1;
for(int I=IA;I<NF;I++){
if((JD[I]-JD[I-1])>=(I-K+1))
M=JD[I]-I+K;
SK[M]=0;
}
U[K]=0;
}
//SOLVE():解线性代数方程组
void SOLVE(){
int IG,MI,IP,JG,MJ,IJ,JI,JK,JJ,I,J,K;
for(I=0;I<NF;I++){
IG=JD[I]-I;
if(I==0)
MI=0;
else
MI=JD[I-1]-IG+1;
for(J=MI;J<=I;J++){
IP=IG+J;
JG=JD[J]-J;
if(J==0) MJ=0;
else
MJ=JD[J-1]-JG+1;
IJ=MI;
if(MJ>MI) IJ=MJ;
JI=J-1;
for(K=IJ;K<=JI;K++){
JK=JD[K];
SK[IP]=SK[IP]-SK[IG+K]*SK[JK]*SK[JG+K];
}
if(I!=J){
JJ=JD[J];
U[I]-=SK[IP]*U[J];
SK[IP]=SK[IP]/SK[JJ];
}
}
JI=JD[I];
U[I]=U[I]/SK[JI];
}
for(I=0;I<NF;I++){
J=NF-I+1;
IG=JD[J]-J;
if(J==0) MI=0;
else
MI=JD[J-1]-IG+1;
JI=J-1;
for(K=MI;K<=JI;K++){
U[K]=U[K]-SK[IG+K]*U[J];
}
}
}
//STRESS():计算单元形心处应力
void STRESS(){
int I,J;
double UE[6],SGM[3],S1,S2,S3,A1,A2,SM1,SM2,CETA;
for(I=0;I<6;I++){
K1=JLL[I];
UE[I]=U[K1];
}
for(I=0;I<3;I++){
K2=I+1;
K1=K2-1;
CM[0][K1]=D1*B[I];
CM[0][K2]=D2*C[I];
CM[1][K1]=D2*B[I];
CM[1][K2]=D1*C[I];
CM[2][K1]=D3*C[I];
CM[2][K2]=D3*B[I];
}
for(I=0;I<3;I++){
SGM[I]=0;
for(J=0;J<6;J++){
SGM[I]=SGM[I]+CM[I][J]*UE[J]/DEL/2;
}
}
S1=SGM[1];
S2=SGM[2];
S3=SGM[3];
A1=sqrt(pow((S1-S2)/2,2)+pow(S3,2));
A2=(S1+S2)/2;
SM1=A2+A1;
SM2=A2-A1;
if(fabs(S3)>1.E-4)
CETA=atan((SM1-S1)/S3)*57.29578;
else{
if(S1>S2)
CETA=0;
else CETA=90;
}
}
//KDisplay():输出总刚度矩阵K
void KDisplay(){
int num=1,L=0,M=0,i=0;
cout<<"总刚度矩阵K:"<<endl;
cout<<SK[0]<<endl;
for(i=1;i<NF;i++){
L=JD[i]-JD[i-1];
M=i-L+1;
for(int j=0;j<=i;j++){
if(j<M)
cout<<setprecision(5)<<setw(7)<<sunyongle<<" ";
else{
cout<<setprecision(5)<<setw(7)<<SK[num]<<" ";
num++;
}
}
cout<<endl;
}
}
//Uisplay():输出位移值
void UDisplay(){
cout<<"结点位移值为:"<<endl;
for(int i=0;i<NN;i++){
if(i!=0&&i%4==0)
cout<<endl;
cout<<i+1<<"("<<U[2*i]<<","<<U[2*i+1]<<")"<<"\t";
}
cout<<endl;
}
//Fisplay():输出结点载荷
void FDisplay(){
cout<<"结点载荷为:"<<endl;
for(int i=0;i<NN;i++){
if(i!=0&&i%4==0)
cout<<endl;
cout<<i+1<<"("<<U[2*i]<<","<<U[2*i+1]<<")"<<"\t";
}
cout<<endl;
}
private:
//NDisplay():输出单元应力转换矩阵
void NDisplay(int I){
cout<<"第"<<I+1<<"个单元应力转换矩阵[B]:"<<endl;
for(I=0;I<3;I++){
for(int J=0;J<6;J++){
cout<<BM[I][J]<<" ";
}
cout<<endl;
}
}
//EKDisplay():输出单元刚度矩阵
void EKDisplay(int I){
cout<<"第"<<I+1<<"个单元刚度矩阵:"<<endl;
for(I=0;I<6;I++){
K1=I*(I+1)/2;
for(int II=0;II<=I;II++)
cout<<setprecision(5)<<setw(7)<<EK[K1+II]<<" ";
cout<<endl;
}
}
};
void main(){
CST test;
//test.dataInput1();
//test.dataInput2();
//test.dataTest1();
test.dataTest2();
if(test.dataCheck())
test.dataTackle();
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -