Fortran Coder

查看: 38|回复: 1

[有限元] 三种有限元差分法求散度

[复制链接]

3

帖子

2

主题

0

精华

新人

F 币
24 元
贡献
12 点
发表于 2017-12-6 20:11:29 | 显示全部楼层 |阅读模式
本帖最后由 云卷云舒 于 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


以u为例,左侧是y,最后一行是x

以u为例,左侧是y,最后一行是x
回复

使用道具 举报

211

帖子

1

主题

0

精华

宗师

F 币
1234 元
贡献
849 点
发表于 2017-12-6 21:44:28 | 显示全部楼层
各种差分格式的结果有差异,这是肯定的,不然就没必要那么多格式了。
后差写错了,符号不对。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|QQ群|Fcode

GMT+8, 2017-12-18 15:08

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表