自编轴对称热弹性有限元程序, 结果总为正确结果的一半

kappi86

大家好！发帖求助有限元c程序的问题，如下：
我采用Chandrupatla和Belegundu的《工程中的有限元方法》（Introduction to Finite Elements in Engineering，第三版，曾攀译）书中的轴对称分析代码AXIQUAD.C， 计算轴对称柱体的热膨胀问题，温升已知(已通过有限差分获得)。计算结果总是刚好为ANSYS结果的一半。

为了验证结果，我给定一均匀温升(10摄氏度)，半径为0.4m，热膨胀系数为1.1e-5，因此，热膨胀应该约为：

    ThermoExpan = 10 * 1.1e-5 * 0.4 = 4.4e-5

ANSYS计算结果正确，代码 AXIQUAD.C的结果总是约为2.2e-5。

程序代码我检查了下，由于我水平有限，没有发现有什么问题。所以请大家帮我看看这代码在哪儿出错了，先谢谢大家了。

AXIQUAD代码，输入与输出文件，ANSYS命令流与结果见附件。

tonnyw

It seems to me a simple question. Here what I would do if I were you.
1. Take one element and calculate the element stiffness and load vector by hand.
2. Compare your hand-calculation results with the ones from the routine-int elstif