using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; namespace WindowsFormsApplication1 { public partial class Form1 : Form { public Form1() { InitializeComponent(); } private void Form1_Load(object sender, EventArgs e) { label1.Dock = DockStyle.Fill; label1.TextAlign = ContentAlignment.MiddleCenter; label3.Dock = DockStyle.Fill; label3.TextAlign = ContentAlignment.MiddleCenter; label5.Dock = DockStyle.Fill; label5.TextAlign = ContentAlignment.MiddleCenter; label7.Dock = DockStyle.Fill; label7.TextAlign = ContentAlignment.MiddleCenter; label9.Dock = DockStyle.Fill; label9.TextAlign = ContentAlignment.MiddleCenter; label2.Dock = DockStyle.Fill; label2.TextAlign = ContentAlignment.MiddleLeft; label4.Dock = DockStyle.Fill; label4.TextAlign = ContentAlignment.MiddleLeft; label6.Dock = DockStyle.Fill; label6.TextAlign = ContentAlignment.MiddleLeft; label8.Dock = DockStyle.Fill; label8.TextAlign = ContentAlignment.MiddleLeft; label10.Dock = DockStyle.Fill; label10.TextAlign = ContentAlignment.MiddleLeft; label1.Text = "入力ファイル"; label3.Text = "出力ファイル"; label5.Text = "開始時刻"; label7.Text = "完了時刻"; label9.Text = "計算時間"; label2.Text = ""; label4.Text = ""; label6.Text = ""; label8.Text = ""; label10.Text = ""; } private void toolStripButton1_Click(object sender, EventArgs e) { main_part(); } private void main_part() { //四角形要素軸対称応力解析 int i, j, k, kk, ne, n, ia, ja; System.IO.StreamReader sr; System.IO.StreamWriter sw; string dat; string[] sbuf; char delim = ','; string fnameR = ""; string fnameW = ""; string strcom; //書出用コメント int NODT; //節点総数 int NELT; //要素総数 int MATEL; //材料種類数 int KOZ; //z方向変位既知節点数 int KOR; //r方向変位既知節点数 int NF; //載荷節点数 int nt; //全自由度(総節点数×2[1節点自由度]) int mm; //既知節点変位処理後自由度(逆行列の寸法) int ib; //バンド幅 int nod = 4; int nhen = 2; double fact0 = 1.0E-30; //0判定 double elarea1; //要素半面積 double elarea2; //要素半面積 double[,] sm = new double[9, 9]; //要素剛性マトリックス double[] fevec = new double[9]; //要素荷重ベクトル int maxita = 5000; int nnn; int ntsum; double dtest, ftest; int icount; int IPR; DateTime sttime; //計算開始時刻 DateTime entime; //計算終了時刻 TimeSpan tspan; //計算時間 //----------------------------------------------------- //入出力データファイル指定 //----------------------------------------------------- openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory(); if (openFileDialog1.ShowDialog() == DialogResult.OK) fnameR = openFileDialog1.FileName; if (saveFileDialog1.ShowDialog() == DialogResult.OK) fnameW = saveFileDialog1.FileName; label2.Text = fnameR; label4.Text = fnameW; //----------------------------------------------------- //時間計測開始 //----------------------------------------------------- sttime = DateTime.Now; label6.Text = sttime.ToString(); //***************************************************** //データ入力 //***************************************************** sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default); //-------------------------------------------------------------------------------------- //書出用コメント入力 //-------------------------------------------------------------------------------------- dat = sr.ReadLine(); sbuf = dat.Split(delim); strcom = sbuf[0]; //------------------------------------------------------------------------------------------ //応力状態[0;平面ひずみ],節点数,要素数,材料種類数,X方向変位既知節点数,Y方向変位既知節点数,載荷点数 //------------------------------------------------------------------------------------------ dat = sr.ReadLine(); sbuf = dat.Split(delim); NODT = int.Parse(sbuf[0]); NELT = int.Parse(sbuf[1]); MATEL = int.Parse(sbuf[2]); KOZ = int.Parse(sbuf[3]); KOR = int.Parse(sbuf[4]); NF = int.Parse(sbuf[5]); IPR = int.Parse(sbuf[6]); //------------------------------------------------------------------------------------------ //配列寸法宣言 //------------------------------------------------------------------------------------------ int[,] kakom = new int[NELT + 1, 5]; //要素構成節点番号 double[] Em = new double[NELT + 1]; //要素弾性係数 double[] po = new double[NELT + 1]; //要素ポアソン比 double[] alpha = new double[NELT + 1]; //要素線膨張係数 double[] ts = new double[NELT + 1]; //要素引張強さ int[] matno = new int[NELT + 1]; //材料種別No int[,] noten = new int[NELT + 1, 5]; //No-tension判定記憶 double[,] wm = new double[MATEL + 1, 5]; //材料物性作業用 double[, ,] strsave = new double[NELT + 1, 5, 5]; //応力記憶用 double[] z = new double[NODT + 1]; //節点z座標 double[] r = new double[NODT + 1]; //節点r座標 double[] deltaT = new double[NODT + 1]; //節点温度変化量[+:温度上昇] int[] nokz = new int[KOZ + 1]; //x方向変位既知節点番号 int[] nokr = new int[KOR + 1]; //y方向変位既知節点番号 nt = NODT * 2; //全自由度 int[] index = new int[nt + 1]; //既知節点変位処理のためのインデックス double[,] tsm = new double[nt + 1, nt + 1]; //全体剛性マトリックス double[] f = new double[nt + 1]; //全体外力ベクトル double[] ftvec = new double[nt + 1]; //全体荷重ベクトル double[] dis = new double[nt + 1]; //全体節点変位増分 double[] rdis = new double[nt + 1]; //既知全体節点変位 double[] wdis = new double[nt + 1]; //全体節点変位作業用記憶 double[] dist = new double[nt + 1]; //全体節点変位 double[] reac = new double[nt + 1]; //作業用配列 //------------------------------------------------------------------------------------------ //材料特性(t,Em,po,gamma,gkh,gkv,alpha,ts) //------------------------------------------------------------------------------------------ for (i = 1; i <= MATEL; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); for (j = 1; j <= 4; j++) { wm[i, j] = double.Parse(sbuf[j - 1]); } } //------------------------------------------------------------------------------------------ //要素構成節点と材料種別(node-1,node-2,node-3,node-4,matno) //------------------------------------------------------------------------------------------ for (ne = 1; ne <= NELT; ne++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); kakom[ne, 1] = int.Parse(sbuf[0]); kakom[ne, 2] = int.Parse(sbuf[1]); kakom[ne, 3] = int.Parse(sbuf[2]); kakom[ne, 4] = int.Parse(sbuf[3]); matno[ne] = int.Parse(sbuf[4]); for (i = 1; i <= MATEL; i++) { if (i == matno[ne]) { Em[ne] = wm[i, 1]; po[ne] = wm[i, 2]; alpha[ne] = wm[i, 3]; ts[ne] = wm[i, 4]; } } } //------------------------------------------------------------------------------------------ //節点座標(x,y),温度変化量 //------------------------------------------------------------------------------------------ for (i = 1; i <= NODT; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); z[i] = double.Parse(sbuf[0]); r[i] = double.Parse(sbuf[1]); deltaT[i] = double.Parse(sbuf[2]); } //------------------------------------------------------------------------------------------ //変位既知節点番号,既知節点変位量 //------------------------------------------------------------------------------------------ for (i = 1; i <= nt; i++) { rdis[i] = 0.0; } if (0 < KOZ) { for (i = 1; i <= KOZ; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokz[i] = int.Parse(sbuf[0]); rdis[2 * nokz[i] - 1] = double.Parse(sbuf[1]); } } if (0 < KOR) { for (i = 1; i <= KOR; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokr[i] = int.Parse(sbuf[0]); rdis[2 * nokr[i]] = double.Parse(sbuf[1]); } } //------------------------------------------------------------------------------------------ //節点番号,x方向荷重,y方向荷重 //------------------------------------------------------------------------------------------ for (i = 1; i <= nt; i++) { f[i] = 0.0; } if (0 < NF) { for (i = 1; i <= NF; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); n = int.Parse(sbuf[0]); f[2 * n - 1] = double.Parse(sbuf[1]); f[2 * n] = double.Parse(sbuf[2]); } } sr.Close(); //***************************************************** //入力確認出力 //***************************************************** INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOZ, KOR, NF, z, r, deltaT, f,rdis, nokz, nokr, kakom, Em, po, alpha, ts, matno,IPR); //要素面積確認と異常終了通知 for (ne = 1; ne <= NELT; ne++) { i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; elarea1 = 0.5 * ((z[k] - z[j]) * r[i] + (z[i] - z[k]) * r[j] + (z[j] - z[i]) * r[k]); if (elarea1 < fact0) IERROR1(fnameW); i = kakom[ne, 1]; j = kakom[ne, 3]; k = kakom[ne, 4]; elarea2 = 0.5 * ((z[k] - z[j]) * r[i] + (z[i] - z[k]) * r[j] + (z[j] - z[i]) * r[k]); if (elarea2 < fact0) IERROR1(fnameW); } //**************************************************** //メイン処理計算 //**************************************************** //----------------------------------------------------- //保存応力クリア //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { for (kk = 1; kk <= 4; kk++) { for (i = 1; i <= 4; i++) { strsave[ne, kk, i] = 0.0; } noten[ne, kk] = 0; } } //----------------------------------------------------- //行列・ベクトルクリア //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftvec[i] = 0.0; reac[i] = 0.0; wdis[i] = 0.0; dist[i] = 0.0; for (j = 1; j <= nt; j++) { tsm[i, j] = 0.0; } } //----------------------------------------------------- //要素剛性行列・荷重ベクトル組立 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { //剛性行列 SM_AX4(ne, kakom, z, r, Em, po, sm, noten); ASMAT(ne, nod, nhen, kakom, sm, tsm); } //----------------------------------------------------- //荷重ベクトルセット //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftvec[i] = f[i]; } //----------------------------------------------------- //変位既知節点自由度処理 //----------------------------------------------------- for (i = 1; i <= nt; i++) { index[i] = i; } if (0 < KOZ) { for (i = 1; i <= KOZ; i++) { n = 2 * nokz[i] - 1; index[n] = 0; } } if (0 < KOR) { for (i = 1; i <= KOR; i++) { n = 2 * nokr[i]; index[n] = 0; } } mm = 0; for (i = 1; i <= nt; i++) { if (index[i] > 0) { 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 (Math.Abs(tsm[i, 1]) < fact0) IERROR2(fnameW); } //----------------------------------------------------- //連立一次方程式(コレスキー法) //----------------------------------------------------- CHBAND10(mm, ib, tsm); //三角行列 //**************************************************** //No-tension収束計算 //**************************************************** nnn = 0; do { nnn = nnn + 1; //----------------------------------------------------- //連立一次方程式(コレスキー法) //----------------------------------------------------- for (i = 1; i <= mm; i++) { reac[i] = ftvec[i]; } CHBAND11(mm, ib, tsm, reac); //前進消去・後退代入 //----------------------------------------------------- //変位増分計算 //----------------------------------------------------- for (i = 1; i <= nt; i++) { dis[i] = 0.0; } for (i = 1; i <= mm; i++) { ia = index[i]; dis[ia] = reac[i]; } for (i = 1; i <= nt; i++) { reac[i] = 0.0; wdis[i] = wdis[i] + dis[i]; dist[i] = wdis[i] + rdis[i]; } //----------------------------------------------------- //応力計算・内力ベクトル組立 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { CALST_AX4TS(ne, kakom, z, r, deltaT, dist, Em, po, alpha, ts, strsave, fevec, noten); ASVEC(ne, nod, nhen, kakom, fevec, reac); //内力ベクトル組立 } //----------------------------------------------------- //不平衡力ベクトル算定 //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftvec[i] = 0.0; } for (i = 1; i <= nt; i++) { ftvec[i] = f[i] - reac[i]; } //----------------------------------------------------- //不平衡力ベクトル拘束節点処理 //----------------------------------------------------- for (i = 1; i <= mm; i++) { ia = index[i]; ftvec[i] = ftvec[ia]; dis[i] = dis[ia]; dist[i] = dist[ia]; } //----------------------------------------------------- //収束判定 //----------------------------------------------------- //No-tension積分点の数 ntsum = 0; for (ne = 1; ne <= NELT; ne++) { for (kk = 1; kk <= 4; kk++) { ntsum = ntsum + Math.Abs(noten[ne, kk]); } } //変位増分絶対値・不平衡力絶対値総和 icount = 0; dtest = 0.0; ftest = 0.0; for (i = 1; i <= mm; i++) { dtest = dtest + Math.Abs(dis[i]); //変位増分絶対値総和 ftest = ftest + Math.Abs(ftvec[i]); //不平衡力絶対値総和 if (Math.Abs(dis[i]/dist[i]) < 1e-6) icount = icount + 1; } //**************************************************** if (2 <= nnn) { if (ntsum == 0) break; if (icount == mm) break; } //**************************************************** } while (nnn <= maxita); // 出力用全変位調整 for (i = 1; i <= nt; i++) { dis[i] = 0.0; } for (i = 1; i <= mm; i++) { ia = index[i]; dis[ia] = dist[i]; } for (i = 1; i <= nt; i++) { dist[i] = dis[i] + rdis[i]; } //----------------------------------------------------- //時間計測完了 //----------------------------------------------------- entime = DateTime.Now; tspan = entime - sttime; label8.Text = entime.ToString(); label10.Text = tspan.TotalSeconds.ToString("0.00sec"); //***************************************************** //結果出力 //***************************************************** OUTRESUL(fnameW, NODT, NELT, nnn, reac, f, z, r, dist, strsave, matno, noten,IPR); sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "NODT=" + NODT.ToString() + " nt=" + nt.ToString() + " mm=" + mm.ToString() + " ib=" + ib.ToString(); sw.WriteLine(dat); dat = "nnn=" + nnn.ToString() + " dtest=" + dtest.ToString("E") + " ftest=" + ftest.ToString("E"); sw.WriteLine(dat); dat = "Calculation time=" + tspan.TotalSeconds.ToString("0.00sec"); sw.WriteLine(dat); dat = entime.ToString(); sw.WriteLine(dat); sw.Close(); MessageBox.Show("計算完了", "通知"); } //------------------------------------------------------------------------------ private void INPDAT(string fnameW, string strcom, int NODT, int NELT, int MATEL, int KOZ, int KOR, int NF, double[] z, double[] r, double[] deltaT, double[] f,double[] rdis, int[] nokz, int[] nokr, int[,] kakom, double[] Em, double[] po, double[] alpha, double[] ts, int[] matno,int IPR) { //入力確認出力 System.IO.StreamWriter sw; string dat; int i, j, k, l, ne; sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default); dat = strcom; sw.WriteLine(dat); dat = "NODT,NELT,MATEL,KOZ,KOR,NF,IPR"; sw.WriteLine(dat); dat = NODT.ToString("0"); dat = dat + "," + NELT.ToString("0"); dat = dat + "," + MATEL.ToString("0"); dat = dat + "," + KOZ.ToString("0"); dat = dat + "," + KOR.ToString("0"); dat = dat + "," + NF.ToString("0"); dat = dat + "," + IPR.ToString("0"); sw.WriteLine(dat); dat = "*node characteristics"; sw.WriteLine(dat); dat = "node,z,r,fz,fr,fix-z,fix-r,rdis-z,rdis-r,deltaT"; sw.WriteLine(dat); for (i = 1; i <= NODT; i++) { k = 0; l = 0; if (0 < KOZ) { for (j = 1; j <= KOZ; j++) { if (i == nokz[j]) k = 1; } } if (0 < KOR) { for (j = 1; j <= KOR; j++) { if (i == nokr[j]) l = 1; } } dat = i.ToString("0"); dat = dat + "," + z[i].ToString("E"); dat = dat + "," + r[i].ToString("E"); dat = dat + "," + f[2 * i - 1].ToString("E"); dat = dat + "," + f[2 * i].ToString("E"); dat = dat + "," + k.ToString("0"); dat = dat + "," + l.ToString("0"); dat = dat + "," + rdis[2 * i - 1].ToString("E"); dat = dat + "," + rdis[2 * i].ToString("E"); dat = dat + "," + deltaT[i].ToString("E"); sw.WriteLine(dat); } dat = "*element characteristics"; sw.WriteLine(dat); dat = "element,node-1,node-2,node-3,node-4,E,po,alpha,ts,matno"; sw.WriteLine(dat); for (ne = 1; ne <= NELT; ne++) { dat = ne.ToString("0"); dat = dat + "," + kakom[ne, 1].ToString("0"); dat = dat + "," + kakom[ne, 2].ToString("0"); dat = dat + "," + kakom[ne, 3].ToString("0"); dat = dat + "," + kakom[ne, 4].ToString("0"); dat = dat + "," + Em[ne].ToString("E"); dat = dat + "," + po[ne].ToString("E"); dat = dat + "," + alpha[ne].ToString("E"); dat = dat + "," + ts[ne].ToString("E"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } sw.Close(); } //------------------------------------------------------------------------------ private void OUTRESUL(string fnameW, int NODT, int NELT, int nnn, double[] reac, double[] f, double[] z, double[] r, double[] dist, double[, ,] strsave, int[] matno, int[,] noten, int IPR) { //結果出力(Appendモード) int i, ne, kk; System.IO.StreamWriter sw; string dat; double sigz, sigr, sigt, tauzr, ps1, ps2, ang; int ntsum; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = nnn.ToString("0") + "step*displacement and forces"; sw.WriteLine(dat); dat = "node,coord-z,coord-r,dis-z,dis-r,reac-z,reac-r,ftvec-z,ftvec-r"; sw.WriteLine(dat); //節点変位 for (i = 1; i <= NODT; i++) { dat = i.ToString("0"); dat = dat + "," + z[i].ToString("E"); dat = dat + "," + r[i].ToString("E"); dat = dat + "," + dist[2 * i - 1].ToString("E"); dat = dat + "," + dist[2 * i].ToString("E"); dat = dat + "," + reac[2 * i - 1].ToString("E"); dat = dat + "," + reac[2 * i].ToString("E"); dat = dat + "," + (f[2 * i - 1] - reac[2 * i - 1]).ToString("E"); dat = dat + "," + (f[2 * i] - reac[2 * i]).ToString("E"); sw.WriteLine(dat); } //要素応力 dat = "*stresses"; sw.WriteLine(dat); dat = "element,kk,sig-z,sig-r,sig-t,tau-zr,ps1,ps2,ang,noten,matno"; sw.WriteLine(dat); switch(IPR){ case 0: for (ne = 1; ne <= NELT; ne++) { for (kk = 1; kk <= 4; kk++) { sigz = strsave[ne, kk, 1]; sigr = strsave[ne, kk, 2]; sigt = strsave[ne, kk, 3]; tauzr = strsave[ne, kk, 4]; PST_AX4(sigz, sigr, tauzr, out ps1, out ps2, out ang); dat = ne.ToString("0"); dat = dat + "," + kk.ToString("0"); dat = dat + "," + sigz.ToString("E"); dat = dat + "," + sigr.ToString("E"); dat = dat + "," + sigt.ToString("E"); dat = dat + "," + tauzr.ToString("E"); dat = dat + "," + ps1.ToString("E"); dat = dat + "," + ps2.ToString("E"); dat = dat + "," + ang.ToString("E"); dat = dat + "," + noten[ne, kk].ToString("0"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } } break; case 1: for (ne = 1; ne <= NELT; ne++) { ntsum = 0; sigz = 0.0; sigr = 0.0; sigt = 0.0; tauzr = 0.0; for (kk = 1; kk <= 4; kk++) { sigz = sigz+strsave[ne, kk, 1]; sigr = sigr+strsave[ne, kk, 2]; sigt = sigt+strsave[ne, kk, 3]; tauzr = tauzr+strsave[ne, kk, 4]; ntsum = ntsum + noten[ne, kk]; } sigz = 0.25 * sigz; sigr = 0.25 * sigr; sigt = 0.25 * sigt; tauzr = 0.25 * tauzr; PST_AX4(sigz, sigr, tauzr, out ps1, out ps2, out ang); kk = 0; dat = ne.ToString("0"); dat = dat + "," + kk.ToString("0"); dat = dat + "," + sigz.ToString("E"); dat = dat + "," + sigr.ToString("E"); dat = dat + "," + sigt.ToString("E"); dat = dat + "," + tauzr.ToString("E"); dat = dat + "," + ps1.ToString("E"); dat = dat + "," + ps2.ToString("E"); dat = dat + "," + ang.ToString("E"); dat = dat + "," + ntsum.ToString("0"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } break; } sw.Close(); } //------------------------------------------------------------------------------ void DMAT_AX4(int ne, double Eme, double pzr, double pot, double[,] d) { // 応力-ひずみ関係行列[De] int i, j; for (i = 1; i <= 4; i++) { for (j = 1; j <= 4; j++) { d[i, j] = 0.0; } } d[1, 1] = 1.0 - pzr; d[1, 2] = pzr; d[1, 3] = pot; d[1, 4] = 0.0; d[2, 1] = pzr; d[2, 2] = 1.0 - pzr; d[2, 3] = pot; d[2, 4] = 0.0; d[3, 1] = pot; d[3, 2] = pot; d[3, 3] = 1.0 - pot; d[3, 4] = 0.0; d[4, 1] = 0.0; d[4, 2] = 0.0; d[4, 3] = 0.0; d[4, 4] = 0.5 * (1.0 - 2.0 * pzr); for (i = 1; i <= 4; i++) { for (j = 1; j <= 4; j++) { d[i, j] = d[i, j] * Eme / (1.0 + pzr) / (1.0 - 2.0 * pzr); } } } //------------------------------------------------------------------------------ void BMAT_AX4(int ne, int[,] kakom, double[] z, double[] r, double[,] bm, out double detJ, double a, double b) { /* ひずみ-変位関係行列[B] */ int i, j, k, l; double N1, N2, N3, N4, rr; double J11, J12, J21, J22; double dn1a, dn2a, dn3a, dn4a; double dn1b, dn2b, dn3b, dn4b; i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; l = kakom[ne, 4]; N1 = 0.25 * (1.0 - a) * (1.0 - b); N2 = 0.25 * (1.0 + a) * (1.0 - b); N3 = 0.25 * (1.0 + a) * (1.0 + b); N4 = 0.25 * (1.0 - a) * (1.0 + b); rr = N1 * r[i] + N2 * r[j] + N3 * r[k] + N4 * r[l]; /* [dN/dadN/db] */ dn1a = -0.25 * (1.0 - b); dn2a = 0.25 * (1.0 - b); dn3a = 0.25 * (1.0 + b); dn4a = -0.25 * (1.0 + b); dn1b = -0.25 * (1.0 - a); dn2b = -0.25 * (1.0 + a); dn3b = 0.25 * (1.0 + a); dn4b = 0.25 * (1.0 - a); /* Jacobiマトリックスとdet[J] */ J11 = dn1a * z[i] + dn2a * z[j] + dn3a * z[k] + dn4a * z[l]; J12 = dn1a * r[i] + dn2a * r[j] + dn3a * r[k] + dn4a * r[l]; J21 = dn1b * z[i] + dn2b * z[j] + dn3b * z[k] + dn4b * z[l]; J22 = dn1b * r[i] + dn2b * r[j] + dn3b * r[k] + dn4b * r[l]; detJ = J11 * J22 - J12 * J21; /*[B]行列要素[dN/dxdN/dy] */ for (i = 1; i <= 4; i++) { for (j = 1; j <= 8; j++) { bm[i, j] = 0.0; } } bm[1, 1] = J22 * dn1a - J12 * dn1b; bm[1, 3] = J22 * dn2a - J12 * dn2b; bm[1, 5] = J22 * dn3a - J12 * dn3b; bm[1, 7] = J22 * dn4a - J12 * dn4b; bm[2, 2] = -J21 * dn1a + J11 * dn1b; bm[2, 4] = -J21 * dn2a + J11 * dn2b; bm[2, 6] = -J21 * dn3a + J11 * dn3b; bm[2, 8] = -J21 * dn4a + J11 * dn4b; bm[3, 2] = N1 / rr * (detJ); bm[3, 4] = N2 / rr * (detJ); bm[3, 6] = N3 / rr * (detJ); bm[3, 8] = N4 / rr * (detJ); bm[4, 1] = -J21 * dn1a + J11 * dn1b; bm[4, 2] = J22 * dn1a - J12 * dn1b; bm[4, 3] = -J21 * dn2a + J11 * dn2b; bm[4, 4] = J22 * dn2a - J12 * dn2b; bm[4, 5] = -J21 * dn3a + J11 * dn3b; bm[4, 6] = J22 * dn3a - J12 * dn3b; bm[4, 7] = -J21 * dn4a + J11 * dn4b; bm[4, 8] = J22 * dn4a - J12 * dn4b; for (i = 1; i <= 4; i++) { for (j = 1; j <= 8; j++) { bm[i, j] = bm[i, j] / (detJ); } } } //------------------------------------------------------------------------------ void IPOINT4(int kk, out double a, out double b) { // Gauss積分点 a = 0.0; b = 0.0; switch (kk) { case 1: a = -0.5773502692; b = -0.5773502692; break; case 2: a = 0.5773502692; b = -0.5773502692; break; case 3: a = 0.5773502692; b = 0.5773502692; break; case 4: a = -0.5773502692; b = 0.5773502692; break; } } //------------------------------------------------------------------------------ void SM_AX4(int ne, int[,] kakom, double[] z, double[] r, double[] Em, double[] po, double[,] sm, int[,] noten) { /* 要素剛性行列 */ int kk, i, j, k; int i1, i2, i3, i4; double[,] d = new double[5, 5]; double[,] bm = new double[5, 9]; double[,] str = new double[5, 9]; double detJ, a, b; double N1, N2, N3, N4, rr; /* 配列クリア */ for (i = 1; i <= 8; i++) { for (j = 1; j <= 8; j++) { sm[i, j] = 0.0; } } /* 剛性行列[sm]=[B]T[D,B]*r*detJ */ DMAT_AX4(ne, Em[ne], po[ne], po[ne], d); i1 = kakom[ne, 1]; i2 = kakom[ne, 2]; i3 = kakom[ne, 3]; i4 = kakom[ne, 4]; for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); BMAT_AX4(ne, kakom, z, r, bm, out detJ, a, b); for (i = 1; i <= 4; i++) { for (j = 1; j <= 8; j++) { str[i, j] = 0.0; for (k = 1; k <= 4; k++) { str[i, j] = str[i, j] + d[i, k] * bm[k, j]; } } } N1 = 0.25 * (1.0 - a) * (1.0 - b); N2 = 0.25 * (1.0 + a) * (1.0 - b); N3 = 0.25 * (1.0 + a) * (1.0 + b); N4 = 0.25 * (1.0 - a) * (1.0 + b); rr = N1 * r[i1] + N2 * r[i2] + N3 * r[i3] + N4 * r[i4]; for (i = 1; i <= 8; i++) { for (j = 1; j <= 8; j++) { for (k = 1; k <= 4; k++) { sm[i, j] = sm[i, j] + bm[k, i] * str[k, j] * rr * detJ; } } } } } //------------------------------------------------------------------------------ void CALST_AX4TS(int ne, int[,] kakom, double[] z, double[] r, double[] deltaT, double[] dist, double[] Em, double[] po, double[] alpha, double[] ts, double[, ,] strsave, double[] fevec, int[,] noten) { /* 応力および内力ベクトル計算 */ int kk, i, j, ia; int i1, i2, i3, i4; double detJ, a, b; double[,] d = new double[5, 5]; double[,] bm = new double[5, 9]; double[] epsilon = new double[5]; double[] eps0 = new double[5]; double[] stress = new double[5]; double[] wd = new double[9]; double N1, N2, N3, N4, rr; double tem, pzr = 0.0, pot = 0.0; double sigz, sigr, sigt, tauzr, ps1, ps2, ang, phi; i1 = kakom[ne, 1]; i2 = kakom[ne, 2]; i3 = kakom[ne, 3]; i4 = kakom[ne, 4]; /* 要素座標節点変位 */ for (i = 1; i <= 4; i++) { ia = kakom[ne, i]; wd[2 * i - 1] = dist[2 * ia - 1]; wd[2 * i] = dist[2 * ia]; } /* 要素応力*/ for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); BMAT_AX4(ne, kakom, z, r, bm, out detJ, a, b); /*ひずみ計算 {epsilon}=[B]{u} */ for (i = 1; i <= 4; i++) { epsilon[i] = 0.0; for (j = 1; j <= 8; j++) { epsilon[i] = epsilon[i] + bm[i, j] * wd[j]; } } /*積分点温度ひずみ*/ N1 = 0.25 * (1.0 - a) * (1.0 - b); N2 = 0.25 * (1.0 + a) * (1.0 - b); N3 = 0.25 * (1.0 + a) * (1.0 + b); N4 = 0.25 * (1.0 - a) * (1.0 + b); tem = N1 * deltaT[i1] + N2 * deltaT[i2] + N3 * deltaT[i3] + N4 * deltaT[i4]; eps0[1] = alpha[ne] * tem; eps0[2] = alpha[ne] * tem; eps0[3] = alpha[ne] * tem; eps0[4] = 0.0; switch (noten[ne, kk]) { case 0: pzr = po[ne]; pot = po[ne]; break; case 1: pzr = 0.0; pot = po[ne]; break; case 2: pzr = 0.0; pot = po[ne]; break; case 3: pzr = po[ne]; pot = 0.0; break; case 4: pzr = 0.0; pot = 0.0; break; case 5: pzr = 0.0; pot = 0.0; break; } DMAT_AX4(ne, Em[ne], pzr, pot, d); /*応力計算 {stress}=[D]{epsilon} */ for (i = 1; i <= 4; i++) { stress[i] = 0.0; for (j = 1; j <= 4; j++) { stress[i] = stress[i] + d[i, j] * (epsilon[j] - eps0[j]); } } /*現状応力保存*/ sigz = stress[1]; sigr = stress[2]; sigt = stress[3]; tauzr = stress[4]; PST_AX4(sigz, sigr, tauzr, out ps1, out ps2, out ang); /*No-tension判定*/ noten[ne, kk] = 0; if (ts[ne] < ps1 && ps2 < ts[ne] && sigt < ts[ne]) noten[ne, kk] = 1; if (ts[ne] < ps1 && ts[ne] < ps2 && sigt < ts[ne]) noten[ne, kk] = 2; if (ps1 < ts[ne] && ps2 < ts[ne] && ts[ne] < sigt) noten[ne, kk] = 3; if (ts[ne] < ps1 && ps2 < ts[ne] && ts[ne] < sigt) noten[ne, kk] = 4; if (ts[ne] < ps1 && ts[ne] < ps2 && ts[ne] < sigt) noten[ne, kk] = 5; /*引張成分除去*/ phi = -ang / 180.0 * Math.PI; if (noten[ne, kk] == 1 || noten[ne, kk] == 2 || noten[ne, kk] == 4 || noten[ne, kk] == 5) { sigz = sigz - 0.5 * ps1 * (1.0 + Math.Cos(2.0 * phi)); sigr = sigr - 0.5 * ps1 * (1.0 - Math.Cos(2.0 * phi)); tauzr = tauzr + 0.5 * ps1 * Math.Sin(2.0 * phi); } if (noten[ne, kk] == 2 || noten[ne, kk] == 5) { sigz = sigz - 0.5 * ps2 * (1.0 - Math.Cos(2.0 * phi)); sigr = sigr - 0.5 * ps2 * (1.0 + Math.Cos(2.0 * phi)); tauzr = tauzr - 0.5 * ps2 * Math.Sin(2.0 * phi); } if (noten[ne, kk] == 3 || noten[ne, kk] == 4 || noten[ne, kk] == 5) { sigt = 0.0; } /*作業用応力保存*/ strsave[ne, kk, 1] = sigz; strsave[ne, kk, 2] = sigr; strsave[ne, kk, 3] = sigt; strsave[ne, kk, 4] = tauzr; } /* 内力ベクトル */ for (i = 1; i <= 8; i++) { fevec[i] = 0.0; } for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); BMAT_AX4(ne, kakom, z, r, bm, out detJ, a, b); N1 = 0.25 * (1.0 - a) * (1.0 - b); N2 = 0.25 * (1.0 + a) * (1.0 - b); N3 = 0.25 * (1.0 + a) * (1.0 + b); N4 = 0.25 * (1.0 - a) * (1.0 + b); rr = N1 * r[i1] + N2 * r[i2] + N3 * r[i3] + N4 * r[i4]; for (i = 1; i <= 8; i++) { for (j = 1; j <= 4; j++) { fevec[i] = fevec[i] + bm[j, i] * strsave[ne, kk, j] * rr * detJ; } } } } //------------------------------------------------------------------------------ void PST_AX4(double sigz, double sigr, double tauzr, out double ps1, out double ps2, out double ang) { //主応力計算 ps1 = 0.5 * (sigz + sigr) + Math.Sqrt(0.25 * (sigz - sigr) * (sigz - sigr) + tauzr * tauzr); ps2 = 0.5 * (sigz + sigr) - Math.Sqrt(0.25 * (sigz - sigr) * (sigz - sigr) + tauzr * tauzr); ang = 0.0; if (sigz == sigr) { if (tauzr > 0) ang = 45.0; if (tauzr < 0) ang = 135.0; if (tauzr == 0.0) ang = 0.0; } else { ang = 180.0 / Math.PI * (0.5 * Math.Atan(2.0 * tauzr / (sigz - sigr))); if (sigz > sigr && tauzr >= 0) ang = ang + 0.0; if (sigz > sigr && tauzr < 0) ang = ang + 180.0; if (sigz < sigr) ang = ang + 90.0; } } //------------------------------------------------------------------------------ private void ASMAT(int ne, int nod, int nhen, int[,] kakom, double[,] sm, double[,] tsm) { //剛性行列組立[nod:1要素節点数,nhen:1節点の自由度] int i, j, k, l, ki, kj, ks, js, kik, kjl, ksk, jsl; for (i = 1; i <= nod; i++) { for (j = 1; j <= nod; j++) { ki = (kakom[ne, i] - 1) * nhen; kj = (kakom[ne, j] - 1) * nhen; ks = (i - 1) * nhen; js = (j - 1) * nhen; for (k = 1; k <= nhen; k++) { for (l = 1; l <= nhen; l++) { kik = ki + k; kjl = kj + l; ksk = ks + k; jsl = js + l; tsm[kik, kjl] = tsm[kik, kjl] + sm[ksk, jsl]; } } } } } //------------------------------------------------------------------------------ private void ASVEC(int ne, int nod, int nhen, int[,] kakom, double[] fevec, double[] ftvec) { //荷重ベクトル組立[nod:1要素節点数,nhen:1節点の自由度] int i, k, ki, ks, kik, ksk; for (i = 1; i <= nod; i++) { ki = (kakom[ne, i] - 1) * nhen; ks = (i - 1) * nhen; for (k = 1; k <= nhen; k++) { kik = ki + k; ksk = ks + k; ftvec[kik] = ftvec[kik] + fevec[ksk]; } } } //------------------------------------------------------------------------------ private void CHBAND10(int n1, int m1, double[,] array) { //三角行列 int i, j, k, mm, mn; double z; array[1, 1] = Math.Sqrt(array[1, 1]); for (j = 2; j <= m1; j++) { array[1, j] = array[1, j] / array[1, 1]; } for (i = 2; i <= n1; i++) { if (n1 - i + 1 > m1) { mm = m1; } else { mm = n1 - i + 1; } for (j = 1; j <= mm; j++) { z = array[i, j]; if (j != m1) { if (m1 - j + 1 < i) { mn = m1 - j + 1; } else { mn = i; } for (k = 2; k <= mn; k++) { z = z - array[i - k + 1, k] * array[i - k + 1, j + k - 1]; } if (j > 1) { array[i, j] = z / array[i, 1]; } else { array[i, 1] = Math.Sqrt(z); } } else { array[i, j] = z / array[i, 1]; } } } } //------------------------------------------------------------------------------ private void CHBAND11(int n1, int m1, double[,] array, double[] vector) { //前進消去・後退代入 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]; } } //------------------------------------------------------------------------------ private int IBAND(int mm, double[,] array, double fact0) { int ib, i, j; //----------------------------------------------------- //バンド幅計算 //----------------------------------------------------- ib = 1; for (i = 1; i <= mm - 1; i++) { for (j = mm; j >= i + 1; j--) { if (fact0 <= Math.Abs(array[i, j])) break; } if (ib <= j - i + 1) ib = j - i + 1; } if (mm < ib) ib = mm; return ib; } //------------------------------------------------------------------------------ private void BAND(int mm, int ib, double[,] array) { int i, j, k, mn; double[] work = new double[mm + 1]; //----------------------------------------------------- //フルマトリックスのバンド化記憶 //----------------------------------------------------- for (i = 2; i <= mm; i++) { if (i + ib - 1 < mm) { mn = ib; } else { mn = mm - i + 1; } for (j = 1; j <= ib; j++) { work[j] = 0.0; } for (k = 1; k <= mn; k++) { work[k] = array[i, i + k - 1]; } for (j = 1; j <= ib; j++) { array[i, j] = work[j]; } } } //------------------------------------------------------------------------------ private void IERROR1(string fnameW) { //結果出力(Appendモード) System.IO.StreamWriter sw; string dat; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "area of element is less or equal to 0"; sw.WriteLine(dat); sw.Close(); MessageBox.Show("要素面積<=0", "異常終了通知"); this.Close(); } //------------------------------------------------------------------------------ private void IERROR2(string fnameW) { //結果出力(Appendモード) System.IO.StreamWriter sw; string dat; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "diagonal value is less or equal to 0"; sw.WriteLine(dat); sw.Close(); MessageBox.Show("対角項<=0", "異常終了通知"); this.Close(); } //------------------------------------------------------------------------------ } }