【Fortran】TDMA算法与一维稳态导热方程
本帖最后由 watermelon 于 2020-10-4 22:05 编辑【注意】：帖子的主要部分是一个探索的过程，最终我和李新宇同学得到了正确的结果，正确的结果可以见该帖子下面的评论.....
复习数值传热学的时候，发现作业里面有一个非常有趣的编程问题，很感兴趣，跳着也要做它！
首先需要非常感谢李新宇同学在数值传热学解法和FORTRAN程序结构上的指导，正是宇哥给我指出了程序中源项在平均的情况是整体节点的平均，才使数值结果最终保持稳定而不会突变！
问题是这个样子的：
那有同学可能要问了，什么是TDMA算法呢？
TDMA算法就是使用Gauss消元法和Thomas算法来求解三对角矩阵的一个特殊方法，用迭代的方法来求解最终的值。在一维导热问题中，使用区域离散方法B对问题进行离散后，经常会得到一个三对角矩阵，这个时候TDMA算法就派上用场了。
对于本问题，有同学可能会有一个比较重要的问题：既然是求解三对角矩阵，那为什么直接用三对角矩阵的LU分解-----“追赶法”来进行求解呢？（追赶法的计算代码见：小弟的上一个帖子的评论中，评论中会有一个Fortran代码实现的追赶法）。
这个问题也是我一开始也所考虑的，但是当我对TDMA和一维导热问题进行复习之后，我们可以看到，在一开始的时候，我们是无法知道三对角矩阵中非零项元素的值的，只有在对原微分方程进行有限体积法（FVM--Finite Volume Method）进行离散，并且对问题的网格进行划分后，我们才能确定可以描述该问题的三对角矩阵中各个非零项元素的值，而这个时候，我们既可以通过三对角矩阵进行追赶法的求解，也可以直接用初始边界条件来直接进行回代求解；显然在这个问题中，我们直接进行回代求解会更加的方便（等一下就可以看出TDMA方法的代码量有多精简了）。
但是，在上面的那道作业题中，我们怎么确定其中的A,B,C,D和温度T呢？
想到这里，我就想自己实现一个简单地一维网格划分，然后利用FVM方法来离散方程，获取各个未知量前面的系数，最后使用TDMA方法来对三对角矩阵进行求解。
由于是第一次，我们把条件设置的简单一些：
1.各个节点的导热率λ（lambda）都是相同的；
2.各个节点的源强都是一个恒定的常数（在此问题中我们可以设置总源强为1）；
3.本程序只讨论稳态问题；
下面是代码实现（其中TDMA算法只有MODULE TDMA才是，其他都是划分网格与求解参数使用的程序）：
运行结果：
1.当设置问题的边界长度为10，节点数为11时，这个时候一个步长为1，我们可以看到最后一部分的运行结果中：是符合一维稳态常物性导热方程的开口向下的抛物线形式的，也就是说温度分布呈现一个左右对称的分布，是符合物理意义的：
2.当设置问题的边界长度为10，节点数为101是，这个时候一个步长为0.1，我们仍然可以看到该温度分布呈现一个左右对称的分布，同样符合物理意义：
但是，这个问题比较简单，可以得到解析解，所以，TDMA和FVM算出来的结果，真的符合解析解么？
请看下面这张图片：
从图中我们可以看到：其中蓝色曲线为我们使用TDMA+FVM计算出的解，灰色和黄色曲线为我们求得的解析解（由于该问题中少边界条件，我们无法确定解析解在Y轴的截距），但是我们仍然可以看出，计算所得结果与解析解的斜率在除了中间部位以外，其他部分的斜率差距还是非常大的。
为此，李新宇同学和我对TDMA的解法进行了分析：在该差分方程中，我们可以看到，A,B,C,D这四个最原始的差分方程中未知数的系数是和划分网格中的界面有关系的，同时，在TDMA解法中，我们可以看出，在求得A,B,C,D之后，我们就可以对方程组进行回代求解；由此，我们可以很明显的看出，依据TDMA的推导过程与差分方程系数，知道了这个问题只和划分网格的过程有关，这样的话，我们的边界条件就会不全，同时，我们离散过后的差分方程虽然最后可以求解，但是会造成计算结果与实际结果不符，虽然符合了物理意义，但是却不符合物理现象，由此可知，在λ为u常数时，差分方程系数只和划分网格的过程有关；另外，TDMA方法的边界条件不全，第一个边界条件只能确定边缘的温度梯度，而且目前并没有找到一个合适的位置引入第二个边界条件。因此，离散过后的差分方程虽然可以求解，符合了物理意义，但是计算值T的梯度与理论解偏差很大，不符合物理现象；还有一个比较重要的方面：导热系数λ，我们知道导热系数λ是温度的函数，在实际的物理情况中，λ不应该是一直不变的，所以如果想要获取更加真实的解，我们需要假定λ，并进行迭代求解真实的λ；但是这里我们的解析解也是在λ很定的情况下解出来的，所以这个方面在计算结果和解析解之间的差距影响不是很大。
这个时候，我们再来看一下之前的那一道作业题目：我理解的题目意思是：我们首先假定一组数据A,B,C（这三项是三对角矩阵中非零项元素），然后再用一个组温度T来求解出D（这里的D是最终方程组等号右边的向量），然后我们在有A,B,C,D的情况下使用TDMA方法求解出T；我们想，当初D项就是用A,B,C,T计算出来的，这次又返回去解方程组，那么接出来的T不是和之前的数值一样么？同学们你们认为的呢？
最后再次感谢李新宇同学在程序中关于源项的问题给指出的错误并且在程序和理论上给我很多其他的指导，同时我们下午一起动手离散方程并且进行深入交流，最后一起分析TDMA+FVM的过程，虽然最后还有很多不懂得要等开学了以后问老师，但是我们度过了一个非常美妙的下午！
参考资料：《数值传热学（第二版）》 陶文铨 编著西安交通大学出版社
西安交通大学数值传热学课件：http://nht.xjtu.edu.cn/NHT_Chapter_3-2_2020.pdf
欢迎大家一起和我们来讨论！请大家多多批评指正！
附件：Visual Studio 2013 + Intel Visual Fortran 2013工程文件，直接双击.sln文件即可打开编译运行。
本帖最后由 watermelon 于 2020-10-8 21:10 编辑
1.在本帖中，在程序大概174行的地方，计算b数组的时候，书上一个面积因子Ap和一个重要的物理参数aP，由于FORTRAN不区分大小写，导致我在当时初次编程的时候把aP写成了程序变量中的AP，结果，我反而在b项里面乘了AP进去，实际上b项的Ap是一个面积因子，在笛卡尔坐标系下为1，所以造成了计算结果的错误，宇哥重新检查了我的程序，非常的细心并且物理感觉很好，及时检查了这个错误，之前我们一直在源项那里兜圈子，后来才发现这个bug。
2.更加重要的是我们在程序中添加了第一类边界条件，这个条件也是非常重要的，这个条件说明温度不是网格所决定的，而是网格参数和边界温度共同组成的，是非常符合物理和数学上的意义的。
3.小弟同样优化了一些程序代码，把单独用于存储温度的数组T取消，而将温度的值放到了节点的结构体数组中，使得程序的空间可以高效利用并且整体结构变得更加整齐。
最后，此次一维稳态导热方程的数值计算和解析解符合的非常好，还是非常感谢宇哥的！
具体的代码见附件，之前的帖子主要内容属于探索的项目，幸好及时改正，希望大家多多批评指正！
设定问题长度为10，分别划分100，1000，10000个控制体来进行计算，结果如下：
100个控制体时：
1000个控制体时：
10000个控制体时：
模块化和type用的挺好
页:
[1]