程序结构力学期中作业

1、程序结构力学作业土木系结02 汪家继清华大学本科生考试试题专用纸考试课程 结构矩阵分析 2013年 1 月 11 日班级 结构02 姓名 汪家继 学号 2010010144 用Fortran90语言编制平面杆系结构自由振动的分析程序：必作：频率计算 (任意结构自由振动的前n阶频率，可求重频率)要求：正确可靠、能通过ELF90的编译、与求解器相接、尽可能优化求解速度。选作：振型计算 (根据完成情况给附加成绩，可以只作单根的振型)注：为了可与求解器教学版接口，若只作频率，振型的杆端位移可输出零值。交卷形式：1) 提交Fortran90的源程序；2) 提交一个程序报告。其中程序报告的内容包括：(1)

2、 简要的程序说明(2) 源程序清单(3) 试题纸(标明结点和单元编码)(4) 输入数据文件(5) 输出结果第一题：求图示结构的前10阶频率。各杆，误差限柱刚度：，梁刚度：， 第 1 页/共 2 页指定数据跨长 (m)层高 (m)165284332第二题：求图示结构的前8阶频率。各杆，误差限，杆长1.2 。LLL1234( 1 )( 2 )( 3 )第 2 页/共 2 页静力分析编程大作业班级 结02 姓名 汪家继 学号 2010010144 一、程序流程设计图如上所示，本次程序结构力学的作业中，要求我们在课程教材的计算代码的基础上补足尚缺少的代码，实现完整的超静定结构的静力分析部分的相应计算。

3、其中，涉及到的各个计算模块流程图如下所示：图一：主程序调用部分 图二：TypeDef部分进入主振型分析主程序定义节点： Use Module Numkind:定义双精度数据TypeDef:定义单元与节点单元长度、刚度、外载SetElement:利用已有信息求解单元所有信息定义单元信息（含主振型计算需要的质量、内插点数量等定义）定义节点荷载（此处无用）定义单元荷载（此处无用）数据定义部分：定义Elem,Joint,Eload,JloadEforce,EDisp数据读入与运算：读入文件，并调用Final_freq完成程序计算流程图三：矩阵变换模块 （Bandmat） 图四：FreqSolve模块调

4、用前述各个模块输入参量：only: typ_ElementNoEStifMat计算动力刚度矩阵SetmatBand计算得出各个列的始行码，并以此为依据开辟存储空间DelMatBand程序结束后释放存储空间GStifMat利用单元定位向量装配整体刚度矩阵NoFinal_J0，Final_JK利用Elem中的信息计算相应的J0，JkVarBandSolv, FreBandSolve两个完成刚度矩阵的LDLT分解的代码，后者只算D矩阵，用于算频率YesFinal_freq使用计数方法完成频率计算输出相应所需的刚度阵调用Final_Disp子过程计算最终的主振型图五：主振型计算部分 （被FreqSol

5、ve调用） Zhenjiao()对于随机产生的一个向量，利用前面算得主振型进行正交化分解调用前述各个模块NoFinal_Disp计算某一给定频率下的位移DEstif计算动力刚度矩阵的求导结果流程中需要使用esh,echFinal_Right计算迭代过程中的右端项结果最终检验迭代过程中的误差No使用插值方法给出更为精确的频率计算结果Fix, Fixend两个子过程分别为计算低于某频率固端振型和介于两个频率之间的振型数目Yes进行文件读写，得出结果附注：内容文件清单1、 SM90.f90为编程中所使用的源文件，内部含有所有代码，使用CVF以及ELF均可通过2、 文本文档：程序结构力学期末作业结构0

6、2汪家继，为一式两份，分别为pdf版本及word版本，两者内容完全一致。3、 计算数据输入文件：汪家继作业题1.inp，汪家继作业题2.inp4、 编程过程中自动生成的Debug文件夹及其workspace文件。5、 编译得到的SM90.exe可执行文件二、 编程总结：通过本次课程期末大作业的编写过程中，我主要使用Fortran语言完成数据读写、荷载动力刚度矩阵的计算与集成、利用Wittrick-Williams的计数方法进行基频的计算、利用Newton法进行主振型的计算等各个环节。并最后对于我们所求解的这两个问题实现了充分的精度。在编程过程中遇到的问题与经验总结如下1、理论推导与编程实际相联

7、系的学习方法：要想真正发挥计算机高超的计算能力，我们既需要对于基本原理的深入了解，也需要对在编程过程中充分考虑到如何通过程序实现我们的构想。理论与实际编程中细微的差距就有可能导致计算结果的完全错误或者完全无法得出收敛的结果。例如在计算刚度矩阵的过程中，由于此类动力刚度矩阵存在cosh(),sinh()之类的函数。这些函数在数学中理论上可以计算，但在计算机的实际应用中会由于计算机存储能力的不足发生数据溢出的现象。为了解决这一问题，老师在课程中提出的方法为使用esh(),ech(),这两个函数用以避免计算其中的指数暴涨的成分。在初期的编程计算中，我处于编写方便的考虑并未使用该做法，后期的实际计算显

8、示出这会使得结构在较高的主振型（100阶以上）出现较多的计算误差，甚至有可能出现软件无法运行通过的情况。为此，我通过细致的排查将刚度矩阵以及刚度矩阵的导数矩阵中的各个元素进行通分与检验，彻底避免了指数爆炸的因式，最终合理地解决了这一问题。此外，在计算主振型时，最为基础的步骤在于给出每一个杆件所对应的固端主振型的数量。在一开始的计算中，我曾经分别考虑了固端轴向振动、弯曲振动且K为偶数、弯曲振动且K为奇数这三种情形并讨论了三者的相互逻辑关系。此后在编程后期我加深了对这一现象的理解。可见，在此处充分地理解计算的原理可以让我们很好地了解到程序中的逻辑关系，避免不必要的浪费工作。最后，在编写完成相应代码

9、并生成EXE文件之后，我又将已有的计算结果与学生版本的结构力学求解器中求得的结果进行互相比对，通过这一过程，我们逐步发现了已有程序中的问题，并逐个排除。2、重频的处理方法本次大作业最为值得一提的是我们对于重频频率的计算方法。我在程序中使用的解决方法是在程序主振型求解的这一子过程Final_Disp（）中引入参量cg。这一参量为重根的拼音缩写，也就是说每当我们调用Final_Disp模块之时，我们就已经是从上级模块中读取出相应的重频次数，并依据这一结果进行计算。而在这里的计算一方面需要使用通用的牛顿算法，另一方面需要确保一开始选取的尝试位移与前面已经求解得出的位移保持互相正交的关系。这一点是通过

10、调用的一个子过程zhenjiao( )而得以解决。对于结构2，当我们将我们的结果与学生版本的结构力学求解器互相比对时，可以发现二者略有出入。例如当2号结构的2阶、3阶主振型对应着同一个频率时，我所编写的代码计算结果如下表所示。通过比对，两个表格计算结果看似不同，但实质一样。因为是重频，相当于我们所求皆的代数方程有两个根，这两个根在空间中可以张成一个 2维子空间。换言之，任何两个解的线性组合为新解。所以可以验证这两个主振型计算结果对应的结果可以互相线性表出，故本质是一样的。此外，学生版本的求解器给出的计算结果较为特殊，恰好显示了单元内部的振动情况，可见这必定是在求解的算法中有更为深度的优化。在未

11、来的学习中我也很希望能学到这里的处理方法。但是在求解工程问题中，给出的这些解答均为可以的，都可以用来进行结构的动力响应分析之中，表格1：程序结构力学求解器（自编计算代码）重频结果第 2 阶频率 = 10.7069610603096振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 2 0.00000000 0.00000000 1.00000000 0

12、.00000000 0.00000000 0.00000000 3 0.00000000 0.00000000 0.75800587 0.00000000 0.00000000 0.00000000-第 3 阶频率 = 10.7069610603096振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 2 0.00000000 0.00000000 -

13、0.75807099 0.00000000 0.00000000 0.00000000 3 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000-表格2：程序结构力学求解器（学生版）重频结果第 2 阶频率 = 10.7070872961049振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 -0.00000000 0.00000000 0.00000000 0

14、.00000000 2 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 3 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000-第 3 阶频率 = 10.7070873015263振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 -0.00000000 0.00000000

15、-0.00000000 0.00000000 2 0.00000000 -0.00000000 -0.00000000 0.00000000 0.00000000 0.00000000 3 0.00000000 -0.00000000 1.00000000 0.00000000 0.00000000 0.000000003、主振型计算中特征值的作用在计算中主振型过程，我们的迭代的每一个步骤都会记录下特征向量中最大元素的数值。但需要注意的是这里的选法：即在这一向量中先找出绝对值最大的元素，此后再取出其数值，求倒数就是相应的特征值。由于我刚开始编程时没有深入理解这里的这一理论依据，所以在选取特征值

16、时选取结果均为绝对值，这就必然导致了计算结果不收敛的情况。在深入研读课本之后最终解决了这一问题，此外，这一特征值在牛顿算法的最后还可以用来估计更为准确的频率。在实际编程中我也注意到了这一点，提高了频率计算的精确程度。4、加速算法的使用：在初期编写程序的过程中，我并未深入意识到期中作业中的加速算法的作用，而是选取了课本中的算法进行LDLT分解，通过实际的编程计算，我发现高效率的加速算法是本程序提高计算速率（收敛速率）的关键所在。通过细致地研究LDLT方法，我在期中作业中已经使用了加速的算法，在此处沿用该算法可以大大提高计算效率。这一点在计算高阶频率及主振型之中显示了重要的功效。可以使得我们很快地

17、求解出各个频率。5、基于模块调用的存储空间优化：通过本次程序设计过程，我深刻地体会到模块化设计的重要作用，在进行各个频率所对应的主振型计算中，我们往往会进行单元的细分以及相应的新结构体系的定义过程。但在这一过程中我们会产生大量的信息。如果不能很好地及时处理这些信息将会较多地占用计算机存储空间，为此我们在计算主振型的时候使用子程序Final_Disp进行计算，每次计算得出一个频率所对应的若干主振型之后就将各个生成的刚度矩阵占用的空间直接释放掉，故可以很好地避免资源浪费。三、程序设计结果文件：3.1输入文件：输入文件内容汪家继作业题1.inp：（其中长度参数与期中作业一样）9结点,1,0,0结点,

18、2,0,5结点,3,6,0结点,4,6,4结点,5,6,8结点,6,6,12结点,7,6,16结点,8,14,0结点,9,14,4结点,10,14,8结点,11,14,12结点,12,14,16结点,13,17,0结点,14,17,2单元,1,2,1,1,0,1,1,1单元,2,5,1,1,1,1,1,0单元,3,4,1,1,1,1,1,1单元,4,5,1,1,1,1,1,1单元,5,6,1,1,1,1,1,1单元,6,7,1,1,1,1,1,0单元,4,9,1,1,1,1,1,1单元,5,10,1,1,1,1,1,1单元,6,11,1,1,1,1,1,1单元,7,12,1,1,0,1,1,1

19、单元,8,9,1,1,0,1,1,1单元,9,10,1,1,1,1,1,1单元,10,11,1,1,1,1,1,1单元,11,12,1,1,1,1,1,1单元,9,14,1,1,0,1,1,1单元,13,14,1,1,1,1,1,1结点支承,1,4,0,0,0结点支承,8,4,0,0,0结点支承,13,6,0,0,0,0结点支承,3,6,0,0,0,0单元材料性质,1,16,100000,0,0,0,-1单元材料性质,1,16,10E5,160000,1,0,-1单元材料性质,2,2,3000000,62500,1,0,-1单元材料性质,7,10,3000000,62500,1,0,-1单元材

20、料性质,15,15,3000000,62500,1,0,-1自振频率参数,10,1,0.00005 图六：作业1各个单元定义图示3.2第一题结构计算结果：表格3：第一题结构主振型计算结果列表振动分析求解数目 = 10起始阶数 = 1误 差 限 = .00005第 1 阶频率 = 9.78617857968478振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 -0.07487937 0.33524056 0.00138386 -0.05

21、183902 2 0.33524056 0.00138386 -0.05183902 0.32909550 0.01481397 0.03109072 3 0.00000000 0.00000000 0.00000000 0.06339755 0.00826514 -0.03780237 4 0.06339755 0.00826514 -0.03780237 0.32909550 0.01481397 -0.07430778 5 0.32909550 0.01481397 -0.07430778 0.66571415 0.01917488 -0.07493468 6 0.66571415 0.

22、01917488 -0.07493468 1.00000000 0.02011379 -0.08720256 7 0.06339755 0.00826514 -0.03780237 0.05895755 0.00280578 -0.03670239 8 0.32909550 0.01481397 -0.07430778 0.32867973 -0.00461660 -0.07489970 9 0.66571415 0.01917488 -0.07493468 0.66660524 -0.00869006 -0.07762486 10 1.00000000 0.02011379 0.033599

23、30 0.99932169 -0.00937935 -0.07523296 11 0.00000000 0.00000000 -0.00376697 0.05895755 0.00280578 -0.03670239 12 0.05895755 0.00280578 -0.03670239 0.32867973 -0.00461660 -0.07489970 13 0.32867973 -0.00461660 -0.07489970 0.66660524 -0.00869006 -0.07762486 14 0.66660524 -0.00869006 -0.07762486 0.999321

24、69 -0.00937935 -0.07523296 15 0.05895755 0.00280578 0.00926991 0.04618604 -0.00595284 -0.03037002 16 0.00000000 0.00000000 0.00000000 0.04618604 -0.00595284 -0.03037002-第 2 阶频率 = 29.2144360829766振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000

25、0.18269664 -0.82776323 0.01048729 0.14110682 2 -0.82776323 0.01048729 0.14110682 -0.82196263 0.01850865 -0.11977042 3 0.00000000 0.00000000 0.00000000 -0.31903516 0.00656450 0.13072736 4 -0.31903516 0.00656450 0.13072736 -0.82196263 0.01850865 0.02320688 5 -0.82196263 0.01850865 0.02320688 -0.350070

26、07 0.03654376 -0.22172980 6 -0.35007007 0.03654376 -0.22172980 1.00000000 0.04501445 -0.39260455 7 -0.31903516 0.00656450 0.13072736 -0.30208041 -0.05527383 0.11328003 8 -0.82196263 0.01850865 0.02320688 -0.81291813 -0.06247923 0.03297092 9 -0.35007007 0.03654376 -0.22172980 -0.34797780 -0.06979240

27、-0.22998789 10 1.00000000 0.04501445 0.23391607 0.99812600 -0.06676205 -0.35853286 11 0.00000000 0.00000000 0.05735434 -0.30208041 -0.05527383 0.11328003 12 -0.30208041 -0.05527383 0.11328003 -0.81291813 -0.06247923 0.03297092 13 -0.81291813 -0.06247923 0.03297092 -0.34797780 -0.06979240 -0.22998789

28、 14 -0.34797780 -0.06979240 -0.22998789 0.99812600 -0.06676205 -0.35853286 15 -0.30208041 -0.05527383 -0.02809682 -0.21739017 0.02593589 0.14542057 16 0.00000000 0.00000000 0.00000000 -0.21739017 0.02593589 0.14542057-第 3 阶频率 = 46.8752273442342振型的杆段位移值 ( 乘子 = 1)- 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-

29、转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 -0.09118607 0.37860201 -0.01233055 -0.05705240 2 0.37860201 -0.01233055 -0.05705240 0.24551133 0.22987842 0.12599372 3 0.00000000 0.00000000 0.00000000 0.27802048 0.12094985 -0.06010630 4 0.27802048 0.12094985 -0.06010630 0.24551133 0.22987842 0.10997

30、477 5 0.24551133 0.22987842 0.10997477 -0.28916387 0.31086834 0.07001858 6 -0.28916387 0.31086834 0.07001858 -0.03624153 0.36763600 -0.13382729 7 0.27802048 0.12094985 -0.06010630 0.27284014 0.18178839 -0.04215179 8 0.24551133 0.22987842 0.10997477 0.23627013 0.32916031 0.07288229 9 -0.28916387 0.31

31、086834 0.07001858 -0.30113630 0.44789698 0.08767315 10 -0.03624153 0.36763600 1.00000000 -0.02340391 0.53142709 -0.32961970 11 0.00000000 0.00000000 -0.08354700 0.27284014 0.18178839 -0.04215179 12 0.27284014 0.18178839 -0.04215179 0.23627013 0.32916031 0.07288229 13 0.23627013 0.32916031 0.07288229

32、 -0.30113630 0.44789698 0.08767315 14 -0.30113630 0.44789698 0.08767315 -0.02340391 0.53142709 -0.32961970 15 0.27284014 0.18178839 -0.04340531 0.13136581 -0.01021976 -0.09575679 16 0.00000000 0.00000000 0.00000000 0.13136581 -0.01021976 -0.09575679-第 4 阶频率 = 63.8454672652027振型的杆段位移值 ( 乘子 = 1)- 杆端 1

33、 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0.00000000 0.00000000 -0.18983557 0.10780922 0.10364491 0.27827249 2 0.10780922 0.10364491 0.27827249 -0.22018283 0.80343592 -0.56500242 3 0.00000000 0.00000000 0.00000000 -0.24223934 0.43962192 0.10952373 4 -0.24223934 0.43962192 0.10952373

34、 -0.22018283 0.80343592 -0.08397042 5 -0.22018283 0.80343592 -0.08397042 0.40188408 0.99360244 0.10870929 6 0.40188408 0.99360244 0.10870929 -0.22421833 1.00000000 0.18232751 7 -0.24223934 0.43962192 0.10952373 -0.22993261 0.17960735 -0.00488373 8 -0.22018283 0.80343592 -0.08397042 -0.21879544 0.373

35、18675 -0.15054371 9 0.40188408 0.99360244 0.10870929 0.42913524 0.47594824 -0.21019826 10 -0.22421833 1.00000000 -0.92813677 -0.21773820 0.44174331 0.45967279 11 0.00000000 0.00000000 0.09313262 -0.22993261 0.17960735 -0.00488373 12 -0.22993261 0.17960735 -0.00488373 -0.21879544 0.37318675 -0.15054371 13 -0.21879544 0.37318675 -0.15054371 0.42913524 0.47594824 -0.21019826 14 0.42913524 0.47594824 -0.21019826 -0.21773820 0.44174331 0.45967279 15 -0.22993261 0.17960735 -0.1630