程序源代码

1 D:\计算力学大作业专区\2\2\Form1.vb Imports System.Math Public Class 平面八节点四边形有限元程序 '以下是定义一些公共变量 Public NELEM, NPOIN, NTYPE, NMATS, JAB, KGASP As Integer Public DJACB As Double Public SHAPE#(8), A#(1000000), V#(1000000), ELCOD#(2, 8), MAXA#(1000000), ESTIF#(136) Public GPCOD#(2, 9), DMATX#(3, 3), BMATX(3, 16), DERIV#(2, 8), CARTD#(2, 8) Public COORD#(2, 1000000), PROPS#(3, 1000000), LNODS#(8, 1000000), IFPRE#(1000000), MATNO# (1000000), NLOAD#(1000000) '***************************************************************************** ' 这是一个平面八节点四边形的有限元主程序。 * '***************************************************************************** Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click Dim LINECHAR, NEWLINECHAR As String Dim locatechar As Integer FileOpen(1, "DATASAVE.txt", OpenMode.Input) FileOpen(6, "FEMout.txt", OpenMode.Output) '以下进行的是从数据文件FEMDATA 中读入数据文件主标信息。 JAB = 1 LINECHAR = LineInput(1) For i = 0 To 80 Print(6, "*") Next PrintLine(6, vbCrLf + LINECHAR) For i = 0 To 80 Print(6, "*") Next LINECHAR = LineInput(1) PrintLine(6, vbCrLf + LINECHAR) '以下进行的是先读入一行字符，然后从字符中读入网格单元总数NELEM。 LINECHAR = LineInput(1) locatechar = LINECHAR.IndexOf("=") locatechar = locatechar + 2 NEWLINECHAR = Mid(LINECHAR, locatechar, Len(LINECHAR)) NELEM = Val(NEWLINECHAR) PrintLine(6, " 单元个数NELEM=" + NEWLINECHAR) ' 以下进行的是先读入一行字符，然后从字符中读入单元节点总数NPOIN。 LINECHAR = LineInput(1) locatechar = LINECHAR.IndexOf("=") locatechar = locatechar + 2 NEWLINECHAR = Mid(LINECHAR, locatechar, Len(LINECHAR)) NPOIN = Val(NEWLINECHAR) PrintLine(6, " 节点总数NPOIN=" + NEWLINECHAR) '以下进行的是先读入一行字符，然后从字符中读入问题类型编号NTYPE。 LINECHAR = LineInput(1) locatechar = LINECHAR.IndexOf("=") locatechar = locatechar + 2 NEWLINECHAR = Mid(LINECHAR, locatechar, Len(LINECHAR)) NTYPE = Val(NEWLINECHAR) PrintLine(6, " 该问题类型为（平面应力问题为1,平面应变问题为2）NTYPE=" + NEWLINECHAR) '以下进行的是先读入一行字符，然后从字符中读入材料类型总数NMATS。 LINECHAR = LineInput(1) locatechar = LINECHAR.IndexOf("=") locatechar = locatechar + 2 NEWLINECHAR = Mid(LINECHAR, locatechar, Len(LINECHAR)) NMATS = Val(NEWLINECHAR) PrintLine(6, " 材料总数NMATS=" + NEWLINECHAR) '以下进行的是调用子程序进行有限元计算。 Call INPUT(COORD, PROPS, LNODS, IFPRE, MATNO) Call STIF(COORD, PROPS, LNODS, MATNO) Call NELOAD(COORD, LNODS, NLOAD) Call ASEM(V, A, LNODS, IFPRE, NLOAD, MAXA) Call SOLVE(V, A, MAXA) 2 D:\计算力学大作业专区\2\2\Form1.vb Call STRE(COORD, PROPS, V, LNODS, MATNO) Call STRE1(COORD, PROPS, V, LNODS, MATNO) FileClose(1) FileClose(6) MsgBox("程序正常运行结束！计算结果保存在FEMout.txt文件中。") Button1.Enabled = False End Sub '***************************************************************************** ' 这是一个进行数据输入并进行数据检验的子程序。 * '***************************************************************************** Sub INPUT(ByVal COORD#(,), ByVal PROPS#(,), ByVal LNODS#(,), ByVal IFPRE#(), ByVal MATNO#()) Dim LINECHAR, NEWLINECHAR As String Dim NOTREADCHAR As Boolean Dim NEROR#(6) Dim NNODE, I, J, K, NUMBERREAD, LOCATECHAR, KELEM, IELEM, IPOIN, INODE, IGASH, IERROR, KEROR, JPOIN, KPOIN, KSTAR, KZERO As Double Dim R, angle, PI, NODST, NDOFN, MIDPT, NODMD As Double NNODE = 8 I = 1 '以下进行的是对材料类型编号赋初值。 Do While I <= NELEM MATNO(I) = 1 I = I + 1 Loop ' 以下进行的是输入材料类型编号信息标题。 LINECHAR = LineInput(1) NEWLINECHAR = LINECHAR.Replace("输入", "输出") PrintLine(6, NEWLINECHAR) ' 以下进行的是输入材料类型编号。 If NMATS > 1 Then NUMBERREAD = 1 Do While (NUMBERREAD <= NELEM) LINECHAR = LineInput(1) LOCATECHAR = LINECHAR.IndexOf(",") I = Val(Mid(LINECHAR, 1, LOCATECHAR)) MATNO(I) = Val(Mid(LINECHAR, LOCATECHAR + 2, 1)) NUMBERREAD = NUMBERREAD + 1 Loop End If '以下进行的是输出材料类型编号。 I = 1 Do While I <= NELEM PrintLine(6, " 单元编号为", I, "材料类型编号为", MATNO(I)) I = I + 1 Loop ' 以下进行的是输入单元节点总体编号标题。 LINECHAR = LineInput(1) NEWLINECHAR = LINECHAR.Replace("输入", "输出") PrintLine(6, NEWLINECHAR) ' 以下进行的是输入单元节点总体编号。 NOTREADCHAR = True Do While NOTREADCHAR LINECHAR = LineInput(1) If IsNumeric(Mid(LINECHAR, 1, 1)) Then LOCATECHAR = LINECHAR.IndexOf(",") IELEM = Val(Mid(LINECHAR, 1, LOCATECHAR)) For I = 1 To NNODE - 1 LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) LOCATECHAR = LINECHAR.IndexOf(",") LNODS(I, IELEM) = Val(Mid(LINECHAR, 1, LOCATECHAR)) Next LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) LNODS(I, IELEM) = Val(Mid(LINECHAR, 1, Len(LINECHAR))) J = 1 3 D:\计算力学大作业专区\2\2\Form1.vb Do While J <= 10 KELEM = IELEM + J If KELEM > NELEM Then GoTo 100 K = 1 Do While K <= 8 LNODS(K, KELEM) = LNODS(K, IELEM) + J K = K + 1 Loop 100: J = J + 1 Loop Else NOTREADCHAR = False End If Loop '输出单元节点总体编号 IELEM = 1 Do While IELEM <= NELEM Print(6, " 单元编号为: " + Format(IELEM, "") + " 节点总体编号为: ") For I = 1 To NNODE Print(6, LNODS(I, IELEM)) Next Print(6, vbCrLf) IELEM = IELEM + 1 Loop '以下进行的是输出单元节点的直角坐标信息标题。 NEWLINECHAR = LINECHAR.Replace("输入", "输出") PrintLine(6, NEWLINECHAR) '以下进行的是输入单元节点的直角坐标信息。 NOTREADCHAR = True Do While NOTREADCHAR LINECHAR = LineInput(1) If IsNumeric(Mid(LINECHAR, 1, 1)) Then LOCATECHAR = LINECHAR.IndexOf(",") IPOIN = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) LOCATECHAR = LINECHAR.IndexOf(",") COORD(1, IPOIN) = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) COORD(2, IPOIN) = Val(Mid(LINECHAR, 1, Len(LINECHAR))) Else NOTREADCHAR = False End If Loop ' 以下进行的是输出单元节点的极坐标信息标题。 NEWLINECHAR = LINECHAR.Replace("输入", "输出") PrintLine(6, NEWLINECHAR) '以下进行的是输入单元节点的极坐标信息。 PI = Atan(1.0) * 4.0 NOTREADCHAR = True Do While NOTREADCHAR LINECHAR = LineInput(1) If IsNumeric(Mid(LINECHAR, 1, 1)) Then LOCATECHAR = LINECHAR.IndexOf(",") IPOIN = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) LOCATECHAR = LINECHAR.IndexOf(",") R = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) angle = Val(Mid(LINECHAR, 1, Len(LINECHAR))) angle = angle / 180.0 * PI COORD(1, IPOIN) = R * Cos(angle) COORD(2, IPOIN) = R * Sin(angle) NUMBERREAD = NUMBERREAD + 1 Else 4 D:\计算力学大作业专区\2\2\Form1.vb NOTREADCHAR = False End If Loop ' 以下进行的是根据输入信息计算出单元节点的坐标。 I = 1 Do While I <= NELEM INODE = 1 Do While (INODE <= NNODE) NODST = LNODS(INODE, I) IGASH = INODE + 2 If IGASH > NNODE Then IGASH = 1 NDOFN = LNODS(IGASH, I) MIDPT = INODE + 1 NODMD = LNODS(MIDPT, I) J = 1 '以下是为了减小误差将节点坐标平均化 Do While (J <= 2) If (Abs(COORD(J, NODMD)) <= 0.000001) Then COORD(J, NODMD) = (COORD(J, NDOFN) + COORD(J, NODST)) / 2.0 End If J = J + 1 Loop INODE = INODE + 2 Loop I = I + 1 Loop ' 以下进行的是输出单元节点的坐标信息。 PrintLine(6, Space(4) + "节点编号" + Space(8) + "X坐标" + Space(8) + "Ｙ坐标") I = 1 Do While (I <= NPOIN) PrintLine(6, Space(6) + Format(I) + Space(20 - Len(CStr(I)) - Len(CStr(Format(COORD(1, I) , "0.00000")))) + Format(COORD(1, I), "0.00000") + Space(13 - Len(CStr(Format(COORD(2, I), "0. 00000")))) + Format(COORD(2, I), "0.00000")) I = I + 1 Loop ' 以下进行的是输出单元节点约束代码信息标题。 NEWLINECHAR = LINECHAR.Replace("输入", "输出") PrintLine(6, NEWLINECHAR) ' 以下进行的是输入单元节点约束代码信息。 NOTREADCHAR = True Do While NOTREADCHAR LINECHAR = LineInput(1) If IsNumeric(Mid(LINECHAR, 1, 1)) Then LOCATECHAR = LINECHAR.IndexOf(",") I = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) IFPRE(I) = Val(Mid(LINECHAR, 1, Len(LINECHAR))) Else NOTREADCHAR = False End If Loop ' 以下进行的是输出单元节点约束代码信息。 PrintLine(6, "节点约束代码信息为一个两位数的阿拉伯数字，其中") PrintLine(6, "十位数示Ｘ方向的约束状况，个位数表示Ｙ方向的约束状况。") PrintLine(6, "具体的，0 表示自由，1 表示固定，2 表示给定非零位移。") I = 1 Do While (I <= NPOIN) PrintLine(6, " 节点编号为", I, "约束代码为", Format(IFPRE(I), "00")) I = I + 1 Loop ' 以下进行的是输入材料的特性参数信息标题。 PrintLine(6, "(7)输出材料特性参数（格式为：材料类型编号NMATS,弹性模量Ｅ,泊松比μ,厚度t）") '以下进行的是输入材料的特性参数。 5 D:\计算力学大作业专区\2\2\Form1.vb NUMBERREAD = 1 Do While NUMBERREAD <= NMATS LINECHAR = LineInput(1) LOCATECHAR = LINECHAR.IndexOf(",") I = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) LOCATECHAR = LINECHAR.IndexOf(",") PROPS(1, I) = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) LOCATECHAR = LINECHAR.IndexOf(",") PROPS(2, I) = Val(Mid(LINECHAR, 1, LOCATECHAR)) LINECHAR = LINECHAR.Remove(0, LOCATECHAR + 1) PROPS(3, I) = Val(Mid(LINECHAR, 1, Len(LINECHAR))) PrintLine(6, "材料类型编号为:", I, " 的特性参数为：", Format(PROPS(0, I), "0.000e+00"), Format(PROPS(1, I), "0.000e+00"), Format(PROPS(2, I), "0.000e+00")) NUMBERREAD = NUMBERREAD + 1 Loop '以下进行的是对输入和计算出的数据进行检验。 IERROR = 1 Do While IERROR <= 6 NEROR(IERROR) = 0 IERROR = IERROR + 1 Loop ' 以下进行的是对输入和计算出的单元节点坐标进行检验。 IPOIN = 1 Do While IPOIN <= NPOIN KPOIN = IPOIN - 1 JPOIN = 1 Do While JPOIN <= KPOIN If COORD(1, IPOIN) = COORD(1, JPOIN) And COORD(2, IPOIN) = COORD(2, JPOIN) Then NEROR(1) = NEROR(1) + 1 PrintLine(6, "坐标错误！节点编号为", IPOIN, "的坐标等于", JPOIN, "的坐标") JPOIN = JPOIN + 1 Else JPOIN = JPOIN + 1 End If Loop IPOIN = IPOIN + 1 Loop ' 以下进行的是对输入和计算出的材料类型编号进行检验。 I = 1 Do While I <= NELEM If MATNO(I) <= 0 Or MATNO(I) > NMATS Then NEROR(2) = NEROR(2) + 1 End If I = I + 1 Loop If NEROR(2) <> 0 Then PrintLine(6, " 材料类型编号错误,错误个数为", NEROR(2), "个") End If ' 以下进行的是对输入和计算出的单元节点总体编号进行检验。 I = 1 Do While I <= NELEM J = 1 Do While J <= 8 If (LNODS(J, I) <= 0 Or LNODS(J, I) > NPOIN) Then NEROR(3) = NEROR(3) + 1 End If J = J + 1 Loop I = I + 1 Loop If (NEROR(3) <> 0) Then PrintLine(6, " 单元节点总体编号错误,错误个数", NEROR(3), "个") 6 D:\计算力学大作业专区\2\2\Form1.vb End If ' 以下进行的是对输入和计算出的单元节点编号进行检验。 I = 1 Do While (I <= NPOIN) KSTAR = 0 IELEM = 1 Do While (IELEM <= NELEM) KZERO = 0 J = 1 Do While (J <= 8) If (LNODS(J, IELEM) = I) Then KZERO = KZERO + 1 KSTAR = IELEM End If J = J + 1 Loop If (KZERO > 1) Then NEROR(4) = NEROR(4) + 1 PrintLine(6, "节点编号", I, " 在单元", IELEM, "中重复。") End If IELEM = IELEM + 1 Loop If (KSTAR = 0) Then NEROR(5) = NEROR(5) + 1 PrintLine(6, "总体节点编号", I, " 在单元从未出现") End If I = I + 1 Loop ' 以下进行的是对输入的单元节点约束代码进行检验。 I = 1 Do While (I <= NPOIN) If (IFPRE(I) < 0) Then NEROR(6) = NEROR(6) + 1 End If If (IFPRE(I) > 22) Then NEROR(6) = NEROR(6) + 1 End If I = I + 1 Loop If (NEROR(6) <> 0) Then PrintLine(6, "总体节点编号为", NEROR(6), "的节点约束代码出现错误。") End If ' 以下进行的是判断出错类型。 KEROR = 0 IERROR = 1 Do While IERROR <= 6 If NEROR(IERROR) <> 0 Then PrintLine(6, "出现第", IERROR, "类错误") KEROR = KEROR + 1 End If IERROR = IERROR + 1 Loop '以下进行的是判断是否出错,如出错,则中断程序运行,否则返回调用主程序。 If KEROR <> 0 Then PrintLine(6, "*****在输入数据文件中发现错误。程序运行中断*****") End '*** 程序运行被中断,这是由于输入数据出现错误引起的。***' End If End Sub '***************************************************************************** ' 这是一个形成单元刚度矩阵的子程序。 * '***************************************************************************** Sub STIF(ByVal COORD#(,), ByVal PROPS#(,), ByVal LNODS#(,), ByVal MATNO#()) Dim POSGP#(3), WEIGP#(3), ESTIF#(136), DBMAT#(3) Dim NSTRE, LPROP, LNODE, KEVAB, KGASP, IELEM As Integer 7 D:\计算力学大作业专区\2\2\Form1.vb Dim YOUNG, POISS, THICK, EXISP, ETASP, DVOLU, BTDBM As Double FileOpen(7, "ESTIF.txt", OpenMode.Output) ' 以下进行的是高斯积分点位置赋值。 POSGP(1) = -Sqrt(0.6) POSGP(2) = 0.0 POSGP(3) = Sqrt(0.6) ' 以下进行的是加权系数赋值。 WEIGP(1) = 5.0 / 9.0 WEIGP(2) = 8.0 / 9.0 WEIGP(3) = 5.0 / 9.0 ' 以下进行的是应力分量个数赋值。 NSTRE = 3 '以下进行的是形成弹性矩阵[D]。 For IELEM = 1 To NELEM PrintLine(6, "网格单元编号为：", IELEM) LPROP = MATNO(IELEM) For I = 1 To 8 LNODE = LNODS(I, IELEM) For J = 1 To 2 ELCOD(J, I) = COORD(J, LNODE) Next Next '以下是形成材料的特性参数 YOUNG = PROPS(1, LPROP) POISS = PROPS(2, LPROP) THICK = PROPS(3, LPROP) Call MODPS(YOUNG, POISS) ' 以下进行的是存放[K]的一维数组赋初值零。 KEVAB = 0 For I = 1 To 16 For J = 1 To I KEVAB = KEVAB + 1 ESTIF(KEVAB) = 0.0 Next Next '以下进行的是计算积分点处形函数及其导数之值以及对整体坐标的偏导数之值。 KGASP = 0 For IGAUS = 1 To 3 For JGAUS = 1 To 3 KGASP = KGASP + 1 EXISP = POSGP(IGAUS) ETASP = POSGP(JGAUS) PrintLine(6, " 高斯积分点编号为：", KGASP) Call SFR2(EXISP, ETASP) Call JACOB2(IELEM, DJACB, KGASP, ELCOD) DVOLU = DJACB * WEIGP(IGAUS) * WEIGP(JGAUS) * THICK Call BMATPS() ' 以下进行的是应力矩阵[D]赋初值零。 KEVAB = 0 For IEVAB = 1 To 16 For ISTRE = 1 To NSTRE DBMAT(ISTRE) = 0.0 For JSTRE = 1 To NSTRE DBMAT(ISTRE) = DBMAT(ISTRE) + BMATX(JSTRE, IEVAB) * DMATX(JSTRE, ISTRE) Next Next For JEVAB = 1 To IEVAB KEVAB = KEVAB + 1 BTDBM = 0.0 For ISTRE = 1 To NSTRE BTDBM = BTDBM + DBMAT(ISTRE) * BMATX(ISTRE, JEVAB) Next ESTIF(KEVAB) = ESTIF(KEVAB) + BTDBM * DVOLU 8 D:\计算力学大作业专区\2\2\Form1.vb Next Next Next Next '以下进行的是将单刚矩阵写入文件中 For KEVAB = 1 To 136 PrintLine(7, ESTIF(KEVAB)) Next Next FileClose(7) End Sub '***************************************************************************** ' 这是一个形成弹性矩阵的子程序。 * '***************************************************************************** Sub MODPS(ByRef Y#, ByRef P#) Dim I, J As Integer '以下是给弹性矩阵赋初值零 For I = 1 To 3 For J = 1 To 3 DMATX(I, J) = 0.0 Next Next '当为平面应力问题时，给出弹性矩阵 If NTYPE = 1 Then DMATX(1, 1) = Y / (1.0 - P * P) DMATX(2, 2) = Y / (1.0 - P * P) DMATX(1, 2) = Y / (1.0 - P * P) * P DMATX(2, 1) = Y / (1.0 - P * P) * P DMATX(3, 3) = Y / (1.0 - P * P) * (1.0 - P) / 2.0 End If '当为平面应变问题时，给出弹性矩阵 If NTYPE = 2 Then DMATX(1, 1) = Y * (1.0 - P) / ((1.0 + P) * (1 - 2.0 * P)) DMATX(2, 2) = Y * (1.0 - P) / ((1.0 + P) * (1 - 2.0 * P)) DMATX(1, 2) = Y * (1.0 - P) / ((1.0 + P) * (1 - 2.0 * P)) * P / (1.0 - P) DMATX(2, 1) = Y * (1.0 - P) / ((1.0 + P) * (1 - 2.0 * P)) * P / (1.0 - P) DMATX(3, 3) = Y / (2.0 * (1.0 + P)) End If End Sub '***************************************************************************** ' 这是一个计算形函数当前积分点及形函数对局部坐标的导数值的子程序。 * '***************************************************************************** Sub SFR2(ByRef S#, ByRef T#) Dim S2, T2, SS, TT, ST, STT, SST, ST2 As Double ' 以下进行的是为简化表达式定义一些变量。 S2 = S * 2.0 T2 = T * 2.0 SS = S * S TT = T * T ST = S * T STT = S * T * T SST = S * S * T ST2 = S * T * 2.0 '以下进行的是给形函数赋值。 SHAPE(1) = (-1.0 + ST + SS + TT - SST - STT) / 4.0 SHAPE(2) = (1.0 - T - SS + SST) / 2.0 SHAPE(3) = (-1.0 - ST + SS + TT - SST + STT) / 4.0 SHAPE(4) = (1.0 + S - TT - STT) / 2.0 SHAPE(5) = (-1.0 + ST + SS + TT + SST + STT) / 4.0 SHAPE(6) = (1.0 + T - SS - SST) / 2.0 SHAPE(7) = (-1.0 - ST + SS + TT + SST - STT) / 4.0 SHAPE(8) = (1.0 - S - TT + STT) / 2.0 ' 以下进行的是给形函数对局部坐标的导数赋值。 DERIV(1, 1) = (T + S2 - ST2 - TT) / 4.0 9 D:\计算力学大作业专区\2\2\Form1.vb DERIV(1, 2) = -S + ST DERIV(1, 3) = (-T + S2 - ST2 + TT) / 4.0 DERIV(1, 4) = (1.0 - TT) / 2.0 DERIV(1, 5) = (T + S2 + ST2 + TT) / 4.0 DERIV(1, 6) = -S - ST DERIV(1, 7) = (-T + S2 + ST2 - TT