三种有限元差分法求散度
本帖最后由 云卷云舒 于 2017-12-6 20:11 编辑各位大神，新手求助求解散度！有经向纬向两组速度数据，（附加图已给出,这是我读出数据后又整理了一下输出的，在我读取源文件是是先读的一行x）
用前差、后差、中心差求解各点的散度（偏u偏x+偏v偏y）（由于边界没法求，有的点没有散度）
但是我求出来三种方法的散度值完全不一样，不知道问题出在那里，求高手指点一下，万分感激！！！
一下是我的程序：（只看计算部分即可，print部分可以不看）
program test11 implicit none integer::n,i,j real::x(17),y(17) real,parameter::d=0.25 real::u(17,17),v(17,17) real::df_ux(16,16),df_vy(16,16),df(16,16) real::da_ux(2:17,2:17),da_vy(2:17,2:17),da(2:17,2:17) real::dm_ux(2:16,2:16),dm_vy(2:16,2:16),dm(2:16,2:16)!----------------一种方法不能只开一个数组，因为两次求导的内外层循环顺序反了！！！！ open(unit=10,file='u.txt',status='old') open(unit=20,file='v.txt',status='old') read(10,*) (x(i),i=1,17)!-----------注意这里 do i=1,17 read(20,'(18f12.7)') y(i),v(i,:) read(10,'(18f12.7)') y(i),u(i,:) end do !============================计算==================================== !=============================前差======================= !---------u/x do j=1,16 do i=1,16 df_ux(i,j)=(u(i+1,j)-u(i,j))/d!----i对应x,j对应y end do end do !---------v/y do i=1,16 do j=1,16 df_vy(i,j)=(v(i,j+1)-v(i,j))/d end do end do df=df_ux+df_vy !=============================后差======================= !---------u/x do j=2,17 do i=17,2,-1 da_ux(i,j)=(u(i-1,j)-u(i,j))/d end do end do !---------v/y do i=2,17 do j=17,2,-1 da_vy(i,j)=(v(i,j-1)-v(i,j))/d end do end do da=da_ux+da_vy !===========================中心差======================= !---------u/x do j=2,16 do i=2,16 dm_ux(i,j)=(u(i+1,j)-u(i-1,j))/(2.0*d) end do end do !---------v/y do j=2,16 do i=2,16 dm_vy(i,j)=(v(i,j+1)-v(i,j-1))/(2.0*d) end do end do dm=dm_ux+dm_vy
!==========================print=========================================== !=============================前差======================= open(unit=11,file='div_f.txt')
do i=16,1,-1 write(11,'(18f12.7)') y(i),df(i,:) end do write(11,'(12x,16f12.7)') (x(i),i=1,16) !------------------------------------- print*,"前差："
do i=16,1,-1 print'(18f10.7)', y(i),df(i,:) end do print'(10x,16f10.7)',(x(i),i=1,16) !=============================后差======================= open(unit=12,file='div_a.txt') do i=17,2,-1 write(12,'(18f12.7)')y(i),da(i,:) end do write(12,'(12x,16f12.7)')(x(i),i=2,17) !--------------------------------------- print*,"后差："
do i=17,2,-1 print'(18f10.7)', y(i),da(i,:) end do
print'(10x,16f10.7)',(x(i),i=2,17) !===========================中心差======================= open(unit=13,file='div_m.txt') do i=16,2,-1 write(13,'(18f12.7)') y(i),dm(i,:) end do write(13,'(12x,15f12.7)') (x(i),i=2,16) !------------------------------- print*,"中心差："
do i=16,2,-1 print'(18f10.7)', y(i),dm(i,:) end do print'(10x,15f10.7)', (x(i),i=2,16) end program
各种差分格式的结果有差异，这是肯定的，不然就没必要那么多格式了。
后差写错了，符号不对。
页:
[1]