Option Explicit On Option Strict On Public Class Form1 Private pi As Double = Math.PI Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load 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 = "" ListView1.Clear() ListView1.View = View.Details ListView1.Columns.Add("nnn", 50, HorizontalAlignment.Right) ListView1.Columns.Add("dtest", 100, HorizontalAlignment.Right) ListView1.Columns.Add("ftest", 100, HorizontalAlignment.Right) End Sub Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click ListView1.Clear() ListView1.View = View.Details ListView1.Columns.Add("nnn", 50, HorizontalAlignment.Right) ListView1.Columns.Add("dtest", 100, HorizontalAlignment.Right) ListView1.Columns.Add("ftest", 100, HorizontalAlignment.Right) Call MAIN() End Sub Private Sub MAIN() '2次元no-tension軸対称応力解析 Dim i As Integer : Dim j As Integer : Dim k As Integer Dim kk As Integer : Dim ne As Integer : Dim n As Integer Dim ia As Integer : Dim ja As Integer Dim sr As System.IO.StreamReader Dim sw As System.IO.StreamWriter Dim dat As String Dim sbuf() As String Dim delim() As Char = {","c} Dim fnameR As String = "" Dim fnameW As String = "" Dim strcom As String '書出用コメント Dim NODT As Integer '節点総数 Dim NELT As Integer '要素総数 Dim MATEL As Integer '材料種類数 Dim KOZ As Integer 'x方向拘束節点数 Dim KOR As Integer 'y方向拘束節点数 Dim NF As Integer '載荷節点数 Dim nt As Integer '全自由度(総節点数×2(1節点自由度)) Dim mm As Integer '拘束節点処理後自由度(逆行列の寸法) Dim kakom(,) As Integer '要素構成接点番号 Dim Em() As Double '要素弾性係数 Dim po() As Double '要素ポアソン比 Dim alpha() As Double '要素線膨張係数 Dim ts() As Double '要素引張強さ Dim matno() As Integer '材料種別No Dim wm(,) As Double '材料物性作業用 Dim z() As Double '節点z座標 Dim r() As Double '節点r座標 Dim deltaT() As Double '節点温度変化量(+:温度上昇) Dim nokz() As Integer 'z方向拘束節点番号 Dim nokr() As Integer 'r方向拘束節点番号 Dim f() As Double '全体外力ベクトル Dim index() As Integer '拘束節点処理のためのインデックス Dim sm(8, 8) As Double '要素剛性マトリックス Dim tsm(,) As Double '全体剛性マトリックス Dim dis() As Double '全体節点変位増分 Dim rdis() As Double Dim wdis() As Double Dim dist() As Double '全体節点累計変位 Dim reac() As Double '作業用配列 Dim fevec(8) As Double '要素荷重ベクトル Dim ftvec() As Double '全体荷重ベクトル Dim strsave(,,) As Double '応力記憶用 Dim noten(,) As Integer 'No-tension積分点識別 Dim ntsum As Integer 'No-tension積分点数 Dim icount As Integer Dim nnn As Integer '収束計算回数制御 Dim maxita As Integer = 5000 '最大繰り返し数 Dim dtest As Double '変位収束状況記憶 Dim ftest As Double '荷重収束状況記憶 Dim ib As Integer 'バンド幅 Dim sttime As Date '計算開始時刻 Dim entime As Date '計算終了時刻 Dim tspan As TimeSpan '計算時間 Dim elarea1 As Double '要素半面積 Dim elarea2 As Double '要素半面積 Dim fact0 As Double = 1.0E-30 '0判定 Dim IPR As Integer '----------------------------------------------------- '入出力データファイル指定 '----------------------------------------------------- OpenFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory() If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameR = OpenFileDialog1.FileName If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameW = SaveFileDialog1.FileName Label2.Text = fnameR Label4.Text = fnameW '----------------------------------------------------- '時間計測開始 '----------------------------------------------------- sttime = DateTime.Now Label6.Text = Now.ToString '***************************************************** 'データ入力 '***************************************************** sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default) '-------------------------------------------------------------------------------------- '書出用コメント入力 '-------------------------------------------------------------------------------------- dat = sr.ReadLine : sbuf = dat.Split(delim) : strcom = sbuf(0) '------------------------------------------------------------------------------------------ '節点数,要素数,材料種類数,z方向拘束点数,r方向拘束点数,載荷点数 '------------------------------------------------------------------------------------------ dat = sr.ReadLine : sbuf = dat.Split(delim) NODT = CInt(sbuf(0)) NELT = CInt(sbuf(1)) MATEL = CInt(sbuf(2)) KOZ = CInt(sbuf(3)) KOR = CInt(sbuf(4)) NF = CInt(sbuf(5)) IPR = CInt(sbuf(6)) '------------------------------------------------------------------------------------------ '配列寸法宣言 '------------------------------------------------------------------------------------------ ReDim kakom(NELT, 4) ReDim Em(NELT) ReDim po(NELT) ReDim alpha(NELT) ReDim ts(NELT) ReDim matno(NELT) ReDim wm(MATEL, 4) ReDim strsave(NELT, 4, 4) ReDim noten(NELT, 4) ReDim z(NODT) ReDim r(NODT) ReDim deltaT(NODT) ReDim nokz(KOZ) ReDim nokr(KOR) nt = NODT * 2 ReDim index(nt) ReDim tsm(nt, nt) ReDim f(nt) ReDim ftvec(nt) ReDim dis(nt) ReDim rdis(nt) ReDim wdis(nt) ReDim dist(nt) ReDim reac(nt) '------------------------------------------------------------------------------------------ '材料特性(Em,po,gamma,alpha,ts) '------------------------------------------------------------------------------------------ For i = 1 To MATEL dat = sr.ReadLine : sbuf = dat.Split(delim) For j = 1 To 4 wm(i, j) = CDbl(sbuf(j - 1)) Next j Next i '------------------------------------------------------------------------------------------ '要素構成節点と材料種別(node-1,node-2,node-3,node-4,matno) '------------------------------------------------------------------------------------------ For ne = 1 To NELT dat = sr.ReadLine : sbuf = dat.Split(delim) kakom(ne, 1) = CInt(sbuf(0)) kakom(ne, 2) = CInt(sbuf(1)) kakom(ne, 3) = CInt(sbuf(2)) kakom(ne, 4) = CInt(sbuf(3)) matno(ne) = CInt(sbuf(4)) For i = 1 To MATEL If i = matno(ne) Then Em(ne) = wm(i, 1) po(ne) = wm(i, 2) alpha(ne) = wm(i, 3) ts(ne) = wm(i, 4) End If Next i Next ne '------------------------------------------------------------------------------------------ '節点座標(z,r),温度変化量 '------------------------------------------------------------------------------------------ For i = 1 To NODT dat = sr.ReadLine : sbuf = dat.Split(delim) z(i) = CDbl(sbuf(0)) r(i) = CDbl(sbuf(1)) deltaT(i) = CDbl(sbuf(2)) Next i '------------------------------------------------------------------------------------------ '拘束節点番号 '------------------------------------------------------------------------------------------ For i = 1 To nt : rdis(i) = 0.0 : Next i If 0 < KOZ Then For i = 1 To KOZ dat = sr.ReadLine : sbuf = dat.Split(delim) nokz(i) = CInt(sbuf(0)) rdis(2 * nokz(i) - 1) = CDbl(sbuf(1)) Next i End If If 0 < KOR Then For i = 1 To KOR dat = sr.ReadLine : sbuf = dat.Split(delim) nokr(i) = CInt(sbuf(0)) rdis(2 * nokr(i)) = CDbl(sbuf(1)) Next i End If '------------------------------------------------------------------------------------------ '節点番号,z方向荷重,r方向荷重 '------------------------------------------------------------------------------------------ For i = 1 To nt : f(i) = 0.0 : Next i If 0 < NF Then For i = 1 To NF dat = sr.ReadLine : sbuf = dat.Split(delim) n = CInt(sbuf(0)) f(2 * n - 1) = CDbl(sbuf(1)) f(2 * n) = CDbl(sbuf(2)) Next i End If sr.Close() '***************************************************** '入力確認出力 '***************************************************** Call 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 To NELT 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 Then Call 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 Then Call IERROR1(fnameW) Next ne '***************************************************** 'メイン処理計算 '***************************************************** '----------------------------------------------------- '保存応力クリア '----------------------------------------------------- For ne = 1 To NELT For i = 1 To 3 : For kk = 1 To 4 : strsave(ne, kk, i) = 0.0 : Next kk, i For kk = 1 To 4 : noten(ne, kk) = 0 : Next kk Next ne '----------------------------------------------------- '全体剛性行列・全体荷重ベクトルクリア '----------------------------------------------------- For i = 1 To nt dist(i) = 0.0 wdis(i) = 0.0 ftvec(i) = 0.0 reac(i) = 0.0 For j = 1 To nt tsm(i, j) = 0.0 Next j, i '----------------------------------------------------- '要素剛性行列・荷重ベクトル組立 '----------------------------------------------------- For ne = 1 To NELT '要素剛性行列 Call SM_AX4(ne, kakom, z, r, Em, po, sm) Call ASMAT(ne, 4, 2, kakom, sm, tsm) Next ne '----------------------------------------------------- '外力ベクトルセット '----------------------------------------------------- For i = 1 To nt : ftvec(i) = 0.0 : Next i '荷重ベクトルクリア For i = 1 To nt ftvec(i) = f(i) Next i '----------------------------------------------------- '拘束節点処理 '----------------------------------------------------- For i = 1 To nt : index(i) = i : Next i If 0 < KOZ Then For i = 1 To KOZ : n = 2 * nokz(i) - 1 : index(n) = 0 : Next i End If If 0 < KOR Then For i = 1 To KOR : n = 2 * nokr(i) : index(n) = 0 : Next i End If mm = 0 For i = 1 To nt If index(i) > 0 Then mm = mm + 1 : index(mm) = i End If Next i For i = 1 To mm ia = index(i) ftvec(i) = ftvec(ia) For j = 1 To mm ja = index(j) tsm(i, j) = tsm(ia, ja) Next j, i '----------------------------------------------------- 'バンド幅計算とバンド行列記憶 '----------------------------------------------------- ib = IBAND(mm, tsm, fact0) Call BAND(mm, ib, tsm) '----------------------------------------------------- '対角項確認と異常終了通知 '----------------------------------------------------- For i = 1 To mm If Math.Abs(tsm(i, 1)) < fact0 Then Call IERROR2(fnameW) Next i '----------------------------------------------------- '連立一次方程式(コレスキー法) '----------------------------------------------------- Call CHBAND10(mm, ib, tsm) '三角行列 '***************************************************** 'No-tension収束計算 '***************************************************** nnn = 0 Do '***************************************************** nnn = nnn + 1 '***************************************************** '----------------------------------------------------- '連立一次方程式(コレスキー法:前進消去・後退代入) '----------------------------------------------------- For i = 1 To mm reac(i) = ftvec(i) Next i Call CHBAND11(mm, ib, tsm, reac) '----------------------------------------------------- '変位増分計算 '----------------------------------------------------- For i = 1 To nt : dis(i) = 0.0 : Next i For i = 1 To mm : ia = index(i) : dis(ia) = reac(i) : Next i For i = 1 To nt wdis(i) = wdis(i) + dis(i) dist(i) = wdis(i) + rdis(i) reac(i) = 0.0 Next i '----------------------------------------------------- '応力計算・内力ベクトル組立 '----------------------------------------------------- For ne = 1 To NELT Call CALST_AX4(ne, kakom, z, r, deltaT, dist, Em, po, alpha, ts, strsave, fevec, noten) Call ASVEC(ne, 4, 2, kakom, fevec, reac) Next ne '----------------------------------------------------- '内力ベクトル保存,不平衡力ベクトル算定 '----------------------------------------------------- For i = 1 To nt : ftvec(i) = 0.0 : Next i For i = 1 To nt ftvec(i) = f(i) - reac(i) Next i '----------------------------------------------------- '不平衡力ベクトル拘束節点処理 '----------------------------------------------------- For i = 1 To mm ia = index(i) ftvec(i) = ftvec(ia) dis(i) = dis(ia) dist(i) = dist(ia) Next i '----------------------------------------------------- '収束判定(変位増分・不平衡力絶対値) '----------------------------------------------------- 'No-tension積分点の数 ntsum = 0 For ne = 1 To NELT For kk = 1 To 4 ntsum = ntsum + Math.Abs(noten(ne, kk)) Next kk Next ne '変位増分絶対値・不平衡力絶対値総和 icount = 0 dtest = 0.0 : ftest = 0.0 For i = 1 To mm dtest = dtest + Math.Abs(dis(i)) '変位増分絶対値総和 ftest = ftest + Math.Abs(ftvec(i)) '不平衡力絶対値総和 If Math.Abs(dis(i) / dist(i)) < 0.000001 Then icount = icount + 1 Next i '収束計算状況のListViewへの書き込み Dim lvitem As ListViewItem = New ListViewItem(nnn.ToString) lvitem.SubItems.Add(dtest.ToString("0.000e+00")) lvitem.SubItems.Add(ftest.ToString("0.000e+00")) ListView1.Items.Add(lvitem) '***************************************************** If 2 <= nnn Then If ntsum = 0 Then Exit Do If icount = mm Then Exit Do End If If maxita <= nnn Then Exit Do '***************************************************** Loop '/* 出力用全変位調整 */ For i = 1 To nt : dis(i) = 0.0 : Next i For i = 1 To mm : ia = index(i) : dis(ia) = dist(i) : Next i For i = 1 To nt : dist(i) = dis(i) + rdis(i) : Next i '----------------------------------------------------- '時間計測完了 '----------------------------------------------------- entime = DateTime.Now tspan = entime - sttime Label8.Text = entime.ToString Label10.Text = tspan.TotalSeconds.ToString("0.00sec") '***************************************************** '結果出力 '***************************************************** Call OUTRESUL(fnameW, NODT, NELT, f, reac, z, r, dist, _ strsave, ts, noten, 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 = "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("計算完了", "通知") End Sub Private Sub INPDAT(ByRef fnameW As String, ByVal strcom As String, _ ByVal NODT As Integer, ByVal NELT As Integer, _ ByVal MATEL As Integer, ByVal KOZ As Integer, ByVal KOR As Integer, _ ByVal NF As Integer, ByVal z() As Double, ByRef r() As Double, _ ByRef deltaT() As Double, ByRef f() As Double, ByRef rdis() As Double, _ ByRef nokz() As Integer, ByRef nokr() As Integer, ByRef kakom(,) As Integer, _ ByRef Em() As Double, ByRef po() As Double, ByRef alpha() As Double, ByRef ts() As Double, ByRef matno() As Integer, _ ByVal IPR As Integer) '入力確認出力 Dim sw As System.IO.StreamWriter Dim dat As String Dim i As Integer Dim k As Integer Dim l As Integer Dim ne As Integer 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,x,y,fx,fy,fix-x,fix-y,rdis-x,rdis-y,deltaT" : sw.WriteLine(dat) For i = 1 To NODT k = 0 : l = 0 If 0 < KOZ Then For j = 1 To KOZ If i = nokz(j) Then k = 1 Next j End If If 0 < KOR Then For j = 1 To KOR If i = nokr(j) Then l = 1 Next j End If 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) Next i 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 To NELT 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 sw.WriteLine(dat) Next ne sw.Close() End Sub Private Sub OUTRESUL(ByVal fnameW As String, ByVal NODT As Integer, ByVal NELT As Integer, _ ByRef f() As Double, ByRef reac() As Double, _ ByRef z() As Double, ByRef r() As Double, ByRef dist() As Double, _ ByRef strsave(,,) As Double, ByRef ts() As Double, _ ByRef noten(,) As Integer, ByRef matno() As Integer, _ ByVal IPR As Integer) '結果出力(Appendモード) Dim i As Integer Dim ne As Integer Dim kk As Integer Dim sw As System.IO.StreamWriter Dim dat As String Dim sigz As Double Dim sigr As Double Dim sigt As Double Dim tauzr As Double Dim ps1 As Double Dim ps2 As Double Dim ang As Double Dim ntsum As Integer sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "*displacements and forces" : sw.WriteLine(dat) dat = "node,coord-z,coord-r,dis-z,dis-r,reac-x,reac-y,ftvec-x,ftvec-y" : sw.WriteLine(dat) '節点変位 For i = 1 To NODT 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) Next i '要素応力 dat = "*stresses" : sw.WriteLine(dat) dat = "element,kk,sig-z,sig-r,sig-t,tau-zr,ps1,ps2,ang,noten,matno" : sw.WriteLine(dat) Select Case IPR Case 0 For ne = 1 To NELT For kk = 1 To 4 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, ps1, ps2, 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 dat = dat & "," & matno(ne).ToString sw.WriteLine(dat) Next kk Next ne Case 1 For ne = 1 To NELT ntsum = 0 sigz = 0.0 : sigr = 0.0 : sigt = 0.0 : tauzr = 0.0 For kk = 1 To 4 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) Next kk sigz = 0.25 * sigz sigr = 0.25 * sigr sigt = 0.25 * sigt tauzr = 0.25 * tauzr PST_AX4(sigz, sigr, tauzr, ps1, ps2, 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 dat = dat & "," & matno(ne).ToString sw.WriteLine(dat) Next ne End Select sw.Close() End Sub Private Sub DMAX_AX4(ByVal ne As Integer, ByVal Eme As Double, _ ByVal pzr As Double, ByVal pot As Double, ByRef d(,) As Double) '応力-ひずみ関係行列(De) Dim i As Integer Dim j As Integer For i = 1 To 4 : For j = 1 To 4 : d(i, j) = 0.0 : Next j, i 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 To 4 For j = 1 To 4 d(i, j) = d(i, j) * Eme / (1.0 + pzr) / (1.0 - 2.0 * pzr) Next j, i End Sub Private Sub BMAT_AX4(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef z() As Double, ByRef r() As Double, ByRef bm(,) As Double, _ ByRef detJ As Double, ByVal a As Double, ByVal b As Double) 'ひずみ-変位関係行列(B) Dim i As Integer : Dim j As Integer : Dim k As Integer : Dim l As Integer Dim J11 As Double : Dim J12 As Double : Dim J21 As Double : Dim J22 As Double Dim dn1a As Double : Dim dn2a As Double : Dim dn3a As Double : Dim dn4a As Double Dim dn1b As Double : Dim dn2b As Double : Dim dn3b As Double : Dim dn4b As Double Dim N1 As Double : Dim N2 As Double : Dim N3 As Double : Dim N4 As Double Dim rr As Double 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 To 4 For j = 1 To 8 bm(i, j) = 0.0 Next j, i 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 To 4 For j = 1 To 8 bm(i, j) = bm(i, j) / detJ Next j, i End Sub Private Sub IPOINT4(ByVal kk As Integer, ByRef a As Double, ByRef b As Double) 'Gauss積分点 Select Case kk Case 1 : a = -0.5773502692 : b = -0.5773502692 Case 2 : a = 0.5773502692 : b = -0.5773502692 Case 3 : a = 0.5773502692 : b = 0.5773502692 Case 4 : a = -0.5773502692 : b = 0.5773502692 End Select End Sub Private Sub SM_AX4(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef z() As Double, ByRef r() As Double, _ ByRef Em() As Double, ByRef po() As Double, ByRef sm(,) As Double) '要素剛性行列 Dim kk As Integer Dim i As Integer : Dim j As Integer : Dim k As Integer Dim d(4, 4) As Double Dim bm(4, 8) As Double Dim str(4, 8) As Double Dim detJ As Double : Dim a As Double : Dim b As Double Dim Eme As Double Dim poe As Double Dim N1 As Double : Dim N2 As Double : Dim N3 As Double : Dim N4 As Double Dim rr As Double Dim i1 As Integer : Dim i2 As Integer : Dim i3 As Integer : Dim i4 As Integer '配列クリア For i = 1 To 8 For j = 1 To 8 sm(i, j) = 0.0 Next j, i '剛性行列(sm)=(B)T(D,B)*r*detJ Eme = Em(ne) poe = po(ne) Call DMAX_AX4(ne, Eme, poe, poe, d) i1 = kakom(ne, 1) : i2 = kakom(ne, 2) : i3 = kakom(ne, 3) : i4 = kakom(ne, 4) For kk = 1 To 4 Call IPOINT4(kk, a, b) Call BMAT_AX4(ne, kakom, z, r, bm, detJ, a, b) For i = 1 To 4 For j = 1 To 8 str(i, j) = 0.0 For k = 1 To 4 str(i, j) = str(i, j) + d(i, k) * bm(k, j) Next k, j, i 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 To 8 For j = 1 To 8 For k = 1 To 4 sm(i, j) = sm(i, j) + bm(k, i) * str(k, j) * rr * detJ Next k, j, i Next kk End Sub Private Sub CALST_AX4(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef z() As Double, ByRef r() As Double, ByRef deltaT() As Double, ByRef dist() As Double, _ ByRef Em() As Double, ByRef po() As Double, ByRef alpha() As Double, ByRef ts() As Double, _ ByRef strsave(,,) As Double, ByRef fevec() As Double, ByRef noten(,) As Integer) '応力計算 Dim kk As Integer Dim i As Integer Dim j As Integer Dim i1 As Integer Dim i2 As Integer Dim i3 As Integer Dim i4 As Integer Dim tem As Double Dim d(4, 4) As Double Dim bm(4, 8) As Double Dim stress(4) As Double : Dim epsilon(4) As Double : Dim eps0(4) As Double Dim ia As Integer Dim sigz As Double : Dim sigr As Double : Dim sigt As Double : Dim tauzr As Double Dim ps1 As Double : Dim ps2 As Double : Dim ang As Double Dim phi As Double Dim wd(8) As Double Dim detJ As Double : Dim a As Double : Dim b As Double Dim Eme As Double : Dim pzr As Double : Dim pot As Double Dim N1 As Double : Dim N2 As Double : Dim N3 As Double : Dim N4 As Double Dim rr As Double '要素座標節点変位増分 For i = 1 To 4 ia = kakom(ne, i) wd(2 * i - 1) = dist(2 * ia - 1) wd(2 * i) = dist(2 * ia) Next i i1 = kakom(ne, 1) : i2 = kakom(ne, 2) : i3 = kakom(ne, 3) : i4 = kakom(ne, 4) For kk = 1 To 4 Call IPOINT4(kk, a, b) Call BMAT_AX4(ne, kakom, z, r, bm, detJ, a, b) 'ひずみ計算 {epsilon}=(B){u} For i = 1 To 4 epsilon(i) = 0.0 For j = 1 To 8 epsilon(i) = epsilon(i) + bm(i, j) * wd(j) Next j, i '積分点温度ひずみ 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 Eme = Em(ne) Select Case noten(ne, kk) Case 0 : pzr = po(ne) : pot = po(ne) Case 1 : pzr = 0.0 : pot = po(ne) Case 2 : pzr = 0.0 : pot = po(ne) Case 3 : pzr = po(ne) : pot = 0.0 Case 4 : pzr = 0.0 : pot = 0.0 Case 5 : pzr = 0.0 : pot = 0.0 End Select Call DMAX_AX4(ne, Eme, pzr, pot, d) '応力計算 {stress}=(D){epsilon} For i = 1 To 4 stress(i) = 0.0 For j = 1 To 4 stress(i) = stress(i) + d(i, j) * (epsilon(j) - eps0(j)) Next j, i '現状応力保存 sigz = stress(1) sigr = stress(2) sigt = stress(3) tauzr = stress(4) Call PST_AX4(sigz, sigr, tauzr, ps1, ps2, ang) 'No-tension判定 noten(ne, kk) = 0 If (ts(ne) < ps1 And ps2 < ts(ne)) And sigt < ts(ne) Then noten(ne, kk) = 1 If (ts(ne) < ps1 And ts(ne) < ps2) And sigt < ts(ne) Then noten(ne, kk) = 2 If (ps1 < ts(ne) And ps2 < ts(ne)) And ts(ne) < sigt Then noten(ne, kk) = 3 If (ts(ne) < ps1 And ps2 < ts(ne)) And ts(ne) < sigt Then noten(ne, kk) = 4 If (ts(ne) < ps1 And ts(ne) < ps2) And ts(ne) < sigt Then noten(ne, kk) = 5 '引張成分除去 phi = -ang / 180.0 * pi If noten(ne, kk) = 1 Or noten(ne, kk) = 2 Or noten(ne, kk) = 4 Then 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) End If If noten(ne, kk) = 2 Then 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) End If If noten(ne, kk) = 3 Or noten(ne, kk) = 4 Or noten(ne, kk) = 5 Then sigt = 0.0 End If '作業用応力保存 strsave(ne, kk, 1) = sigz strsave(ne, kk, 2) = sigr strsave(ne, kk, 3) = sigt strsave(ne, kk, 4) = tauzr Next kk '内力ベクトル For i = 1 To 8 : fevec(i) = 0.0 : Next i For kk = 1 To 4 Call IPOINT4(kk, a, b) Call BMAT_AX4(ne, kakom, z, r, bm, 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 To 8 For j = 1 To 4 fevec(i) = fevec(i) + bm(j, i) * strsave(ne, kk, j) * rr * detJ Next j, i Next kk End Sub Private Sub PST_AX4(ByVal sigz As Double, ByVal sigr As Double, ByVal tauzr As Double, _ ByRef ps1 As Double, ByRef ps2 As Double, ByRef ang As Double) '主応力計算 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) If sigz = sigr Then If tauzr > 0 Then ang = 45.0 If tauzr < 0 Then ang = 135.0 If tauzr = 0.0 Then ang = 0.0 Else ang = 0.5 * Math.Atan(2.0 * tauzr / (sigz - sigr)) : ang = 180.0 / 3.141592654 * ang If (sigz > sigr) And (tauzr >= 0) Then ang = ang If (sigz > sigr) And (tauzr < 0) Then ang = ang + 180.0 If sigz < sigr Then ang = ang + 90.0 End If End Sub Private Sub ASMAT(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _ ByRef kakom(,) As Integer, ByRef sm(,) As Double, ByRef tsm(,) As Double) '剛性行列組立(nod:1要素節点数,nhen:1節点の自由度) Dim i As Integer : Dim j As Integer : Dim k As Integer : Dim l As Integer Dim ki As Integer : Dim kj As Integer : Dim ks As Integer : Dim js As Integer Dim kik As Integer : Dim kjl As Integer : Dim ksk As Integer : Dim jsl As Integer For i = 1 To nod For j = 1 To nod ki = (kakom(ne, i) - 1) * nhen kj = (kakom(ne, j) - 1) * nhen ks = (i - 1) * nhen js = (j - 1) * nhen For k = 1 To nhen For l = 1 To nhen kik = ki + k : kjl = kj + l : ksk = ks + k : jsl = js + l tsm(kik, kjl) = tsm(kik, kjl) + sm(ksk, jsl) Next l, k Next j, i End Sub Private Sub ASVEC(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _ ByRef kakom(,) As Integer, ByRef fevec() As Double, ByRef ftvec() As Double) '荷重ベクトル組立(nod:1要素節点数,nhen:1節点の自由度) Dim i As Integer : Dim k As Integer Dim ki As Integer : Dim ks As Integer : Dim kik As Integer : Dim ksk As Integer For i = 1 To nod ki = (kakom(ne, i) - 1) * nhen ks = (i - 1) * nhen For k = 1 To nhen kik = ki + k : ksk = ks + k ftvec(kik) = ftvec(kik) + fevec(ksk) Next k Next i End Sub Private Sub CHBAND10(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double) Dim i As Integer Dim j As Integer Dim k As Integer Dim mm As Integer Dim mn As Integer Dim z As Double array(1, 1) = Math.Sqrt(array(1, 1)) For j = 2 To m1 array(1, j) = array(1, j) / array(1, 1) Next j For i = 2 To n1 If n1 - i + 1 > m1 Then mm = m1 Else mm = n1 - i + 1 End If For j = 1 To mm z = array(i, j) If j <> m1 Then If m1 - j + 1 < i Then mn = m1 - j + 1 Else mn = i End If For k = 2 To mn z = z - array(i - k + 1, k) * array(i - k + 1, j + k - 1) Next k If j > 1 Then array(i, j) = z / array(i, 1) Else array(i, 1) = Math.Sqrt(z) End If Else array(i, j) = z / array(i, 1) End If Next j, i End Sub Private Sub CHBAND11(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double, ByRef vector() As Double) Dim i As Integer Dim j As Integer Dim mm As Integer Dim z As Double vector(1) = vector(1) / array(1, 1) For i = 2 To n1 z = vector(i) If i > m1 Then mm = m1 Else mm = i End If For j = 2 To mm z = z - array(i - j + 1, j) * vector(i - j + 1) Next j vector(i) = z / array(i, 1) Next i vector(n1) = vector(n1) / array(n1, 1) For i = n1 - 1 To 1 Step -1 z = vector(i) For j = 2 To m1 If i + j - 1 > n1 Then Exit For z = z - array(i, j) * vector(i + j - 1) Next j vector(i) = z / array(i, 1) Next i End Sub Private Function IBAND(ByVal mm As Integer, ByRef array(,) As Double, ByVal fact0 As Double) As Integer Dim i As Integer Dim j As Integer '----------------------------------------------------- 'バンド幅計算 '----------------------------------------------------- IBAND = 1 For i = 1 To mm - 1 For j = mm To i + 1 Step -1 If fact0 <= Math.Abs(array(i, j)) Then Exit For Next j If IBAND <= j - i + 1 Then IBAND = j - i + 1 Next i If mm < IBAND Then IBAND = mm End Function Private Sub BAND(ByVal mm As Integer, ByVal ib As Integer, ByRef array(,) As Double) Dim i As Integer Dim j As Integer Dim mn As Integer Dim work() As Double ReDim work(mm) '----------------------------------------------------- 'フルマトリックスのバンド化記憶 '----------------------------------------------------- For i = 2 To mm If i + ib - 1 < mm Then mn = ib Else mn = mm - i + 1 End If For j = 1 To ib work(j) = 0.0 Next j For k = 1 To mn work(k) = array(i, i + k - 1) Next k For j = 1 To ib array(i, j) = work(j) Next j Next i End Sub Private Sub IERROR1(ByVal fnameW As String) '結果出力(Appendモード) Dim sw As System.IO.StreamWriter Dim dat As String 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", "異常終了通知") Me.Close() End Sub Private Sub IERROR2(ByVal fnameW As String) '結果出力(Appendモード) Dim sw As System.IO.StreamWriter Dim dat As String 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", "異常終了通知") Me.Close() End Sub End Class