平面杆系结构有限元编制程序

1、平面杆系结构有限元编制程序 斜腿刚构计算一、程序编制说明 本程序全部采用Visual Basic语言编写，可以进行简单杆系结构的受力分析，能够精确计算出结构体系的结点位移和内力，并采用文本方式输入和输出数据，方便实用。例子是针对斜腿刚构的计算，全桥共划分为6个单元，计算简图如下。二、子程序划分DATA：数据的准备与输入ZG：总体刚体矩阵DG：单元刚度矩阵TRANSLATE：坐标转换矩阵RESULT：杆端位移与杆端力计算MULT：矩阵相乘MULT1：矩阵与向量相乘GAS：高斯消元法BOUNDARY：边界条件LOAD：外荷载三、标识符说明：1、整数型变量NN：结点总个数ELE：单元总个数GD：固定

2、支座个数JZ：铰支座个数ZC：支承个数MS:系统总自由度数目（总刚阶数）PN：荷载作用的结点数（结点荷载个数）2、一维数组变量EI：杆件抗弯刚度EA：杆件抗拉刚度Cn：cos值Sn：sin值LN：单元左端结点标号数组RN：单元右端结点标号数组FX：水平方向结点外力FY：垂直方向结点外力MOMENT：结点弯矩Z：总体坐标系下结点位移存储P：总体坐标系下结点力存储3、二维数组变量BKE：局部坐标单元刚度矩阵存储TKE：体系总体刚度矩阵存储T：坐标转换矩阵存储四、源程序代码Dim I, J As DoubleDim NN As Integer, ELE As Integer, MS As Integ

3、erDim GD As Integer, JZ As Integer, ZC As Integer, PN As IntegerDim LRN(1 To 12, 1 To 12) As Double, L(1 To 7) As Double, G(1 To 21, 1 To 21) As DoubleDim EA(1 To 7) As Double, EI(1 To 7) As Double, Cn As Double, Sn As DoubleDim X(1 To 7) As Double, Y(1 To 7) As Double, P(1 To 21) As DoubleDim TKE(1

4、 To 21, 1 To 21) As Double, TKE1(1 To 21, 1 To 21) As DoubleDim BKE(1 To 6, 1 To 6) As Double, T(1 To 6, 1 To 6) As DoubleDim TZ(1 To 6, 1 To 6) As Double, Z(1 To 21) As DoubleDim LN(1 To 6) As Integer, RN(1 To 6) As IntegerDim FX(1 To 7) As Double, FY(1 To 7) As Double, MOMENT(1 To 7) As DoublePriv

5、ate Sub Form_Click() Open e:homeworkgginput.txt For Input As #1 Open e:homeworkggoutput.txt For Output As #2 Call Data 输入数据 Call ZG 总刚度矩阵 Call Boundary 输入边界条件 Call load 输入荷载 Call Result 计算结果（杆端力） Close #1 Close #2End SubPublic Sub Data() 数据的准备与输入 Print #2, 平面杆系有限元程序斜腿刚构计算 Print #2, = Print #2, 原始数据的

6、输入 Print #2, Print #2, 结点数; Spc(3); 单元数; Spc(3); 固定支座数; Spc(3); 铰支座数; Spc(3); 支承数; Spc(3); 结点荷载数; Spc(3); 系统总自由度数 Input #1, NN, GD, JZ, ZC, PN MS = 3 * NN: ELE = NN - 1 Print #2, NN; Spc(6); ELE; Spc(6); GD; Spc(6); JZ; Spc(6); ZC; Spc(6); PN; Spc(6); MS Print #2, Print #2, = Print #2, 结点坐标的输入 Print

7、 #2, Print #2, 结点; Spc(3); X; Spc(6); Y; Spc(3); 单位：米 For I = 1 To NN Input #1, I, X(I), Y(I) Print #2, I; Spc(3); X(I); Spc(6); Y(I) Next I Print #2, Print #2, = Print #2, 杆件相关信息的输入 Print #2, Print #2, 杆件号; Spc(3); 左结点号; Spc(3); 右结点号; Spc(3); 抗弯刚度; Spc(3); 抗拉刚度; Spc(3); 杆件长度; Spc(3); E=2.06E+11N/m2

8、 For I = 1 To ELE Input #1, I, LN(I), RN(I), EI(I), EA(I) Next I For I = 1 To ELE L(I) = Sqr(Y(RN(I) - Y(LN(I) 2 + (X(RN(I) - X(LN(I) 2) Print #2, I; Spc(8); LN(I); Spc(8); RN(I); Spc(8); Format(EI(I), 0.000E+00); Spc(8); Format(EA(I), 0.000E+00); Spc(8); (Format(L(I), 0.000) Next I 输入结点编号 LRN(1, 1)

9、 = 3#: LRN(1, 2) = 3#: LRN(1, 3) = 4# LRN(2, 1) = 1#: LRN(2, 2) = 2#: LRN(2, 3) = 3# LRN(1, 4) = 5#: LRN(1, 5) = 7#: LRN(1, 6) = 6# LRN(2, 4) = 4#: LRN(2, 5) = 5#: LRN(2, 6) = 5# Print #2, Print #2, = Print #2, 输入结点荷载 Print #2, Print #2, 结点号; Spc(3); FX; Spc(4); FY; Spc(4); MOMENT; Spc(4); 单位：N For

10、I = 1 To NN Input #1, I, FX(I), FY(I), MOMENT(I) Next I For I = 1 To NN Print #2, I; Spc(3); Format(FX(I), 0.0E+); Spc(4); Format(FY(I), 0.0E+); Spc(4); Format(MOMENT(I), 0.0E+) Next I Print #2, Print #2, = Print #2, 输入约束条件 End SubPublic Sub DG(K) 形成单元刚度矩阵For K = 1 To ELE For I = 1 To 6 For J = 1 To

11、 6 BKE(I, J) = 0# Next J Next I BKE(1, 1) = EA(K) / L(K) BKE(2, 1) = 0# BKE(2, 2) = 12 * EI(K) / L(K) * L(K) * L(K) BKE(3, 1) = 0# BKE(3, 2) = 6 * EI(K) / L(K) * L(K) BKE(3, 3) = 4 * EI(K) / L(K) BKE(4, 1) = -EA(K) / L(K) BKE(4, 2) = 0# BKE(4, 3) = 0# BKE(4, 4) = EA(K) / L(K) BKE(5, 1) = 0# BKE(5, 2

12、) = -12 * EI(K) / L(K) * L(K) * L(K) BKE(5, 3) = -6 * EI(K) / L(K) * L(K) BKE(5, 4) = 0# BKE(5, 5) = 12 * EI(K) / L(K) * L(K) * L(K) BKE(6, 1) = 0# BKE(6, 2) = 6 * EI(K) / L(K) * L(K) BKE(6, 3) = 2 * EI(K) / L(K) BKE(6, 4) = 0# BKE(6, 5) = -6 * EI(K) / L(K) * L(K) BKE(6, 6) = 4 * EI(K) / L(K) For I

13、= 1 To 6 For J = 1 To 6 BKE(I, J) = BKE(J, I) Next J Next I Print #2, K, 单刚 For I = 1 To ELE For J = 1 To ELE Print #2, (Format(BKE(I, J), 0.000); Spc(2); Next J Print #2, Next I Next KEnd SubPublic Sub Translate(K) 形成坐标变换矩阵 For K = 1 To ELE Cn = (X(RN(K) - X(LN(K) / L(K) Sn = (Y(RN(K) - Y(LN(K) / L

14、(K) For I = 1 To 6 For J = 1 To 6 T(I, J) = 0# Next J Next I T(1, 1) = Cn T(1, 2) = Sn T(1, 3) = 0# T(2, 1) = -Sn T(2, 2) = Cn T(2, 3) = 0# T(3, 1) = 0# T(3, 2) = 0# T(3, 3) = 1# For I = 1 To 3 For J = 1 To 3 T(I + 3, J + 3) = T(I, J) Next J Next I Print #2, K; 单元, 坐标变换矩阵 For I = 1 To ELE For J = 1

15、To ELE Print #2, (Format(T(I, J), 0.000); Spc(2); Next J Print #2, Next I Next K 形成坐标变换矩阵的转置矩阵(逆矩阵） For I = 1 To 6 For J = 1 To 6 TZ(I, J) = T(J, I) Next J Next IEnd SubPublic Sub ZG() 形成总体刚度矩阵Dim M As Integer, M1 As Integer, H As Integer, H1 As IntegerDim I As Integer, J As Integer, W As Integer, I

16、I As Integer, JJ As Integer For I = 1 To MS For J = 1 To MS TKE(I, J) = 0# Next J Next I For W = 1 To ELE 形成总体坐标的单元刚度矩阵 Call Translate(W) Call DG(W) Call MULT(BKE, T, G, 6, 6, 6) Call MULT(TZ, G, TKE1, 6, 6, 6) 形成总体刚度矩阵 For I = 1 To 2 For II = 1 To 3 M = 3 * (I - 1) + II M1 = 3 * (LRN(I, W) - 1) + I

17、 For J = 1 To 2 For JJ = 1 To 3 H = 3 * (J - 1) + JJ H1 = 3 * (LRN(J, W) - 1) + JJ TKE(M1, H1) = TKE(M1, H1) + TKE1(M, H) Next JJ Next J Next II Next I Next WPrint #2, 总刚For I = 1 To MS For J = 1 To MS Print #2, (Format(TKE(I, J), 0.0E+00); Spc(2); Next J Print #2,Next IEnd SubPublic Sub MULT(A, B,

18、C, N1, N2, N3) 矩阵相乘运算Dim A(N1, N2) As Double, B(N2, N3) As Double, C(N1, N3) As Double For I = 1 To N1 For J = 1 To N3 C(I, J) = 0# Next J Next I For I = 1 To N1 For J = 1 To N3 For K = 1 To N2 C(I, J) = C(I, J) + A(I, K) * B(K, J) Next K Next J Next IEnd SubPublic Sub GAS(A, B, N, X) 高斯消元法Dim K As

19、Double, H As DoubleFor K = 1 To N - 1 For J = K + 1 To N A(K, J) = A(K, J) / A(K, K) Next J B(K) = B(K) / A(K, K) For I = K + 1 To N For J = K + 1 To N A(I, J) = A(I, J) - A(I, K) * A(K, J) Next J B(I) = B(I) - A(I, K) * B(K) Next INext KX(N) = B(N) / A(N, N)For I = N - 1 To 1 Step -1 H = 0 For J =

20、I + 1 To N H = H + A(I, J) * X(J) Next J X(I) = B(I) - HNext IEnd SubPublic Sub load() 形成结点荷载向量 For I = 1 To MS P(I) = 0# Next I P(1) = 0# P(2) = 0# P(3) = 0# P(4) = 0# P(5) = 0# P(6) = 0# P(7) = 0# P(8) = 0# P(9) = 0# P(10) = 0# P(11) = - P(12) = 0# P(13) = 35000 P(14) = 0# P(15) = 0# P(16) = 0# P(

21、17) = 0# P(18) = 0# P(19) = 0# P(20) = 0# P(21) = 0# Print #2, 结点荷载列向量 For I = 1 To MS Print #2, (Format(P(I), #.#E+) Next IEnd SubPublic Sub Boundary() 边界条件处理 For I = 4 To 5 For J = 1 To MS TKE(I, J) = 0# Next J Next I For J = 4 To 5 For I = 1 To MS TKE(I, J) = 0# Next I Next J For I = 19 To 20 For

22、 J = 1 To MS TKE(I, J) = 0# Next J Next I For J = 19 To 20 For I = 1 To MS TKE(I, J) = 0# Next I Next J For J = 1 To MS TKE(2, J) = 0# Next J For I = 1 To MS TKE(I, 2) = 0# Next I For J = 1 To MS TKE(17, J) = 0# Next J For I = 1 To MS TKE(I, 17) = 0# Next I TKE(2, 2) = 1# TKE(4, 4) = 1# TKE(5, 5) =

23、1# TKE(17, 17) = 1# TKE(19, 19) = 1# TKE(20, 20) = 1#End SubPublic Sub Result() 计算杆端位移和杆端力Dim Q(1 To 6) As Double, F(1 To 6) As Double Call GAS(TKE, P, MS, Z) Call MULT(TKE, Z, P, MS, MS) Print #2, = Print #2, 计算数据输出 Print #2, = Print #2, 杆端位移 Print #2, 结点号, 水平方向u, 垂直方向v, 转角fai Print #2, = Print #2,

24、 杆端力 Print #2, 杆件号, U（左端); Spc(4); V（左端); Spc(4); M（左端); Spc(4); U(右端); Spc(4); V（右端）; Spc(4); M（右端） For I = 1 To MS Print #2, Format(Z(I), 0.00E-00) Next I Print #2, Print #2, For I = 1 To MS Print #2, Format(P(I), 0.00000E+00) Next IEnd SubPublic Sub MULT1(A, B, C, N1, N2) 矩阵与向量相乘 For I = 1 To N1

25、C(I) = 0# Next I For I = 1 To N1 For J = 1 To N2 C(I) = C(I) + A(I, J) * B(J) Next J Next IEnd Sub五、数据的输入与结果输出数据的输入7,0,4,6,21,-0.8,3.62,0,03,4.2,3.64,11.9,3.65,23.1,3.66,39.3,3.67,38.5,01,1,3,0.618E+11,2.06E+112,2,3,0.618E+09,0.309E+113,3,4,0.618E+11,2.06E+114,4,5,0.618E+11,2.06E+115,5,7,0.618E+09,0

26、.309E+116,5,6,0.618E+11,2.06E+111,0,0,02,0,0,03,0,0,04,0,-,05,35000,0,06,0,0,07,0,0,0平面杆系有限元程序斜腿刚构计算=原始数据的输入结点数 单元数 固定支座数 铰支座数 支承数 结点荷载数 系统总自由度数 7 6 0 4 6 2 21 =结点坐标的输入结点 X Y 1 -.8 3.6 2 0 0 3 4.2 3.6 4 11.9 3.6 5 23.1 3.6 6 39.3 3.6 7 38.5 0 =杆件相关信息的输入杆件号 左结点号 右结点号 抗弯刚度 抗拉刚度 杆件长度 E=2.06E+11N/m2 1 1

27、 3 6.180E+10 2.060E+11 5.000 2 2 3 6.180E+08 3.090E+10 5.532 3 3 4 6.180E+10 2.060E+11 7.700 4 4 5 6.180E+10 2.060E+11 11.200 5 5 7 6.180E+08 3.090E+10 15.815 6 5 6 6.180E+10 2.060E+11 16.200=输入结点荷载结点号 FX FY MOMENT 单位：N 1 0.0E+0 0.0E+0 0.0E+0 2 0.0E+0 0.0E+0 0.0E+0 3 0.0E+0 0.0E+0 0.0E+0 4 0.0E+0 -1

28、.0E+5 0.0E+0 5 3.5E+4 0.0E+0 0.0E+0 6 0.0E+0 0.0E+0 0.0E+0 7 0.0E+0 0.0E+0 0.0E+0=1 单元 坐标变换矩阵 2 单元 坐标变换矩阵1.000 0.000 0.000 0.000 0.000 0.000 0.759 0.651 0.000 0.000 0.000 0.0000.000 1.000 0.000 0.000 0.000 0.000 -0.651 0.759 0.000 0.000 0.000 0.0000.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.00

29、0 0.000 0.000 0.0000.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.759 0.651 0.0000.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 -0.651 0.759 0.0000.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.000 3 单元 坐标变换矩阵 4 单元 坐标变换矩阵1.000 0.000 0.000 0.000 0.000 0.000 1.000 0

30、.000 0.000 0.000 0.000 0.0000.000 1.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.0000.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.0000.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.0000.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000

31、 1.000 0.0000.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.0005 单元 坐标变换矩阵 6 单元 坐标变换矩阵0.974 -0.228 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.0000.228 0.974 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.0000.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.0

32、00 0.000 0.000 0.0000.000 0.000 0.000 0.974 -0.228 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.228 0.974 0.000 0.000 0.000 0.000 0.000 1.000 0.0000.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.0001单元 单刚 2单元 单刚0.200 0.000 0.000 -0.200 0.000 0.000 0.027 0.000 0.0

33、00 -0.027 0.000 0.0000.000 18.000 1.800 0.000 -18.000 1.800 0.000 0.199 0.018 0.000 -0.199 0.0180.000 1.800 0.240 0.000 -1.800 0.120 0.000 0.018 0.002 0.000 -0.018 0.001-0.200 0.000 0.000 0.200 0.000 0.000 -0.027 0.000 0.000 0.027 0.000 0.000 0.000 -18.000 -1.800 0.000 18.000 -1.800 0.000 -0.199 -0.

34、018 0.000 0.199 -0.0180.000 1.800 0.120 0.000 -1.800 0.240 0.000 0.018 0.001 0.000 -0.018 0.002 3单元 单刚 4单元 单刚0.130 0.000 0.000 -0.130 0.000 0.000 0.089 0.000 0.000 -0.089 0.000 0.0000.000 27.720 1.800 0.000 -27.720 1.800 0.000 40.320 1.800 0.000 -40.320 1.8000.000 1.800 0.156 0.000 -1.800 0.078 0.00

35、0 1.800 0.107 0.000 -1.800 0.054-0.130 0.000 0.000 0.130 0.000 0.000 -0.089 0.000 0.000 0.089 0.000 0.0000.000 -27.720 -1.800 0.000 27.720 -1.800 0.000 -40.320 -1.800 0.000 40.320 -1.8000.000 1.800 0.078 0.000 -1.800 0.156 0.000 1.800 0.054 0.000 -1.800 0.1075单元 单刚 6单元 单刚0.009 0.000 0.000 -0.009 0.0

36、00 0.000 0.062 0.000 0.000 -0.062 0.000 0.0000.000 0.569 0.018 0.000 -0.569 0.018 0.000 58.320 1.800 0.000 -58.320 1.8000.000 0.018 0.001 0.000 -0.018 0.000 0.000 1.800 0.074 0.000 -1.800 0.037-0.009 0.000 0.000 0.009 0.000 0.000 -0.062 0.000 0.000 0.062 0.000 0.0000.000 -0.569 -0.018 0.000 0.569 -0

37、.018 0.000 -58.320 -1.800 0.000 58.320 -1.8000.000 0.018 0.000 0.000 -0.018 0.001 0.000 1.800 0.037 0.000 -1.800 0.074计算数据输出=杆端位移（单位；米）结点号 水平方向u 垂直方向v 转角fai1 2.70E-06 0.00E00 1.01E-062 0.00E00 0.00E00 1.54E-053 2.70E-06 -1.90E-07 8.42E-074 3.45E-06 -2.31E-07 -2.88E-075 4.53E-06 -1.35E-07 -1.66E-066

38、4.53E-06 0.00E00 -2.46E-067 0.00E00 0.00E00 2.21E-05=杆端力 （单位：N）杆件号 U（左端) V（左端) M（左端) U(右端) V（右端） M（右端）1 1.11344E+05 0.00000E+00 4.97608E+04 -5.56955E+04 -8.54568E+05 1.14291E+052 0.00000E+00 0.00000E+00 6.86927E+03 -5.56955E+04 -8.54568E+05 1.14291E+053 -5.56955E+04 -8.54568E+05 1.14291E+05 1.08313E

39、+04 -6.70673E+05 1.47854E+054 1.08313E+04 -6.70673E+05 1.47854E+05 4.99074E+04 3.87023E+05 1.16445E+055 4.99074E+04 3.87023E+05 1.16445E+05 0.00000E+00 0.00000E+00 -3.02171E+036 4.99074E+04 3.87023E+05 1.16445E+05 -3.10727E+04 0.00000E+00 5.12607E+04总体刚度矩阵4.1E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E

40、+00 -4.1E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 3.7E+12 -3.7E+11 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -3.7E+12 -3.7E+11 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E

41、+00 -3.7E+11 4.9E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 3.7E+11 2.5E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 2.1E+10 -1.8E+10 2.4E+09 -2.1E+10 1.8E+10 2.4E+09 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+

42、00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -1.8E+10 2.6E+10 -2.8E+09 1.8E+10 -2.6E+10 -2.8E+09 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 2.4E+09 -2.8E+09 4.5E+08 -2.4E+09 2.8E+09 2.2E+08 0.0E+00 0.0E+00 0.0

43、E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -4.1E+10 0.0E+00 0.0E+00 -2.1E+10 1.8E+10 -2.4E+09 8.9E+10 -1.8E+10 -2.4E+09 -2.7E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -3.7E+12 3.7E+11 1.8E+10 -2.6E+10 2.8E+09

44、 -1.8E+10 9.4E+12 2.8E+09 0.0E+00 -5.7E+12 -3.7E+11 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -3.7E+11 2.5E+10 2.4E+09 -2.8E+09 2.2E+08 -2.4E+09 2.8E+09 8.2E+10 0.0E+00 3.7E+11 1.6E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+

45、00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -2.7E+10 0.0E+00 0.0E+00 4.5E+10 0.0E+00 0.0E+00 -1.8E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -5.7E+12 3.7E+11 0.0E+00 1.4E+13 0.0E+00 0.0E+00 -8.3E+12 -3.7E+11 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 -3.7E+11 1.6E+10 0.0E+00 0.0E+00 5.4E+10 0.0E+00 3.7E+11 1.1E+10 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E