/*===========================*/ /* gccFEM3nodTS.c */ /*===========================*/ #include #include #include #include #include #define NDmax 5000 #define NEmax 3000 #define MAmax 10 #define NTmax NDmax*2 int main_part(char fnameR[50],char fnameW[50]); double CALA_PL3(int ne,int kakom[NEmax+1][4],double x[NDmax+1],double y[NDmax+1]); void BMAT_PL3(int ne,int kakom[NEmax+1][4],double x[NDmax+1],double y[NDmax+1],double bm[4][7]); void DMAT_PL(int NSTRES,double Eme,double poe,double d[4][4]); void SM_PL3(int ne,int kakom[NEmax+1][4],double x[NDmax+1],double y[NDmax+1], int NSTRES,double Em[NEmax+1],double po[NEmax+1],double t[NEmax+1],double sm[7][7]); void MV_PL3(int ne,int kakom[NEmax+1][4],double x[NDmax+1],double y[NDmax+1], double t[NEmax+1],double gamma[NEmax+1],double gkh[NEmax+1],double gkv[NEmax+1],double fevec[7]); void CALST_PL3TS(int ne,int kakom[NEmax+1][4],double x[NDmax+1],double y[NDmax+1],double deltaT[NDmax+1],double dist[NTmax+1], int NSTRES, double Em[NEmax+1],double po[NEmax+1],double t[NEmax+1],double alpha[NEmax+1],double ts[NEmax+1], double strsave[NEmax+1][4],double fevec[7],int noten[NEmax+1]); void PST_PL(double sigx,double sigy,double tauxy,double *ps1,double *ps2,double *ang); void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][4],double sm[7][7],double tsm[NTmax+1][NTmax+1]); void ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][4],double fevec[7],double ftvec[NTmax+1]); void CHBAND10(int n1,int m1,double array[NTmax+1][NTmax+1]); void CHBAND11(int n1,int m1,double array[NTmax+1][NTmax+1],double vector[NTmax+1]); int IBAND(int mm,double a[NTmax+1][NTmax+1],double fact0); void BAND(int mm,int ib,double a[NTmax+1][NTmax+1]); /*------------------------------------------------------------------------------*/ int main(int argc, char* argv[]) { char fnameR[50]; char fnameW[50]; strcpy(fnameR,argv[1]); strcpy(fnameW,argv[2]); main_part(fnameR,fnameW); return 0; } /*------------------------------------------------------------------------------*/ int main_part(char fnameR[50],char fnameW[50]) { /*2次元三角形要素応力解析 */ FILE *fin,*fout; char dat[256]=""; char sv[20]=""; int i,j,k,l,kk,ne,n,ia,ja,mn; char strcom[256]=""; /* 書出用コメント */ int NSTRES; /* 応力状態(0;平面歪み,1;平面応力) */ int NODT; /* 節点総数 */ int NELT; /* 要素総数 */ int MATEL; /* 材料種類数 */ int KOX; /* x方向変位既知節点数 */ int KOY; /* y方向変位既知節点数 */ int NF; /* 載荷節点数 */ int nt; /* 全自由度(総節点数×2[1節点自由度]) */ int nod=3; /* 1要素節点数 */ int nhen=2; /* 1節点自由度数 */ int mm; /* 既知節点変位処理後自由度(逆行列の寸法) */ int ib; /* バンド幅 */ double fact0=1.0E-30; /* 0判定 */ double elarea1; /* 要素面積 */ double sigx,sigy,tauxy,ps1,ps2,ang; clock_t t_sta,t_end; double cal_time; time_t jikoku; struct tm *lt; char *p; /*--------------------------------------------------------------------------------------*/ /*配列寸法宣言*/ /*--------------------------------------------------------------------------------------*/ int kakom[NEmax+1][4]; /* 要素構成節点番号 */ double Em[NEmax+1]; /* 要素弾性係数 */ double po[NEmax+1]; /* 要素ポアソン比 */ double t[NEmax+1]; /* 要素板厚 */ double gamma[NEmax+1]; /* 要素単位体積重量[質量ではない] */ double gkh[NEmax+1]; /* 要素の水平方向物体力係数[重力加速度に対する比率] */ double gkv[NEmax+1]; /* 要素の鉛直方向物体力係数[重力加速度に対する比率] */ double alpha[NEmax+1]; /* 要素線膨張係数 */ double ts[NEmax]; /* 材料引張強度 */ int matno[NEmax+1]; /* 材料種別No */ double wm[MAmax+1][9]; /* 材料物性作業用 */ double strsave[NEmax+1][4]; /* 応力記憶用 */ int noten[NEmax]; /* No-tension積分点識別 */ int ntsum; /* No-tension積分点数 */ int nnn; /* 収束計算回数制御 */ double dtest; /* 変位委増分総和 */ double ftest; /* 不平衡力総和 */ int icount; /* 規定以下変位増分自由度総和 */ int maxita=2000; /* 収束計算回数上限 */ double x[NDmax+1]; /* 節点x座標 */ double y[NDmax+1]; /* 節点y座標 */ double deltaT[NDmax+1]; /* 節点温度変化量[+:温度上昇] */ int nokx[NDmax+1]; /* x方向変位既知節点番号 */ int noky[NDmax+1]; /* y方向変位既知節点番号 */ int index[NTmax+1]; /* 既知節点変位処理のためのインデックス */ static double tsm[NTmax+1][NTmax+1]; /* 全体剛性マトリックス */ double f[NTmax+1]; /* 全体外力ベクトル */ double ftvec[NTmax+1]; /* 全体荷重ベクトル */ double fmass[NTmax+1]; /* 全体物体力ベクトル */ double dis[NTmax+1]; /* 全体節点変位増分 */ double rdis[NTmax+1]; /* 既知全体節点変位 */ double wdis[NTmax+1]; /* 全体節点変位作業用記憶 */ double dist[NTmax+1]; /* 全体節点変位 */ double reac[NTmax+1]; /* 作業用配列 */ double sm[7][7]; /* 要素剛性マトリックス */ double fevec[7]; /* 要素荷重ベクトル */ /*時間計測開始*/ t_sta=clock(); /******************************************************/ /*データ入力*/ /******************************************************/ fin=fopen(fnameR,"r"); /*--------------------------------------------------------------------------------------*/ /*書出用コメント入力*/ /*--------------------------------------------------------------------------------------*/ fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1); /*--------------------------------------------------------------------------------------*/ /*応力状態[0;平面ひずみ],節点数,要素数,材料種類数,X方向変位既知節点数,Y方向変位既知節点数,載荷点数*/ /*--------------------------------------------------------------------------------------*/ fgets(dat,sizeof dat,fin); NSTRES=atoi(strtok(dat,",")); NODT =atoi(strtok(NULL,",")); NELT =atoi(strtok(NULL,",")); MATEL =atoi(strtok(NULL,",")); KOX =atoi(strtok(NULL,",")); KOY =atoi(strtok(NULL,",")); NF =atoi(strtok(NULL,",")); nt=NODT*2; /* 全自由度 */ /*--------------------------------------------------------------------------------------*/ /*材料特性(t,Em,po,gamma,gkh,gkv,alpha)*/ /*--------------------------------------------------------------------------------------*/ for(i=1;i<=MATEL;i++){ fgets(dat,sizeof dat,fin); wm[i][1]=atof(strtok(dat,",")); wm[i][2]=atof(strtok(NULL,",")); wm[i][3]=atof(strtok(NULL,",")); wm[i][4]=atof(strtok(NULL,",")); wm[i][5]=atof(strtok(NULL,",")); wm[i][6]=atof(strtok(NULL,",")); wm[i][7]=atof(strtok(NULL,",")); wm[i][8]=atof(strtok(NULL,",")); } /*--------------------------------------------------------------------------------------*/ /*要素構成節点と材料種別(node-1,node-2,node-3,matno)*/ /*--------------------------------------------------------------------------------------*/ for(ne=1;ne<=NELT;ne++){ fgets(dat,sizeof dat,fin); kakom[ne][1]=atoi(strtok(dat,",")); kakom[ne][2]=atoi(strtok(NULL,",")); kakom[ne][3]=atoi(strtok(NULL,",")); matno[ne]=atoi(strtok(NULL,",")); for(i=1;i<=MATEL;i++){ if(i==matno[ne]){ t[ne] =wm[i][1]; Em[ne] =wm[i][2]; po[ne] =wm[i][3]; gamma[ne]=wm[i][4]; gkh[ne] =wm[i][5]; gkv[ne] =wm[i][6]; alpha[ne]=wm[i][7]; ts[ne] =wm[i][7]; } } } /*--------------------------------------------------------------------------------------*/ /*節点座標(x,y),温度変化量*/ /*--------------------------------------------------------------------------------------*/ for(i=1;i<=NODT;i++){ fgets(dat,sizeof dat,fin); x[i] =atof(strtok(dat,",")); y[i] =atof(strtok(NULL,",")); deltaT[i]=atof(strtok(NULL,",")); } /*--------------------------------------------------------------------------------------*/ /*変位既知節点番号,既知節点変位量*/ /*--------------------------------------------------------------------------------------*/ for(i=1;i<=nt;i++){rdis[i]=0.0;} if(00){ mm=mm+1;index[mm]=i; } } for(i=1;i<=mm;i++){ ia=index[i]; ftvec[i]=ftvec[ia]; for(j=1;j<=mm;j++){ ja=index[j]; tsm[i][j]=tsm[ia][ja]; } } /*-----------------------------------------------------*/ /*バンド幅計算とフルマトリックスのバンド化記憶 /*-----------------------------------------------------*/ ib=IBAND(mm,tsm,fact0); BAND(mm,ib,tsm); /*-----------------------------------------------------*/ /*対角項確認と異常終了通知*/ /*-----------------------------------------------------*/ for(i=1;i<=mm;i++){ if(fabs(tsm[i][1])0) *ang=45.0; if(tauxy<0) *ang=135.0; if(tauxy==0.0) *ang=0.0; }else{ *ang=180.0/M_PI*(0.5*atan(2.0*tauxy/(sigx-sigy))); if(sigx>sigy&&tauxy>=0) *ang=*ang+0.0; if(sigx>sigy&&tauxy<0) *ang=*ang+180.0; if(sigxm1){ mm=m1; }else{ mm=n1-i+1; } for(j=1;j<=mm;j++){ z=array[i][j]; if(j!=m1){ if(m1-j+11){ array[i][j]=z/array[i][1]; }else{ array[i][1]=sqrt(z); } }else{ array[i][j]=z/array[i][1]; } } } } /*--------------------------------------------------------------------------------*/ void CHBAND11(int n1,int m1,double array[NTmax+1][NTmax+1],double vector[NTmax+1]) { /*前進消去・後退代入*/ int i,j,mm; double z; vector[1]=vector[1]/array[1][1]; for(i=2;i<=n1;i++){ z=vector[i]; if(i>m1){ mm=m1; }else{ mm=i; } for(j=2;j<=mm;j++){ z=z-array[i-j+1][j]*vector[i-j+1]; } vector[i]=z/array[i][1]; } vector[n1]=vector[n1]/array[n1][1]; for(i=n1-1;i>=1;i--){ z=vector[i]; for(j=2;j<=m1;j++){ if(i+j-1>n1)break; z=z-array[i][j]*vector[i+j-1]; } vector[i]=z/array[i][1]; } } /*--------------------------------------------------------------------------------*/ int IBAND(int mm,double array[NTmax+1][NTmax+1],double fact0) { int ib,i,j; /*-----------------------------------------------------*/ /*バンド幅計算 /*-----------------------------------------------------*/ ib=1; for(i=1;i<=mm-1;i++){ for(j=mm;j>=i+1;j--){ if(fact0