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 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 NSTRES As Integer '応力状態(0:平面歪み,1:平面応力) Dim NODT As Integer '節点総数 Dim NELT As Integer '要素総数 Dim MATEL As Integer '材料種別数 Dim KOX As Integer 'x方向変位既知節点数 Dim KOY As Integer 'y方向変位既知節点数 Dim NF As Integer '載荷節点数 Dim nt As Integer '全自由度 Dim mm As Integer '変位既知節点処理後の自由度 Dim kakom(,) As Integer '要素構成接点番号 Dim Em() As Double '要素弾性係数 Dim po() As Double '要素ポアソン比 Dim t() As Double '要素板厚 Dim gamma() As Double '要素単位体積重量(質量ではない) Dim gkh() As Double '要素の水平方向物体力係数(重力加速度に対する比率) Dim gkv() As Double '要素の鉛直方向物体力係数(重力加速度に対する比率) Dim alpha() As Double '要素線膨張係数 Dim matno() As Integer '材料種別No Dim x() As Double '節点x座標 Dim y() As Double '節点y座標 Dim deltaT() As Double '節点温度変化量(+:温度上昇) Dim nokx() As Integer 'x方向変位既知節点番号 Dim noky() As Integer 'y方向変位既知節点番号 Dim f() As Double '全体節点外力ベクトル Dim ftemp() As Double '全体温度ひずみ相当荷重ベクトル Dim fmass() As Double '全体物体力ベクトル Dim rdis(i) As Double '既知節点変位量 Dim fdis(i) As Double '既知節点変位による荷重成分格納 Dim index() As Integer '既知節点変位処理のためのインデックス Dim sm(6, 6) As Double '要素剛性マトリックス Dim tsm(,) As Double '全体剛性マトリックス Dim dis() As Double '全体節点変位 Dim reac() As Double '作業用配列 Dim fevec(6) As Double '要素荷重ベクトル Dim ftvec() As Double '全体荷重ベクトル Dim strsave(,) As Double '応力記憶用 Dim wm(,) As Double '材料種別作業用 Dim sttime As Date '計算開始時刻 Dim entime As Date '計算終了時刻 Dim tspan As TimeSpan '計算時間 Dim ib As Integer 'バンド幅 Dim fact0 As Double = 1.0E-30 '0判定 '----------------------------------------------------- '入出力データファイル指定 '----------------------------------------------------- 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 = 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 = CInt(sbuf(0)) NODT = CInt(sbuf(1)) NELT = CInt(sbuf(2)) MATEL = CInt(sbuf(3)) KOX = CInt(sbuf(4)) KOY = CInt(sbuf(5)) NF = CInt(sbuf(6)) '-------------------------------------------------------------------------------------- '配列寸法宣言 '-------------------------------------------------------------------------------------- ReDim kakom(NELT, 3) ReDim Em(NELT) ReDim po(NELT) ReDim t(NELT) ReDim gamma(NELT) ReDim gkh(NELT) ReDim gkv(NELT) ReDim alpha(NELT) ReDim matno(NELT) ReDim wm(MATEL, 7) ReDim strsave(NELT, 3) ReDim x(NODT) ReDim y(NODT) ReDim deltaT(NODT) ReDim nokx(KOX) ReDim noky(KOY) nt = NODT * 2 ReDim index(nt) ReDim tsm(nt, nt) ReDim f(nt) ReDim ftvec(nt) ReDim ftemp(nt) ReDim fmass(nt) ReDim fdis(nt) ReDim rdis(nt) ReDim dis(nt) ReDim reac(nt) '------------------------------------------------------------------------------------------ '材料特性(t,Em,po,gamma,gkh,gkv,alpha) '------------------------------------------------------------------------------------------ For i = 1 To MATEL dat = sr.ReadLine : sbuf = dat.Split(delim) For j = 1 To 7 wm(i, j) = CDbl(sbuf(j - 1)) Next j Next i '------------------------------------------------------------------------------------------ '要素構成節点と材料種別(node-1,node-2,node-3,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)) matno(ne) = CInt(sbuf(3)) For i = 1 To MATEL If i = matno(ne) Then 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) End If Next i Next ne '-------------------------------------------------------------------------------------- '節点座標(x,y),温度変化量(deltaT) '-------------------------------------------------------------------------------------- For i = 1 To NODT dat = sr.ReadLine : sbuf = dat.Split(delim) x(i) = CDbl(sbuf(0)) y(i) = CDbl(sbuf(1)) deltaT(i) = CDbl(sbuf(2)) Next i '-------------------------------------------------------------------------------------- '変位既知節点番号,既知節点変位量 '-------------------------------------------------------------------------------------- For i = 1 To nt : rdis(i) = 0.0 : Next i If 0 < KOX Then For i = 1 To KOX dat = sr.ReadLine : sbuf = dat.Split(delim) nokx(i) = CInt(sbuf(0)) rdis(2 * nokx(i) - 1) = CDbl(sbuf(1)) Next i End If If 0 < KOY Then For i = 1 To KOY dat = sr.ReadLine : sbuf = dat.Split(delim) noky(i) = CInt(sbuf(0)) rdis(2 * noky(i)) = CDbl(sbuf(1)) Next i End If '-------------------------------------------------------------------------------------- '節点番号,x方向荷重,y方向荷重 '-------------------------------------------------------------------------------------- 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, NSTRES, NODT, NELT, MATEL, KOX, KOY, NF, x, y, f, _ deltaT, nokx, noky, kakom, Em, po, t, gamma, gkh, gkv, alpha, matno, rdis) '要素面積確認と異常終了通知 For ne = 1 To NELT If CALA_PL3(ne, kakom, x, y) < fact0 Then Call IERROR1(fnameW) Next ne '***************************************************** 'メイン処理計算 '***************************************************** '----------------------------------------------------- '全体剛性行列・全体荷重ベクトルクリア '----------------------------------------------------- For i = 1 To nt ftemp(i) = 0.0 fmass(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_PL3(ne, kakom, x, y, NSTRES, Em, po, t, sm) Call ASMAT(ne, 3, 2, kakom, sm, tsm) '温度ひずみ相当ベクトル,温度ひずみ相当応力保存(strsave) Call TV_PL3(ne, kakom, x, y, deltaT, NSTRES, Em, po, t, alpha, fevec, strsave) Call ASVEC(ne, 3, 2, kakom, fevec, ftemp) '物体力ベクトル Call MV_PL3(ne, kakom, x, y, t, gamma, gkh, gkv, fevec) Call ASVEC(ne, 3, 2, kakom, fevec, fmass) Next ne '----------------------------------------------------- '既知節点変位による荷重成分ベクトル作成 '----------------------------------------------------- For i = 1 To nt fdis(i) = 0.0 For j = 1 To nt fdis(i) = fdis(i) + rdis(j) * tsm(i, j) Next j, i '----------------------------------------------------- '荷重ベクトルセット '----------------------------------------------------- For i = 1 To nt ftvec(i) = f(i) + fmass(i) + ftemp(i) - fdis(i) Next i '----------------------------------------------------- '変位既知節点自由度処理 '----------------------------------------------------- For i = 1 To nt index(i) = i reac(i) = ftvec(i) Next i If 0 < KOX Then For i = 1 To KOX : n = 2 * nokx(i) - 1 : index(n) = 0 : Next i End If If 0 < KOY Then For i = 1 To KOY : n = 2 * noky(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) reac(i) = reac(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) '三角行列 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 : reac(i) = 0.0 : dis(i) = dis(i) + rdis(i) : Next i '----------------------------------------------------- '応力計算 {strsave}={sig}-{strsave} '----------------------------------------------------- For ne = 1 To NELT Call CALST_PL3(ne, kakom, x, y, dis, NSTRES, Em, po, t, strsave, fevec) Call ASVEC(ne, 3, 2, kakom, fevec, reac) '内力ベクトル組立 Next ne '----------------------------------------------------- '時間計測完了 '----------------------------------------------------- entime = DateTime.Now tspan = entime - sttime Label8.Text = entime.ToString Label10.Text = tspan.TotalSeconds.ToString("0.00sec") Label11.Text = "NODT=" & NODT.ToString() & " nt=" & nt.ToString() & " mm=" & mm.ToString() & " ib=" & ib.ToString() '***************************************************** '結果出力 '***************************************************** Call OUTRESUL(fnameW, NODT, NELT, ftvec, reac, x, y, dis, strsave, matno) 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("計算完了", "通知") End Sub Private Sub INPDAT(ByRef fnameW As String, ByVal strcom As String, _ ByVal NSTRES As Integer, ByVal NODT As Integer, ByVal NELT As Integer, _ ByVal MATEL As Integer, ByVal KOX As Integer, ByVal KOY As Integer, ByVal NF As Integer, _ ByVal x() As Double, ByRef y() As Double, ByRef f() As Double, ByRef deltaT() As Double, _ ByRef nokx() As Integer, ByRef noky() As Integer, ByRef kakom(,) As Integer, _ ByRef Em() As Double, ByRef po() As Double, ByRef t() As Double, ByRef gamma() As Double, _ ByRef gkh() As Double, ByRef gkv() As Double, ByRef alpha() As Double, ByRef matno() As Integer, _ ByRef rdis() As Double) '入力確認出力 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 = "NSTRES,NODT,NELT,MATEL,KOX,KOY,NF" : 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") 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 < KOX Then For j = 1 To KOX If i = nokx(j) Then k = 1 Next j End If If 0 < KOY Then For j = 1 To KOY If i = noky(j) Then l = 1 Next j End If 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) Next i dat = "*element characteristics" : sw.WriteLine(dat) dat = "element,node-1,node-2,node-3,E,po,t,gamma,kh,kv,alpha,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 & "," & 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) Next ne sw.Close() End Sub Private Sub OUTRESUL(ByVal fnameW As String, ByVal NODT As Integer, ByVal NELT As Integer, _ ByRef ftvec() As Double, ByRef reac() As Double, _ ByRef x() As Double, ByRef y() As Double, ByRef dis() As Double, _ ByRef strsave(,) As Double, ByRef matno() As Integer) '結果出力(Appendモード) Dim i As Integer Dim ne As Integer Dim sw As System.IO.StreamWriter Dim dat As String Dim sigx As Double Dim sigy As Double Dim tauxy As Double Dim ps1 As Double Dim ps2 As Double Dim ang As Double sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "*displacements and forces" : sw.WriteLine(dat) dat = "node,coord-x,coord-y,dir-x,dir-y,reac-x,reac-y,ftvec-x,ftvec-y" : sw.WriteLine(dat) For i = 1 To NODT 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) Next i dat = "*stresses" : sw.WriteLine(dat) dat = "element,sig-x,sig-y,tau-xy,ps1,ps2,ang,matno" : sw.WriteLine(dat) For ne = 1 To NELT '主応力計算 sigx = strsave(ne, 1) sigy = strsave(ne, 2) tauxy = strsave(ne, 3) Call PST_PL(sigx, sigy, tauxy, ps1, ps2, ang) dat = ne.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) Next ne sw.Close() End Sub Private Function CALA_PL3(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double) As Double Dim i As Integer Dim j As Integer Dim k As Integer i = kakom(ne, 1) : j = kakom(ne, 2) : k = kakom(ne, 3) CALA_PL3 = 0.5 * ((x(k) - x(j)) * y(i) + (x(i) - x(k)) * y(j) + (x(j) - x(i)) * y(k)) End Function Private Sub BMAT_PL3(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, ByRef b(,) As Double) Dim i As Integer : Dim j As Integer : Dim k As Integer Dim xi As Double : Dim xj As Double : Dim xk As Double Dim yi As Double : Dim yj As Double : Dim yk As Double Dim area As Double 'ひずみ-変位関係行列[B] area = CALA_PL3(ne, kakom, x, y) i = kakom(ne, 1) : j = kakom(ne, 2) : k = kakom(ne, 3) xi = x(i) : xj = x(j) : xk = x(k) yi = y(i) : yj = y(j) : yk = y(k) For i = 1 To 3 For j = 1 To 6 b(i, j) = 0.0 Next j, i b(1, 1) = yj - yk b(1, 3) = yk - yi b(1, 5) = yi - yj b(2, 2) = xk - xj b(2, 4) = xi - xk b(2, 6) = xj - xi b(3, 1) = b(2, 2) b(3, 2) = b(1, 1) b(3, 3) = b(2, 4) b(3, 4) = b(1, 3) b(3, 5) = b(2, 6) b(3, 6) = b(1, 5) For i = 1 To 3 For j = 1 To 6 b(i, j) = b(i, j) * 0.5 / area Next j, i End Sub Private Sub DMAT_PL(ByVal ne As Integer, ByVal NSTRES As Integer, _ ByRef Em() As Double, ByRef po() As Double, ByRef t() As Double, ByRef d(,) As Double) '応力-ひずみ関係行列[D] For i = 1 To 3 : For j = 1 To 3 : d(i, j) = 0.0 : Next j, i Select Case 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)) 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)) End Select End Sub Private Sub SM_PL3(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _ ByVal NSTRES As Integer, ByRef Em() As Double, ByRef po() As Double, ByRef t() As Double, _ ByRef sm(,) As Double) '要素剛性行列 Dim i As Integer : Dim j As Integer : Dim k As Integer Dim b(3, 6) As Double : Dim d(3, 3) As Double : Dim str(3, 6) As Double Dim area As Double area = CALA_PL3(ne, kakom, x, y) Call BMAT_PL3(ne, kakom, x, y, b) Call DMAT_PL(ne, NSTRES, Em, po, t, d) '剛性行列[sm]=A*t*[B]T[D][B] For i = 1 To 3 For j = 1 To 6 str(i, j) = 0.0 For k = 1 To 3 str(i, j) = str(i, j) + d(i, k) * b(k, j) Next k, j, i For i = 1 To 6 For j = 1 To 6 sm(i, j) = 0.0 For k = 1 To 3 sm(i, j) = sm(i, j) + b(k, i) * str(k, j) Next k, j, i If NSTRES = 0 Then t(ne) = 1.0 For i = 1 To 6 For j = 1 To 6 sm(i, j) = sm(i, j) * area * t(ne) Next j, i End Sub Private Sub TV_PL3(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _ ByRef deltaT() As Double, ByVal NSTRES As Integer, _ ByRef Em() As Double, ByRef po() As Double, ByRef t() As Double, _ ByRef alpha() As Double, ByRef fevec() As Double, ByRef strsave(,) As Double) '温度ひずみ相当ベクトル Dim i As Integer : Dim j As Integer : Dim k As Integer Dim b(3, 6) As Double : Dim d(3, 3) As Double : Dim w1(6) As Double : Dim w2(6) As Double Dim area As Double Dim Te As Double i = kakom(ne, 1) j = kakom(ne, 2) k = kakom(ne, 3) Te = (deltaT(i) + deltaT(j) + deltaT(k)) / 3.0 area = CALA_PL3(ne, kakom, x, y) Call BMAT_PL3(ne, kakom, x, y, b) Call DMAT_PL(ne, NSTRES, Em, po, t, d) '{f}=area*t*[B][D]{ε0} If NSTRES = 0 Then t(ne) = 1.0 Select Case NSTRES Case 0 '平面歪み w2(1) = (1.0 + po(ne)) * alpha(ne) * Te Case 1 '平面応力 w2(1) = alpha(ne) * Te End Select w2(2) = w2(1) w2(3) = 0.0 For i = 1 To 3 w1(i) = 0.0 For j = 1 To 3 w1(i) = w1(i) + d(i, j) * w2(j) Next j Next i '温度ひずみ相当応力 strsave(ne, 1) = w1(1) strsave(ne, 2) = w1(2) strsave(ne, 3) = w1(3) '温度ひずみ相当荷重ベクトル For i = 1 To 6 w2(i) = 0.0 For j = 1 To 3 w2(i) = w2(i) + b(j, i) * w1(j) Next j, i For i = 1 To 6 fevec(i) = w2(i) * area * t(ne) Next i End Sub Private Sub MV_PL3(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _ ByRef t() As Double, ByRef gamma() As Double, ByRef gkh() As Double, ByRef gkv() As Double, ByRef fevec() As Double) '物体力ベクトル Dim i As Integer : Dim j As Integer Dim mmat(6, 6) As Double : Dim w(6) As Double Dim area As Double area = CALA_PL3(ne, kakom, x, y) mmat(1, 1) = 0.5 : mmat(1, 2) = 0.0 : mmat(1, 3) = 0.25 : mmat(1, 4) = 0.0 : mmat(1, 5) = 0.25 : mmat(1, 6) = 0.0 mmat(2, 1) = 0.0 : mmat(2, 2) = 0.5 : mmat(2, 3) = 0.0 : mmat(2, 4) = 0.25 : mmat(2, 5) = 0.0 : mmat(2, 6) = 0.25 mmat(3, 1) = 0.25 : mmat(3, 2) = 0.0 : mmat(3, 3) = 0.5 : mmat(3, 4) = 0.0 : mmat(3, 5) = 0.25 : mmat(3, 6) = 0.0 mmat(4, 1) = 0.0 : mmat(4, 2) = 0.25 : mmat(4, 3) = 0.0 : mmat(4, 4) = 0.5 : mmat(4, 5) = 0.0 : mmat(4, 6) = 0.25 mmat(5, 1) = 0.25 : mmat(5, 2) = 0.0 : mmat(5, 3) = 0.25 : mmat(5, 4) = 0.0 : mmat(5, 5) = 0.5 : mmat(5, 6) = 0.0 mmat(6, 1) = 0.0 : mmat(6, 2) = 0.25 : mmat(6, 3) = 0.0 : mmat(6, 4) = 0.25 : mmat(6, 5) = 0.0 : mmat(6, 6) = 0.5 For i = 1 To 6 For j = 1 To 6 mmat(i, j) = mmat(i, j) * gamma(ne) * area * t(ne) / 3.0 Next j, i For i = 1 To 3 w(i * 2 - 1) = gkh(ne) w(i * 2) = gkv(ne) Next For i = 1 To 6 fevec(i) = 0.0 For j = 1 To 6 fevec(i) = fevec(i) + mmat(i, j) * w(j) Next j, i End Sub Private Sub CALST_PL3(ByVal ne As Integer, ByRef kakom(,) As Integer, _ ByRef x() As Double, ByRef y() As Double, ByRef dis() As Double, _ ByVal NSTRES As Integer, ByRef Em() As Double, ByRef po() As Double, ByRef t() As Double, _ ByRef strsave(,) As Double, ByRef fevec() As Double) '応力計算 Dim i As Integer : Dim j As Integer Dim b(3, 6) As Double : Dim d(3, 3) As Double : Dim epsilon(3) As Double : Dim area As Double Dim ia As Integer Dim wd(6) As Double : Dim ws(3) As Double '要素節点変位 For i = 1 To 3 ia = kakom(ne, i) wd(2 * i - 1) = dis(2 * ia - 1) wd(2 * i) = dis(2 * ia) Next i '節点変位による応力 If NSTRES = 0 Then t(ne) = 1.0 area = CALA_PL3(ne, kakom, x, y) Call BMAT_PL3(ne, kakom, x, y, b) Call DMAT_PL(ne, NSTRES, Em, po, t, d) For i = 1 To 3 epsilon(i) = 0.0 For j = 1 To 6 epsilon(i) = epsilon(i) + b(i, j) * wd(j) Next j, i For i = 1 To 3 ws(i) = 0.0 For j = 1 To 3 ws(i) = ws(i) + d(i, j) * epsilon(j) Next j, i '応力保存(温度応力差引き) strsave(ne, 1) = ws(1) - strsave(ne, 1) strsave(ne, 2) = ws(2) - strsave(ne, 2) strsave(ne, 3) = ws(3) - strsave(ne, 3) '内力ベクトル For i = 1 To 6 : fevec(i) = 0.0 : Next i For i = 1 To 6 For j = 1 To 3 fevec(i) = fevec(i) + b(j, i) * strsave(ne, j) * area * t(ne) Next j, i End Sub Private Sub PST_PL(ByVal sigx As Double, ByVal sigy As Double, ByVal tauxy As Double, _ ByRef ps1 As Double, ByRef ps2 As Double, ByRef ang As Double) '主応力計算 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) If sigx = sigy Then If tauxy > 0 Then ang = 45.0 If tauxy < 0 Then ang = 135.0 If tauxy = 0.0 Then ang = 0.0 Else ang = 0.5 * Math.Atan(2.0 * tauxy / (sigx - sigy)) : ang = 180.0 / 3.141592654 * ang If (sigx > sigy) And (tauxy >= 0) Then ang = ang If (sigx > sigy) And (tauxy < 0) Then ang = ang + 180.0 If sigx < sigy 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