Option Explicit On Option Strict On Public Class Form1 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 = "" Label11.Text = "nt,mm,ib" End Sub Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click Call MAIN() End Sub Private Sub MAIN() '2次元熱伝導解析 Dim i As Integer : Dim j As Integer : Dim k As Integer Dim n As Integer : Dim nnn As Integer Dim ia As Integer Dim ja As Integer Dim ne As Integer Dim it 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 fwnel As String = "" Dim ftinp As String = "" Dim strcom As String '書出用コメント Dim NODT As Integer '節点総数 Dim NELT As Integer '要素総数 Dim MATEL As Integer '材料種類数 Dim KOT As Integer '温度指定節点数 Dim KOC As Integer '熱伝達指定辺数 Dim nt As Integer '全自由度(総節点数×1[1節点自由度]) Dim mm As Integer '温度固定点処理後自由度 Dim ib As Integer 'バンド幅 Dim kakom(,) As Integer '要素構成接点番号 Dim Ak() As Double '要素熱伝導率 Dim Ac() As Double '要素比熱 Dim Arho() As Double '要素単位堆積重量 Dim Tk() As Double '要素最終温度上昇量 Dim Al() As Double '要素発熱速度係数 Dim matno() As Integer '材料種別No Dim wm(,) As Double '作業用材料物性記憶 Dim x() As Double '節点x座標 Dim y() As Double '節点y座標 Dim nokt() As Integer '温度指定節点番号 Dim nekc(,) As Integer '熱伝達境界要素・辺番号 Dim alphac() As Double '指定辺熱伝達率 Dim Tinp() As Double '指定温度入力 Dim Tcin1() As Double '指定熱伝達境界温度入力(前回) Dim Tcin2() As Double '指定熱伝達境界温度入力(現在) Dim alpc As Double '作業用熱伝達率 Dim tc As Double '作業用熱伝達境界温度 Dim kchen As Integer '作業用熱伝達境界辺番号 Dim dotq As Double '作業用発熱率 Dim tck1(,) As Double '全体剛性行列(前回) Dim tck2(,) As Double '全体剛性行列(現在) Dim rtck2(,) As Double '全体剛性行列(現在)拘束節点列保存 Dim ck1(4, 4) As Double '要素剛性行列(前回) Dim ck2(4, 4) As Double '要素剛性行列(現在) Dim ht(4, 4) As Double '要素熱伝達行列 Dim fq(4) As Double '要素発熱率ベクトル Dim fc(4) As Double '要素熱伝達ベクトル Dim ftvec() As Double '全体熱流束ベクトル Dim tempe0() As Double '初期全体節点温度ベクトル Dim tempe() As Double '全体節点温度ベクトル Dim reac() As Double '作業用節点温度ベクトル Dim index() As Integer '温度固定点処理用 Dim delta As Double '時間刻み Dim ttime As Double '時刻 Dim itmax As Integer '時刻歴ステップ数 Dim nout As Integer '出力節点数 Dim noout() As Integer '出力節点番号 Dim ntout As Integer '全節点温度出力回数 Dim nntout(0) As Integer '全節点温度出力ステップ Dim tempend(0, 0) As Double '出力用節点温度 Dim sttime As Date '計算開始時刻 Dim entime As Date '計算終了時刻 Dim tspan As TimeSpan '計算時間 Dim fact0 As Double = 1.0E-30 '0判定 Dim elarea1 As Double Dim elarea2 As Double '----------------------------------------------------- '基本データ読み込みファイル・出力ファイル指定 '----------------------------------------------------- OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory 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 = sttime.ToString '***************************************************** 'データ入力 '***************************************************** sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default) '-------------------------------------------------------------------------------------- '書出用コメント入力 '-------------------------------------------------------------------------------------- dat = sr.ReadLine : sbuf = dat.Split(delim) : strcom = sbuf(0) '-------------------------------------------------------------------------------------- '時刻歴読み込みファイル指定 '-------------------------------------------------------------------------------------- dat = sr.ReadLine : sbuf = dat.Split(delim) ftinp = System.IO.Path.GetDirectoryName(fnameR) & "\" & sbuf(0) '------------------------------------------------------------------------------------------ '節点数,要素数,材料種類数,温度固定点数,熱伝達境界辺数,時間刻み '------------------------------------------------------------------------------------------ dat = sr.ReadLine : sbuf = dat.Split(delim) NODT = CInt(sbuf(0)) NELT = CInt(sbuf(1)) MATEL = CInt(sbuf(2)) KOT = CInt(sbuf(3)) KOC = CInt(sbuf(4)) delta = CDbl(sbuf(5)) '------------------------------------------------------------------------------------------ '配列寸法宣言 '------------------------------------------------------------------------------------------ ReDim kakom(NELT, 4) ReDim Ak(NELT) ReDim Ac(NELT) ReDim Arho(NELT) ReDim Tk(NELT) ReDim Al(NELT) ReDim matno(NELT) ReDim wm(MATEL, 5) ReDim x(NODT) ReDim y(NODT) ReDim nokt(KOT) ReDim nekc(KOC, 2) ReDim alphac(KOC) ReDim Tinp(KOT) ReDim Tcin1(KOC) ReDim Tcin2(KOC) nt = NODT * 1 ReDim index(nt) ReDim tck1(nt, nt) ReDim tck2(nt, nt) ReDim rtck2(nt, KOT) ReDim tempe0(nt) ReDim tempe(nt) ReDim reac(nt) ReDim ftvec(nt) '------------------------------------------------------------------------------------------ '材料特性(Ak,Ac,Arho,Tk,Al) '------------------------------------------------------------------------------------------ For i = 1 To MATEL dat = sr.ReadLine : sbuf = dat.Split(delim) For j = 1 To 5 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 Ak(ne) = wm(i, 1) Ac(ne) = wm(i, 2) Arho(ne) = wm(i, 3) Tk(ne) = wm(i, 4) Al(ne) = wm(i, 5) End If Next i Next ne '------------------------------------------------------------------------------------------ '節点座標(x,y),初期節点温度 '------------------------------------------------------------------------------------------ For i = 1 To NODT dat = sr.ReadLine : sbuf = dat.Split(delim) x(i) = CDbl(sbuf(0)) y(i) = CDbl(sbuf(1)) tempe0(i) = CDbl(sbuf(2)) Next i '------------------------------------------------------------------------------------------ '温度指定節点番号 '------------------------------------------------------------------------------------------ If 0 < KOT Then For i = 1 To KOT dat = sr.ReadLine : sbuf = dat.Split(delim) : nokt(i) = CInt(sbuf(0)) Next i End If '------------------------------------------------------------------------------------------ '熱伝達境界指定要素・辺番号・熱伝達率 '------------------------------------------------------------------------------------------ If 0 < KOC Then For ne = 1 To KOC dat = sr.ReadLine : sbuf = dat.Split(delim) nekc(ne, 1) = CInt(sbuf(0)) '要素番号 nekc(ne, 2) = CInt(sbuf(1)) '辺番号 alphac(ne) = CDbl(sbuf(2)) '熱伝達率 Next ne End If '------------------------------------------------------------------------------------------ '出力節点数,出力節点番号 '------------------------------------------------------------------------------------------ dat = sr.ReadLine : sbuf = dat.Split(delim) nout = CInt(sbuf(0)) ReDim noout(nout) dat = sr.ReadLine : sbuf = dat.Split(delim) For i = 1 To nout noout(i) = CInt(sbuf(i - 1)) Next i '------------------------------------------------------------------------------------------ '全節点温度出力回数,出力ステップ番号 '------------------------------------------------------------------------------------------ dat = sr.ReadLine : sbuf = dat.Split(delim) ntout = CInt(sbuf(0)) If 0 < ntout Then ReDim nntout(ntout) ReDim tempend(NODT, ntout) dat = sr.ReadLine : sbuf = dat.Split(delim) For i = 1 To ntout nntout(i) = CInt(sbuf(i - 1)) Next i End If sr.Close() '--------------------------------------------------- '時刻歴入力ファイル空読み '--------------------------------------------------- sr = New System.IO.StreamReader(ftinp, System.Text.Encoding.Default) dat = sr.ReadLine() ' //コメント Do While sr.Peek() >= 0 dat = sr.ReadLine() : sbuf = dat.Split(delim) itmax = Integer.Parse(sbuf(0)) Loop sr.Close() '***************************************************** '入力確認出力 '***************************************************** Call INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOT, KOC, delta, nout, ntout, itmax, _ x, y, tempe, nokt, nekc, kakom, Ak, Ac, Arho, Tk, Al, matno) '要素面積確認と異常終了通知 For ne = 1 To NELT 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 Then Call 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 Then Call IERROR1(fnameW) Next ne '***************************************************** 'メイン処理計算 '***************************************************** '----------------------------------------------------- '全体行列クリア '----------------------------------------------------- For i = 1 To nt For j = 1 To nt tck1(i, j) = 0.0 tck2(i, j) = 0.0 Next j, i '----------------------------------------------------- '全体行列組立 '----------------------------------------------------- For ne = 1 To NELT Call STIFF(ne, kakom, x, y, Ak, Ac, Arho, ck1, ck2, delta) Call ASMAT(ne, 4, 1, kakom, ck1, tck1) '前回剛性-0.5*[K]+1/dt*[C] Call ASMAT(ne, 4, 1, kakom, ck2, tck2) '現在剛性+0.5*[K]+1/dt*[C] '熱伝達境界処理 If 0 < KOC Then For k = 1 To KOC If ne = nekc(k, 1) Then kchen = nekc(k, 2) alpc = alphac(k) Call HTMAT(ne, kchen, kakom, x, y, alpc, ht) '熱伝達行列 For i = 1 To 4 For j = 1 To 4 ht(i, j) = 0.5 * ht(i, j) Next j, i Call ASMAT(ne, 4, 1, kakom, ht, tck2) '現在剛性 For i = 1 To 4 For j = 1 To 4 ht(i, j) = -ht(i, j) Next j, i Call ASMAT(ne, 4, 1, kakom, ht, tck1) '前回剛性 End If Next k End If Next ne '----------------------------------------------------- '全体行列温度指定列要素保存 '----------------------------------------------------- If 0 < KOT Then For k = 1 To KOT For i = 1 To nt rtck2(i, k) = tck2(i, nokt(k)) Next i Next k End If '----------------------------------------------------- '温度指定節点処理 '----------------------------------------------------- For i = 1 To nt : index(i) = i : Next i If 0 < KOT Then For i = 1 To KOT : n = nokt(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) For j = 1 To mm ja = index(j) tck2(i, j) = tck2(ia, ja) Next j, i '----------------------------------------------------- 'バンド幅計算とフルマトリックスのバンド化記憶 '----------------------------------------------------- ib = IBAND(mm, tck2, fact0) Call BAND(mm, ib, tck2) '----------------------------------------------------- '対角項確認と異常終了通知 '----------------------------------------------------- For i = 1 To mm If Math.Abs(tck2(i, 1)) < fact0 Then Call IERROR2(fnameW) Next i '----------------------------------------------------- 'コレスキー法(三角行列) '----------------------------------------------------- Call CHBAND10(mm, ib, tck2) '三角行列 '----------------------------------------------------- '時刻歴出力準備 '----------------------------------------------------- '初期温度ベクトルセット For i = 1 To nt tempe(i) = tempe0(i) Next i '熱伝達境界初期値設定 If (0 < KOC) Then For k = 1 To KOC ne = nekc(k, 1) Select Case nekc(k, 2) Case 1 ia = kakom(ne, 1) ja = kakom(ne, 2) Case 2 ia = kakom(ne, 2) ja = kakom(ne, 3) Case 3 ia = kakom(ne, 3) ja = kakom(ne, 4) Case 4 ia = kakom(ne, 4) ja = kakom(ne, 1) End Select Tcin1(k) = 0.5 * (tempe(ia) + tempe(ja)) Next k End If sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "*temperature records of selected nodes" : sw.WriteLine(dat) '====================== it = 0 : ttime = 0.0 '====================== dat = "it,time" For i = 1 To nout dat = dat & ",Node-" & noout(i).ToString Next i sw.WriteLine(dat) dat = it.ToString("0") dat = dat & "," & ttime.ToString("E") For i = 1 To nout dat = dat & "," & tempe(noout(i)).ToString("0.000") Next i sw.WriteLine(dat) sr = New System.IO.StreamReader(ftinp, System.Text.Encoding.Default) dat = sr.ReadLine 'コメント '***************************************************** '時刻歴計算(時刻歴入力ファイル:ftinpの最終行まで) '***************************************************** it = 0 : nnn = 0 Do While sr.Peek() >= 0 it = it + 1 ttime = delta * CDbl(it) '----------------------------------------------------- '時刻歴入力(指定温度はその都度入力) '----------------------------------------------------- dat = sr.ReadLine : sbuf = dat.Split(delim) '1列目はダミー i = 0 '温度指定節点温度 If 0 < KOT Then For k = 1 To KOT : i = i + 1 : Tinp(k) = CDbl(sbuf(i)) : Next k End If '熱伝達境界辺温度 If 0 < KOC Then For k = 1 To KOC : i = i + 1 : Tcin2(k) = CDbl(sbuf(i)) : Next k End If '----------------------------------------------------- '熱流束ベクトル組立 '----------------------------------------------------- '前回温度ベクトル For i = 1 To nt reac(i) = 0.0 For j = 1 To nt reac(i) = reac(i) + tck1(i, j) * tempe(j) Next j Next i '温度指定境界処理 If 0 < KOT Then For k = 1 To KOT For i = 1 To nt reac(i) = reac(i) - Tinp(k) * rtck2(i, k) Next i Next k End If '熱伝達・発熱率ベクトル For i = 1 To nt ftvec(i) = 0.0 Next i For ne = 1 To NELT '発熱率ベクトル dotq = Arho(ne) * Ac(ne) * Tk(ne) * Al(ne) * Math.Exp(-Al(ne) * ttime) dotq = dotq + Arho(ne) * Ac(ne) * Tk(ne) * Al(ne) * Math.Exp(-Al(ne) * (ttime - delta)) dotq = 0.5 * dotq Call FQVEC(ne, kakom, x, y, dotq, fq) Call ASVEC(ne, 4, 1, kakom, fq, ftvec) '熱伝達ベクトル If 0 < KOC Then For k = 1 To KOC If ne = nekc(k, 1) Then kchen = nekc(k, 2) alpc = alphac(k) tc = 0.5 * (Tcin1(k) + Tcin2(k)) Call HTVEC(ne, kchen, kakom, x, y, alpc, tc, fc) Call ASVEC(ne, 4, 1, kakom, fc, ftvec) End If Next k End If Next ne '全体ベクトル For i = 1 To nt reac(i) = reac(i) + ftvec(i) Next i '----------------------------------------------------- '温度指定節点処理 '----------------------------------------------------- For i = 1 To mm ia = index(i) reac(i) = reac(ia) Next i '----------------------------------------------------- 'コレスキー法(前進消去・後退代入) '----------------------------------------------------- Call CHBAND11(mm, ib, tck2, reac) '----------------------------------------------------- '現在温度算定 '----------------------------------------------------- For i = 1 To nt : tempe(i) = 0.0 : Next i For i = 1 To mm : ia = index(i) : tempe(ia) = reac(i) : Next i For k = 1 To KOT tempe(nokt(k)) = Tinp(k) Next k '----------------------------------------------------- '記憶更新 '----------------------------------------------------- If 0 < KOC Then For k = 1 To KOC Tcin1(k) = Tcin2(k) Next k End If '----------------------------------------------------- '結果出力 '----------------------------------------------------- dat = it.ToString("0") dat = dat & "," & ttime.ToString("E") For i = 1 To nout dat = dat & "," & tempe(noout(i)).ToString("0.000") Next i sw.WriteLine(dat) '----------------------------------------------------- '節点温度保存 '----------------------------------------------------- If 0 < ntout Then For n = 1 To ntout If it = nntout(n) Then nnn = nnn + 1 For i = 1 To NODT tempend(i, nnn) = tempe(i) Next i End If Next n End If Loop sr.Close() '----------------------------------------------------- '節点温度出力 '----------------------------------------------------- If 0 < ntout Then dat = "*temperature of nodes for selected times" : sw.WriteLine(dat) dat = "Node,init" For i = 1 To ntout dat = dat & "," & nntout(i).ToString Next i sw.WriteLine(dat) For i = 1 To NODT dat = i.ToString dat = dat & "," & tempe0(i).ToString("0.000") For j = 1 To ntout dat = dat & "," & tempend(i, j).ToString("0.000") Next j sw.WriteLine(dat) Next i End If '----------------------------------------------------- '時間計測完了 '----------------------------------------------------- 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() 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 & " (VB)" sw.WriteLine(dat) sw.Close() 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 KOT As Integer, ByVal KOC As Integer, ByVal delta As Double, _ ByVal nout As Integer, ByVal ntout As Integer, ByVal itmax As Integer, _ ByVal x() As Double, ByRef y() As Double, ByRef tempe() As Double, _ ByRef nokt() As Integer, ByRef nekc(,) As Integer, _ ByRef kakom(,) As Integer, _ ByRef Ak() As Double, ByRef Ac() As Double, ByRef Arho() As Double, _ ByRef Tk() As Double, ByRef Al() As Double, ByRef matno() As Integer) '入力確認出力 Dim sw As System.IO.StreamWriter Dim dat As String Dim i As Integer Dim j As Integer Dim k As Integer Dim k1 As Integer Dim k2 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,KOT,KOC,delta,nout,ntout,itmax" : sw.WriteLine(dat) dat = NODT.ToString("0") dat = dat & "," & NELT.ToString("0") dat = dat & "," & MATEL.ToString("0") dat = dat & "," & KOT.ToString("0") dat = dat & "," & KOC.ToString("0") dat = dat & "," & delta.ToString("E") dat = dat & "," & nout.ToString("0") dat = dat & "," & ntout.ToString("0") dat = dat & "," & itmax.ToString("0") sw.WriteLine(dat) dat = "*node characteristics" : sw.WriteLine(dat) dat = "node,x,y,fix-T,Tini" : sw.WriteLine(dat) For i = 1 To NODT k = 0 If 0 < KOT Then For j = 1 To KOT If i = nokt(j) Then k = 1 : Exit For End If Next j End If dat = i.ToString("0") dat = dat & "," & x(i).ToString("E") dat = dat & "," & y(i).ToString("E") dat = dat & "," & k.ToString("0") dat = dat & "," & tempe(i).ToString("E") sw.WriteLine(dat) Next i dat = "*element characteristics" : sw.WriteLine(dat) dat = "element,node-1,node-2,node-3,node-4,k1,k2,Ak,Ac,Arho,Tk,Al,matno" : sw.WriteLine(dat) For ne = 1 To NELT k1 = 0 : k2 = 0 If 0 < KOC Then k1 = 0 : k2 = 0 For k = 1 To KOC If ne = nekc(k, 1) Then k1 = 1 '熱伝達境界を有する要素(k1=1) k2 = nekc(k, 2) '熱伝達境界の要素辺(1 or 2 or 3 or 4) End If Next k End If 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 & "," & k1.ToString("0") dat = dat & "," & k2.ToString("0") dat = dat & "," & Ak(ne).ToString("E") dat = dat & "," & Ac(ne).ToString("E") dat = dat & "," & Arho(ne).ToString("E") dat = dat & "," & Tk(ne).ToString("E") dat = dat & "," & Al(ne).ToString("E") dat = dat & "," & matno(ne).ToString sw.WriteLine(dat) Next ne sw.Close() End Sub Private Sub NMAT(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, _ ByRef dnk(,) As Double, ByRef dnc(,) As Double, _ ByRef detJ As Double, ByVal a As Double, ByVal b As Double) 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 dn1x As Double : Dim dn2x As Double : Dim dn3x As Double : Dim dn4x As Double Dim dn1y As Double : Dim dn2y As Double : Dim dn3y As Double : Dim dn4y As Double Dim nm1 As Double : Dim nm2 As Double : Dim nm3 As Double : Dim nm4 As Double 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 '[dN/dx dN/dy] dn1x = (J22 * dn1a - J12 * dn1b) / detJ dn2x = (J22 * dn2a - J12 * dn2b) / detJ dn3x = (J22 * dn3a - J12 * dn3b) / detJ dn4x = (J22 * dn4a - J12 * dn4b) / detJ dn1y = (-J21 * dn1a + J11 * dn1b) / detJ dn2y = (-J21 * dn2a + J11 * dn2b) / detJ dn3y = (-J21 * dn3a + J11 * dn3b) / detJ dn4y = (-J21 * dn4a + J11 * dn4b) / detJ '[N]行列要素 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) dnk(1, 1) = dn1x * dn1x + dn1y * dn1y dnk(1, 2) = dn1x * dn2x + dn1y * dn2y dnk(1, 3) = dn1x * dn3x + dn1y * dn3y dnk(1, 4) = dn1x * dn4x + dn1y * dn4y dnk(2, 1) = dn2x * dn1x + dn2y * dn1y dnk(2, 2) = dn2x * dn2x + dn2y * dn2y dnk(2, 3) = dn2x * dn3x + dn2y * dn3y dnk(2, 4) = dn2x * dn4x + dn2y * dn4y dnk(3, 1) = dn3x * dn1x + dn3y * dn1y dnk(3, 2) = dn3x * dn2x + dn3y * dn2y dnk(3, 3) = dn3x * dn3x + dn3y * dn3y dnk(3, 4) = dn3x * dn4x + dn3y * dn4y dnk(4, 1) = dn4x * dn1x + dn4y * dn1y dnk(4, 2) = dn4x * dn2x + dn4y * dn2y dnk(4, 3) = dn4x * dn3x + dn4y * dn3y dnk(4, 4) = dn4x * dn4x + dn4y * dn4y dnc(1, 1) = nm1 * nm1 dnc(1, 2) = nm1 * nm2 dnc(1, 3) = nm1 * nm3 dnc(1, 4) = nm1 * nm4 dnc(2, 1) = nm2 * nm1 dnc(2, 2) = nm2 * nm2 dnc(2, 3) = nm2 * nm3 dnc(2, 4) = nm2 * nm4 dnc(3, 1) = nm3 * nm1 dnc(3, 2) = nm3 * nm2 dnc(3, 3) = nm3 * nm3 dnc(3, 4) = nm3 * nm4 dnc(4, 1) = nm4 * nm1 dnc(4, 2) = nm4 * nm2 dnc(4, 3) = nm4 * nm3 dnc(4, 4) = nm4 * nm4 End Sub Private Sub IPOINT(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 STIFF(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, _ ByRef Ak() As Double, ByRef Ac() As Double, ByRef Arho() As Double, _ ByRef ck1(,) As Double, ByRef ck2(,) As Double, ByVal delta As Double) '要素剛性行列 Dim kk As Integer Dim i As Integer Dim j As Integer Dim dnk(4, 4) As Double Dim dnc(4, 4) As Double Dim ww1(4, 4) As Double Dim ww2(4, 4) As Double Dim detJ As Double : Dim a As Double : Dim b As Double For i = 1 To 4 For j = 1 To 4 ww1(i, j) = 0.0 ww2(i, j) = 0.0 ck1(i, j) = 0.0 ck2(i, j) = 0.0 Next j, i For kk = 1 To 4 Call IPOINT(kk, a, b) Call NMAT(ne, kakom, x, y, dnk, dnc, detJ, a, b) For i = 1 To 4 For j = 1 To 4 ww1(i, j) = ww1(i, j) + Ak(ne) * dnk(i, j) * detJ ww2(i, j) = ww2(i, j) + Arho(ne) * Ac(ne) * dnc(i, j) * detJ Next j, i Next kk For i = 1 To 4 For j = 1 To 4 ck1(i, j) = -0.5 * ww1(i, j) + 1.0 / delta * ww2(i, j) ck2(i, j) = +0.5 * ww1(i, j) + 1.0 / delta * ww2(i, j) Next j, i End Sub Private Sub LPOINT(ByVal kk As Integer, ByVal kchen As Integer, _ ByRef a As Double, ByRef b As Double) '要素辺上積分におけるGauss積分点 Select Case kchen Case 1 b = -1.0 If kk = 1 Then a = -0.5773502692 If kk = 2 Then a = 0.5773502692 Case 2 a = 1.0 If kk = 1 Then b = -0.5773502692 If kk = 2 Then b = 0.5773502692 Case 3 b = 1.0 If kk = 1 Then a = -0.5773502692 If kk = 2 Then a = 0.5773502692 Case 4 a = -1.0 If kk = 1 Then b = -0.5773502692 If kk = 2 Then b = 0.5773502692 End Select End Sub Private Sub HTMAT(ByVal ne As Integer, ByVal kchen As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, _ ByVal alpc As Double, ByRef ht(,) As Double) '要素熱伝達マトリックス Dim i As Integer Dim j As Integer Dim kk As Integer Dim a As Double Dim b As Double Dim s As Double Dim vn(4) As Double Select Case kchen Case 1 : i = kakom(ne, 1) : j = kakom(ne, 2) Case 2 : i = kakom(ne, 2) : j = kakom(ne, 3) Case 3 : i = kakom(ne, 3) : j = kakom(ne, 4) Case 4 : i = kakom(ne, 4) : j = kakom(ne, 1) End Select s = Math.Sqrt((x(j) - x(i)) ^ 2 + (y(j) - y(i)) ^ 2) For i = 1 To 4 For j = 1 To 4 ht(i, j) = 0.0 Next j, i For kk = 1 To 2 Call LPOINT(kk, kchen, a, b) '[N]行列要素 vn(1) = 0.25 * (1.0 - a) * (1.0 - b) vn(2) = 0.25 * (1.0 + a) * (1.0 - b) vn(3) = 0.25 * (1.0 + a) * (1.0 + b) vn(4) = 0.25 * (1.0 - a) * (1.0 + b) For i = 1 To 4 For j = 1 To 4 ht(i, j) = ht(i, j) + 0.5 * s * alpc * vn(i) * vn(j) Next j, i Next kk End Sub Private Sub FQVEC(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, _ ByVal dotq As Double, ByRef fq() As Double) '発熱率ベクトル Dim i As Integer Dim j As Integer Dim k As Integer Dim l As Integer Dim kk As Integer Dim a As Double Dim b As Double Dim vn(4) As Double Dim detJ As Double 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 For i = 1 To 4 fq(i) = 0.0 Next i For kk = 1 To 4 Call IPOINT(kk, a, b) 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 '[N]行列要素 vn(1) = 0.25 * (1.0 - a) * (1.0 - b) vn(2) = 0.25 * (1.0 + a) * (1.0 - b) vn(3) = 0.25 * (1.0 + a) * (1.0 + b) vn(4) = 0.25 * (1.0 - a) * (1.0 + b) For i = 1 To 4 fq(i) = fq(i) + dotq * vn(i) * detJ Next i Next kk End Sub Private Sub HTVEC(ByVal ne As Integer, ByVal kchen As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, _ ByVal alpc As Double, ByVal tc As Double, ByRef fc() As Double) '熱伝達ベクトル Dim i As Integer Dim j As Integer Dim kk As Integer Dim vn(4) As Double Dim a As Double Dim b As Double Dim s As Double Select Case kchen Case 1 : i = kakom(ne, 1) : j = kakom(ne, 2) Case 2 : i = kakom(ne, 2) : j = kakom(ne, 3) Case 3 : i = kakom(ne, 3) : j = kakom(ne, 4) Case 4 : i = kakom(ne, 4) : j = kakom(ne, 1) End Select s = Math.Sqrt((x(j) - x(i)) ^ 2 + (y(j) - y(i)) ^ 2) For i = 1 To 4 fc(i) = 0.0 Next i For kk = 1 To 2 Call LPOINT(kk, kchen, a, b) '[N]行列要素 vn(1) = 0.25 * (1.0 - a) * (1.0 - b) vn(2) = 0.25 * (1.0 + a) * (1.0 - b) vn(3) = 0.25 * (1.0 + a) * (1.0 + b) vn(4) = 0.25 * (1.0 - a) * (1.0 + b) For i = 1 To 4 fc(i) = fc(i) + 0.5 * s * alpc * tc * vn(i) Next i Next kk End Sub Private Function CALtempeNEL(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef tempe() As Double) As Double Dim i As Integer Dim j As Integer Dim k As Integer Dim l As Integer Dim nm1 As Double Dim nm2 As Double Dim nm3 As Double Dim nm4 As Double Dim a As Double Dim b As Double Dim tem As Double i = kakom(ne, 1) j = kakom(ne, 2) k = kakom(ne, 3) l = kakom(ne, 4) tem = 0.0 For kk = 1 To 4 Call IPOINT(kk, a, b) '[N]行列要素 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) tem = tem + nm1 * tempe(i) + nm2 * tempe(j) + nm3 * tempe(k) + nm4 * tempe(l) Next kk CALtempeNEL = 0.25 * tem End Function 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