一维的热传导方程向前差分法

热传导方程

一维问题
- 问题背景
- 问题抽象
- 离散方法
- - 1 向前差分
  - - 迭代格式
    - 初值设定
    - 取参数
    - C语言实现
    - fortran实现
    - 画图

hexo博客为啥用不了了，将就着先用这个博客吧

一维问题

问题背景

摘抄自网络 ：在距离为L的两个半无限长壁面之间有传热的流体。假设整个流场初始时刻具有温度T=T1(常数，本题中取10摄氏度)，并处于平衡状态(即：初始时刻壁面和流体的温度都为10摄氏度)。两个壁面的初始温度Tw1=Tw2=T1（本题中即为10摄氏度）。现假设在t=0时刻，右边的壁面温度突然增加到Tw2=T2并保持在T2（本题中取T2为20摄氏度），而左边的壁面温度保持在Tw1=T1（即10摄氏度）。求时间t1、t2和时间趋于无穷时壁面间流体温度的分布？

问题抽象

ut=α⋅∂2u∂x2初值：u(x,0)=u0(x)边值：u(0,t)=μ0(t),u(l,t)=μ1(t)u_t=\alpha\cdot \frac{\partial^2u}{\partial x^2} \\ 初值：u(x,0)=u_0(x)\\ 边值：u(0,t)=\mu_0(t),u(l,t)=\mu_1(t) ut=α⋅∂x2∂2u初值：u(x,0)=u0(x)边值：u(0,t)=μ0(t),u(l,t)=μ1(t)

离散方法

1 向前差分

u(xi,t+Δt)−u(xi,t)Δt=α⋅u(xi+1,t)−2u(xi,t)+u(xi−1,t)Δx2\frac{u(x_i,t+\Delta t)-u(x_i,t)}{\Delta t}=\alpha\cdot\frac{u(x_{i+1},t)-2u(x_{i},t)+u(x_{i-1},t)}{\Delta x^2} Δtu(xi,t+Δt)−u(xi,t)=α⋅Δx2u(xi+1,t)−2u(xi,t)+u(xi−1,t)

迭代格式

u(xi,t+Δt)=ru(xi+1,t)+(1−2r)u(xi,t)+u(xi−1,t)u(x_i,t+\Delta t)=ru(x_{i+1},t)+(1-2r)u(x_i,t)+u(x_{i-1},t)u(xi,t+Δt)=ru(xi+1,t)+(1−2r)u(xi,t)+u(xi−1,t)
其中r=aΔt(Δx)2r=\frac{a\Delta t}{(\Delta x)^2}r=(Δx)2aΔt

初值设定

T1=10T2=20L=0.2mt=10sT_1=10\\ T_2=20\\ L=0.2m\\ t=10sT1=10T2=20L=0.2mt=10s

取参数

dx=0.01dt=0.001dx=0.01\\ dt=0.001dx=0.01dt=0.001

C语言实现

#include <stdio.h>
#include <stdlib.h>
#include<math.h>
int main()
{float dx,dt,t,L=0.2,a[100000],b[100000],alpha=1.42857*10e-3;/*dx表示delta(x)，dt表示delta(t)，alpha为热扩散率。由于数组大小随输入值而变化，所以取得比较大*/int i,n1,n2,n=0,total=1;FILE *F;//文件指针F=fopen("data.txt","w");dx=0.01;dt=0.001;t=10;n1=(int)(L/dx);n2=(int)(t/dt);//printf("n1=%d  n2=%d\n",n1,n2);for(i=0;i<n1;i++){     /*初始赋值*/if(i<n1-1){a[i]=10;}else a[i]=20;fprintf(F,"%8.4f",a[i]);//写入t=0时刻的温度}fprintf(F,"\n");    //换行while(n<n2){    /*时间到了，停止循环*/for(i=0;i<n1;i++){/*计算*/if(i==0){b[i]=10;a[i]=b[i];}     /*表示左壁面的温度始终保持在10摄氏度*/else if(i==n1-1){b[i]=20;a[i]=b[i];}     /*表示右壁面的温度始终保持在20摄氏度*/else{b[i]=a[i]+alpha*dt*(a[i+1]-2*a[i]+a[i-1])/(dx*dx);   /*差分方程的表达式*/a[i]=b[i];}n++;fprintf(F,"%8.4f",a[i]);//printf("%8.4f",a[i]);if(total%n1==0){             /*输出data文件*///printf("\n");fprintf(F,"\n");}total++;}}fclose(F);return 0;
}

输出结果为data.txt

fortran实现

program f1
implicit noneinteger::i,n1,n2,n,total,counreal::a(100000),b(100000)real::dx,dt,t,L,alphaopen(unit=1,file='data-f.txt')  !打开文件n=1total=0dx=0.01dt=0.001t=10 !时间L=0.2   !长度alpha=1.42857*10e-3n1=int(L/dx)n2=int(t/dt)do coun=1,n1        !类似于C语言的for循环if (coun<n1) thena(coun)=10      !赋初始值elsea(coun)=20end ifend dodo while(n<n2)!时间限制do i=1,n1if (i==1) thenb(i)=10         !表示左壁面的温度始终保持在10摄氏度a(i)=b(i)else if(i==n1) then     !表示右壁面的温度始终保持在20摄氏度b(i)=20a(i)=b(i)else                    !其他情况需要用差分方程计算b(i)=a(i)+alpha*dt*(a(i+1)-2*a(i)+a(i-1))/(dx*dx)a(i)=b(i)end ifwrite(1,"(F8.4)",advance='no')a(i)!advancce=n0使得每次读写结束，位置不会自动向下移动一行if(mod(n,n1)==0) thenwrite(1,*)                      !每一行够n1个之后换行end ifn=n+1end doend dodo coun=1,n1write(*,"(F8.4)")a(coun) !控制台输出t=10s的温度分布end doclose(unit=1)                   !关闭文件
end program

数据储存在文件data-f.txt中

画图

想画随时间变化的动图，不知道怎么画