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 vcsFEM4nodCH { 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 = ""; label11.Text = "nt,mm,ib"; } //------------------------------------------------------------------------------ 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 NSTRES; //応力状態（0;平面歪み,1;平面応力） int NODT; //節点総数 int NELT; //要素総数 int MATEL; //材料種類数 int KOX; //x方向変位既知節点数 int KOY; //y方向変位既知節点数 int NF; //載荷節点数 int nt; //全自由度（総節点数×２[１節点自由度]） int mm; //既知節点変位処理後自由度（逆行列の寸法） int ib; //バンド幅 double fact0 = 1.0E-30; //0判定 double elarea1; //要素半面積 double elarea2; //要素半面積 int IPR; double[,] sm = new double[9, 9]; //要素剛性マトリックス double[] fevec = new double[9]; //要素荷重ベクトル 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); NSTRES = int.Parse(sbuf[0]); NODT = int.Parse(sbuf[1]); NELT = int.Parse(sbuf[2]); MATEL = int.Parse(sbuf[3]); KOX = int.Parse(sbuf[4]); KOY = int.Parse(sbuf[5]); NF = int.Parse(sbuf[6]); IPR = int.Parse(sbuf[7]); //------------------------------------------------------------------------------------------ //配列寸法宣言 //------------------------------------------------------------------------------------------ int[,] kakom = new int[NELT + 1, 5]; //要素構成節点番号 double[] Em = new double[NELT + 1]; //要素弾性係数 double[] po = new double[NELT + 1]; //要素ポアソン比 double[] t = new double[NELT + 1]; //要素板厚 double[] gamma = new double[NELT + 1]; //要素単位体積重量[質量ではない] double[] gkh = new double[NELT + 1]; //要素の水平方向物体力係数[重力加速度に対する比率] double[] gkv = new double[NELT + 1]; //要素の鉛直方向物体力係数[重力加速度に対する比率] double[] alpha = new double[NELT + 1]; //要素線膨張係数 int[] matno = new int[NELT + 1]; //材料種別No double[,] wm = new double[MATEL + 1, 9]; //材料物性作業用 double[, ,] strsave = new double[NELT + 1, 5, 4]; //応力記憶用 double[] x = new double[NODT + 1]; //節点x座標 double[] y = new double[NODT + 1]; //節点y座標 double[] deltaT = new double[NODT + 1]; //節点温度変化量[＋：温度上昇] int[] nokx = new int[KOX + 1]; //x方向変位既知節点番号 int[] noky = new int[KOY + 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[] ftemp = new double[nt + 1]; //全体温度ひずみ相当荷重ベクトル double[] fmass = new double[nt + 1]; //全体物体力ベクトル double[] fdis = new double[nt + 1]; //既知節点変位相当外力ベクトル double[] rdis = new double[nt + 1]; //既知節点変位による荷重成分格納 double[] dis = new double[nt + 1]; //全体節点変位増分 double[] reac = new double[nt + 1]; //作業用配列 //------------------------------------------------------------------------------------------ //材料特性（t,Em,po,gamma,gkh,gkv,alpha） //------------------------------------------------------------------------------------------ for (i = 1; i <= MATEL; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); for (j = 1; j <= 7; 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]) { 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]; } } } //------------------------------------------------------------------------------------------ //節点座標（x,y），温度変化量 //------------------------------------------------------------------------------------------ for (i = 1; i <= NODT; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); x[i] = double.Parse(sbuf[0]); y[i] = double.Parse(sbuf[1]); deltaT[i] = double.Parse(sbuf[2]); } //------------------------------------------------------------------------------------------ //変位既知節点番号，既知節点変位量 //------------------------------------------------------------------------------------------ for (i = 1; i <= nt; i++) { rdis[i] = 0.0; } if (0 < KOX) { for (i = 1; i <= KOX; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokx[i] = int.Parse(sbuf[0]); rdis[2 * nokx[i] - 1] = double.Parse(sbuf[1]); } } if (0 < KOY) { for (i = 1; i <= KOY; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); noky[i] = int.Parse(sbuf[0]); rdis[2 * noky[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, NSTRES, NODT, NELT, MATEL, KOX, KOY, NF, x, y, deltaT, f, nokx, noky, kakom, Em, po, t, gamma, gkh, gkv, alpha, matno, rdis,IPR); //要素面積確認と異常終了通知 for (ne = 1; ne <= NELT; ne++) { i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; elarea1 = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]); if (elarea1 < fact0) IERROR1(fnameW); i = kakom[ne, 1]; j = kakom[ne, 3]; k = kakom[ne, 4]; elarea2 = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]); if (elarea2 < fact0) IERROR1(fnameW); } //***************************************************** //メイン処理計算 //***************************************************** //----------------------------------------------------- //保存応力クリア //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { for (i = 1; i <= 3; i++) { for (kk = 1; kk <= 4; kk++) { strsave[ne, kk, i] = 0.0; } } } //----------------------------------------------------- //全体剛性行列・全体荷重ベクトルクリア //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftemp[i] = 0.0; fmass[i] = 0.0; ftvec[i] = 0.0; reac[i] = 0.0; for (j = 1; j <= nt; j++) { tsm[i, j] = 0.0; } } //----------------------------------------------------- //要素剛性行列・荷重ベクトル組立 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { //剛性行列 SM_PL4(ne, kakom, x, y, NSTRES, Em, po, t, sm); ASMAT(ne, 4, 2, kakom, sm, tsm); //温度歪み相当荷重ベクトル，温度歪み相当応力[strsave] TV_PL4(ne, kakom, x, y, deltaT, NSTRES, Em, po, t, alpha, fevec, strsave); ASVEC(ne, 4, 2, kakom, fevec, ftemp); //物体力ベクトル MV_PL4(ne, kakom, x, y, t, gamma, gkh, gkv, fevec); ASVEC(ne, 4, 2, kakom, fevec, fmass); } //----------------------------------------------------- //既知節点変位による荷重成分ベクトル作成 //----------------------------------------------------- for (i = 1; i <= nt; i++) { fdis[i] = 0.0; for (j = 1; j <= nt; j++) { fdis[i] = fdis[i] + rdis[j] * tsm[i, j]; } } //----------------------------------------------------- //荷重ベクトルセット //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftvec[i] = f[i] + fmass[i] + ftemp[i] - fdis[i]; } //----------------------------------------------------- //変位既知節点自由度処理 //----------------------------------------------------- for (i = 1; i <= nt; i++) { index[i] = i; reac[i] = ftvec[i]; } if (0 < KOX) { for (i = 1; i <= KOX; i++) { n = 2 * nokx[i] - 1; index[n] = 0; } } if (0 < KOY) { for (i = 1; i <= KOY; i++) { n = 2 * noky[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]; reac[i] = reac[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); //三角行列 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; dis[i] = dis[i] + rdis[i]; } for (ne = 1; ne <= NELT; ne++) { CALST_PL4(ne, kakom, x, y, dis, NSTRES, Em, po, t, strsave, fevec); ASVEC(ne, 4, 2, kakom, fevec, reac); //内力ベクトル組立 } //----------------------------------------------------- //時間計測完了 //----------------------------------------------------- entime = DateTime.Now; tspan = entime - sttime; label8.Text = entime.ToString(); label10.Text = tspan.TotalSeconds.ToString("0.00sec"); label11.Text = "nt=" + nt.ToString() + " " + "mm=" + mm.ToString() + " " + "ib=" + ib.ToString(); //***************************************************** //結果出力 //***************************************************** OUTRESUL(fnameW, NODT, NELT, ftvec, reac, x, y, dis, strsave, matno,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 = "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 NSTRES, int NODT, int NELT, int MATEL, int KOX, int KOY, int NF, double[] x, double[] y, double[] deltaT, double[] f, int[] nokx, int[] noky, int[,] kakom, double[] Em, double[] po, double[] t, double[] gamma, double[] gkh, double[] gkv, double[] alpha, int[] matno, double[] rdis,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 = "NSTRES,NODT,NELT,MATEL,KOX,KOY,NF,IPR"; sw.WriteLine(dat); dat = NSTRES.ToString("0"); dat = dat + "," + NODT.ToString("0"); dat = dat + "," + NELT.ToString("0"); dat = dat + "," + MATEL.ToString("0"); dat = dat + "," + KOX.ToString("0"); dat = dat + "," + KOY.ToString("0"); dat = dat + "," + NF.ToString("0"); dat = dat + "," + IPR.ToString("0"); sw.WriteLine(dat); dat = "*node characteristics"; sw.WriteLine(dat); dat = "node,x,y,fx,fy,fix-x,fix-y,rdis-x,rdis-y,deltaT"; sw.WriteLine(dat); for (i = 1; i <= NODT; i++) { k = 0; l = 0; if (0 < KOX) { for (j = 1; j <= KOX; j++) { if (i == nokx[j]) k = 1; } } if (0 < KOY) { for (j = 1; j <= KOY; j++) { if (i == noky[j]) l = 1; } } dat = i.ToString("0"); dat = dat + "," + x[i].ToString("E"); dat = dat + "," + y[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 = "*node characteristics"; sw.WriteLine(dat); dat = "element,node-1,node-2,node-3,node-4,E,po,t,gamma,kh,kv,alpha,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 + "," + t[ne].ToString("E"); dat = dat + "," + gamma[ne].ToString("E"); dat = dat + "," + gkh[ne].ToString("E"); dat = dat + "," + gkv[ne].ToString("E"); dat = dat + "," + alpha[ne].ToString("E"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } sw.Close(); } //------------------------------------------------------------------------------ private void OUTRESUL(string fnameW, int NODT, int NELT, double[] ftvec, double[] reac, double[] x, double[] y, double[] dis, double[, ,] strsave, int[] matno,int IPR) { //結果出力（Appendモード） int i, ne, kk; System.IO.StreamWriter sw; string dat; double sigx, sigy, tauxy, ps1, ps2, ang; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "*displacements and forces"; sw.WriteLine(dat); dat = "node,coord-x,coord-y,dis-x,dis-y,reac-x,reac-y,ftvec-x,ftvec-y"; sw.WriteLine(dat); //節点変位 for (i = 1; i <= NODT; i++) { dat = i.ToString("0"); dat = dat + "," + x[i].ToString("E"); dat = dat + "," + y[i].ToString("E"); dat = dat + "," + dis[2 * i - 1].ToString("E"); dat = dat + "," + dis[2 * i].ToString("E"); dat = dat + "," + reac[2 * i - 1].ToString("E"); dat = dat + "," + reac[2 * i].ToString("E"); dat = dat + "," + ftvec[2 * i - 1].ToString("E"); dat = dat + "," + ftvec[2 * i].ToString("E"); sw.WriteLine(dat); } //要素応力 dat = "*stresses"; sw.WriteLine(dat); dat = "element,kk,sig-x,sig-y,tau-xy,ps1,ps2,ang,matno"; sw.WriteLine(dat); switch(IPR){ case 0: for (ne = 1; ne <= NELT; ne++) { for (kk = 1; kk <= 4; kk++) { sigx = strsave[ne, kk, 1]; sigy = strsave[ne, kk, 2]; tauxy = strsave[ne, kk, 3]; PST_PL(sigx, sigy, tauxy, out ps1, out ps2, out ang); dat = ne.ToString("0"); dat = dat + "," + kk.ToString("0"); dat = dat + "," + sigx.ToString("E"); dat = dat + "," + sigy.ToString("E"); dat = dat + "," + tauxy.ToString("E"); dat = dat + "," + ps1.ToString("E"); dat = dat + "," + ps2.ToString("E"); dat = dat + "," + ang.ToString("E"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } } break; case 1: for (ne = 1; ne <= NELT; ne++) { sigx = 0.0; sigy = 0.0; tauxy = 0.0; for (kk = 1; kk <= 4; kk++) { sigx = sigx+strsave[ne, kk, 1]; sigy = sigy+strsave[ne, kk, 2]; tauxy = tauxy+strsave[ne, kk, 3]; } sigx = 0.25 * sigx; sigy = 0.25 * sigy; tauxy = 0.25 * tauxy; PST_PL(sigx, sigy, tauxy, out ps1, out ps2, out ang); dat = ne.ToString("0"); dat = dat + "," + kk.ToString("0"); dat = dat + "," + sigx.ToString("E"); dat = dat + "," + sigy.ToString("E"); dat = dat + "," + tauxy.ToString("E"); dat = dat + "," + ps1.ToString("E"); dat = dat + "," + ps2.ToString("E"); dat = dat + "," + ang.ToString("E"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } break; } sw.Close(); } //------------------------------------------------------------------------------ private void DMAT_PL(int ne, int NSTRES, double[] Em, double[] po, double[,] d) { //応力-ひずみ関係行列[De] int i, j; for (i = 1; i <= 3; i++) { for (j = 1; j <= 3; j++) { d[i, j] = 0.0; } } switch (NSTRES) { //平面ひずみ case 0: d[1, 1] = Em[ne] * (1.0 - po[ne]) / ((1.0 + po[ne]) * (1.0 - 2.0 * po[ne])); d[1, 2] = Em[ne] * po[ne] / ((1.0 + po[ne]) * (1.0 - 2.0 * po[ne])); d[2, 1] = d[1, 2]; d[2, 2] = d[1, 1]; d[3, 3] = Em[ne] * 0.5 / (1.0 + po[ne]); break; //平面応力 case 1: d[1, 1] = Em[ne] / (1.0 - po[ne] * po[ne]); d[1, 2] = po[ne] * d[1, 1]; d[2, 1] = d[1, 2]; d[2, 2] = d[1, 1]; d[3, 3] = Em[ne] * 0.5 / (1.0 + po[ne]); break; } } //------------------------------------------------------------------------------ private void BMAT_PL4(int ne, int[,] kakom, double[] x, double[] y, double[,] bm, out double detJ, double a, double b) { //ひずみ-変位関係行列[B] int i, j, k, l; 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]; //[dN/da dN/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 * x[i] + dn2a * x[j] + dn3a * x[k] + dn4a * x[l]; J12 = dn1a * y[i] + dn2a * y[j] + dn3a * y[k] + dn4a * y[l]; J21 = dn1b * x[i] + dn2b * x[j] + dn3b * x[k] + dn4b * x[l]; J22 = dn1b * y[i] + dn2b * y[j] + dn3b * y[k] + dn4b * y[l]; detJ = J11 * J22 - J12 * J21; //[B]行列要素[dN/dx dN/dy] for (i = 1; i <= 3; 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, 1] = -J21 * dn1a + J11 * dn1b; bm[3, 2] = J22 * dn1a - J12 * dn1b; bm[3, 3] = -J21 * dn2a + J11 * dn2b; bm[3, 4] = J22 * dn2a - J12 * dn2b; bm[3, 5] = -J21 * dn3a + J11 * dn3b; bm[3, 6] = J22 * dn3a - J12 * dn3b; bm[3, 7] = -J21 * dn4a + J11 * dn4b; bm[3, 8] = J22 * dn4a - J12 * dn4b; for (i = 1; i <= 3; i++) { for (j = 1; j <= 8; j++) { bm[i, j] = bm[i, j] / detJ; } } } //------------------------------------------------------------------------------ private 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; } } //------------------------------------------------------------------------------ private void SM_PL4(int ne, int[,] kakom, double[] x, double[] y, int NSTRES, double[] Em, double[] po, double[] t, double[,] sm) { //要素剛性行列 int kk, i, j, k; double[,] d = new double[4, 4]; double[,] bm = new double[4, 9]; double[,] str = new double[4, 9]; double detJ, a, b; //配列クリア for (i = 1; i <= 8; i++) { for (j = 1; j <= 8; j++) { sm[i, j] = 0.0; } } //剛性行列[sm]=[B]T[D][B]*t*detJ if (NSTRES == 0) t[ne] = 1.0; DMAT_PL(ne, NSTRES, Em, po, d); for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); BMAT_PL4(ne, kakom, x, y, bm, out detJ, a, b); for (i = 1; i <= 3; i++) { for (j = 1; j <= 8; j++) { str[i, j] = 0.0; for (k = 1; k <= 3; k++) { str[i, j] = str[i, j] + d[i, k] * bm[k, j]; } } } for (i = 1; i <= 8; i++) { for (j = 1; j <= 8; j++) { for (k = 1; k <= 3; k++) { sm[i, j] = sm[i, j] + bm[k, i] * str[k, j] * t[ne] * detJ; } } } } } //------------------------------------------------------------------------------ private void TV_PL4(int ne, int[,] kakom, double[] x, double[] y, double[] deltaT, int NSTRES, double[] Em, double[] po, double[] t, double[] alpha, double[] fevec, double[, ,] strsave) { //温度ひずみ相当荷重ベクトル int kk, i, j, k, l; double[,] d = new double[4, 4]; double[,] bm = new double[4, 9]; double[] w1 = new double[4]; double[] w2 = new double[9]; double tempa; double detJ, a, b; double tm1, tm2, tm3, tm4; double nm1, nm2, nm3, nm4; i = kakom[ne, 1]; tm1 = deltaT[i]; j = kakom[ne, 2]; tm2 = deltaT[j]; k = kakom[ne, 3]; tm3 = deltaT[k]; l = kakom[ne, 4]; tm4 = deltaT[l]; //配列クリア for (i = 1; i <= 8; i++) { fevec[i] = 0.0; } if (NSTRES == 0) t[ne] = 1.0; DMAT_PL(ne, NSTRES, Em, po, d); for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); nm1 = 0.25 * (1.0 - a) * (1.0 - b); nm2 = 0.25 * (1.0 + a) * (1.0 - b); nm3 = 0.25 * (1.0 + a) * (1.0 + b); nm4 = 0.25 * (1.0 - a) * (1.0 + b); tempa = nm1 * tm1 + nm2 * tm2 + nm3 * tm3 + nm4 * tm4; switch (NSTRES) { case 0: //平面歪み w2[1] = tempa * (1.0 + po[ne]) * alpha[ne]; w2[2] = w2[1]; w2[3] = 0.0; break; case 1: //平面応力 w2[1] = tempa * alpha[ne]; w2[2] = w2[1]; w2[3] = 0.0; break; } for (i = 1; i <= 3; i++) { w1[i] = 0.0; for (j = 1; j <= 3; j++) { w1[i] = w1[i] + d[i, j] * w2[j]; } } //温度ひずみ相当応力 strsave[ne, kk, 1] = w1[1]; strsave[ne, kk, 2] = w1[2]; strsave[ne, kk, 3] = w1[3]; //温度ひずみ相当荷重ベクトル{f}=[B][D]{ε0}*t*detJ BMAT_PL4(ne, kakom, x, y, bm, out detJ, a, b); for (i = 1; i <= 8; i++) { w2[i] = 0.0; for (j = 1; j <= 3; j++) { w2[i] = w2[i] + bm[j, i] * w1[j]; } } for (i = 1; i <= 8; i++) { fevec[i] = fevec[i] + w2[i] * t[ne] * detJ; } } } //------------------------------------------------------------------------------ private void MV_PL4(int ne, int[,] kakom, double[] x, double[] y, double[] t, double[] gamma, double[] gkh, double[] gkv, double[] fevec) { //物体力ベクトル int i, j, k, l, ii, jj, kk, ll; double[,] nm = new double[3, 9]; double[,] ntn = new double[9, 9]; double[] w = new double[9]; double detJ, a, b; double J11, J12, J21, J22; double dn1a, dn2a, dn3a, dn4a, dn1b, dn2b, dn3b, dn4b; for (ii = 1; ii <= 8; ii++) { for (jj = 1; jj <= 8; jj++) { ntn[ii, jj] = 0.0; } } i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; l = kakom[ne, 4]; for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); //[dN/da dN/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 * x[i] + dn2a * x[j] + dn3a * x[k] + dn4a * x[l]; J12 = dn1a * y[i] + dn2a * y[j] + dn3a * y[k] + dn4a * y[l]; J21 = dn1b * x[i] + dn2b * x[j] + dn3b * x[k] + dn4b * x[l]; J22 = dn1b * y[i] + dn2b * y[j] + dn3b * y[k] + dn4b * y[l]; detJ = J11 * J22 - J12 * J21; for (ii = 1; ii <= 2; ii++) { for (jj = 1; jj <= 8; jj++) { nm[ii, jj] = 0.0; } } nm[1, 1] = 0.25 * (1.0 - a) * (1.0 - b); nm[2, 2] = nm[1, 1]; nm[1, 3] = 0.25 * (1.0 + a) * (1.0 - b); nm[2, 4] = nm[1, 3]; nm[1, 5] = 0.25 * (1.0 + a) * (1.0 + b); nm[2, 6] = nm[1, 5]; nm[1, 7] = 0.25 * (1.0 - a) * (1.0 + b); nm[2, 8] = nm[1, 7]; for (ii = 1; ii <= 8; ii++) { for (jj = 1; jj <= 8; jj++) { for (ll = 1; ll <= 2; ll++) { ntn[ii, jj] = ntn[ii, jj] + nm[ll, ii] * nm[ll, jj] * gamma[ne] * t[ne] * detJ; } } } } //物体力{ρg} w[1] = gkh[ne]; w[2] = gkv[ne]; w[3] = gkh[ne]; w[4] = gkv[ne]; w[5] = gkh[ne]; w[6] = gkv[ne]; w[7] = gkh[ne]; w[8] = gkv[ne]; for (i = 1; i <= 8; i++) { fevec[i] = 0.0; for (j = 1; j <= 8; j++) { fevec[i] = fevec[i] + ntn[i, j] * w[j]; } } } //------------------------------------------------------------------------------ private void CALST_PL4(int ne, int[,] kakom, double[] x, double[] y, double[] dis, int NSTRES, double[] Em, double[] po, double[] t, double[, ,] strsave, double[] fevec) { //応力計算 int kk, i, j, ia; double detJ, a, b; double[,] d = new double[4, 4]; double[,] bm = new double[4, 9]; double[] epsilon = new double[4]; double[] stress = new double[4]; double[] wd = new double[9]; //要素座標節点変位 for (i = 1; i <= 4; i++) { ia = kakom[ne, i]; wd[2 * i - 1] = dis[2 * ia - 1]; wd[2 * i] = dis[2 * ia]; } //要素応力 [sigma]=[D][B]{u} if (NSTRES == 0) t[ne] = 1.0; DMAT_PL(ne, NSTRES, Em, po, d); for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); BMAT_PL4(ne, kakom, x, y, bm, out detJ, a, b); for (i = 1; i <= 3; i++) { epsilon[i] = 0.0; for (j = 1; j <= 8; j++) { epsilon[i] = epsilon[i] + bm[i, j] * wd[j]; } } for (i = 1; i <= 3; i++) { stress[i] = 0.0; for (j = 1; j <= 3; j++) { stress[i] = stress[i] + d[i, j] * epsilon[j]; } } //応力算定(変位からの算定応力−温度歪み相当応力) strsave[ne, kk, 1] = stress[1] - strsave[ne, kk, 1]; strsave[ne, kk, 2] = stress[2] - strsave[ne, kk, 2]; strsave[ne, kk, 3] = stress[3] - strsave[ne, kk, 3]; } //内力ベクトル for (i = 1; i <= 8; i++) { fevec[i] = 0.0; } for (kk = 1; kk <= 4; kk++) { IPOINT4(kk, out a, out b); BMAT_PL4(ne, kakom, x, y, bm, out detJ, a, b); for (i = 1; i <= 8; i++) { for (j = 1; j <= 3; j++) { fevec[i] = fevec[i] + bm[j, i] * strsave[ne, kk, j] * t[ne] * detJ; } } } } //------------------------------------------------------------------------------ private void PST_PL(double sigx, double sigy, double tauxy, out double ps1, out double ps2, out double ang) { //主応力計算 ps1 = 0.5 * (sigx + sigy) + Math.Sqrt(0.25 * (sigx - sigy) * (sigx - sigy) + tauxy * tauxy); ps2 = 0.5 * (sigx + sigy) - Math.Sqrt(0.25 * (sigx - sigy) * (sigx - sigy) + tauxy * tauxy); ang = 0.0; if (sigx == sigy) { if (tauxy > 0) ang = 45.0; if (tauxy < 0) ang = 135.0; if (tauxy == 0.0) ang = 0.0; } else { ang = 0.5 * Math.Atan(2.0 * tauxy / (sigx - sigy)); ang = 180.0 / 3.141592654 * ang; if (sigx > sigy && tauxy >= 0) ang = ang + 0.0; if (sigx > sigy && tauxy < 0) ang = ang + 180.0; if (sigx < sigy) 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(); } //------------------------------------------------------------------------------ } }