蛙跳差分格式（非恒定流）

最新推荐文章于 2022-12-02 13:45:58 发布
最新推荐文章于 2022-12-02 13:45:58 发布
阅读量7k 收藏 17
点赞数 3
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_43843738/article/details/90369078
版权

1.1蛙跳格式
蛙跳格式是一个很常用的差分格式，它在时间和空间上都采用中心差分。假设离散化如下图所示：
在这里插入图片描述
在这里插入图片描述(2-1)
蛙跳格式是时间和空间上均为二阶精度的，三层一步显式差分格式。因而，为了计算(n+1)层的值，需要(n-1)层和n层的已知值。另外，还可以看出，后一偶数层的值实际上等于前一偶数层的值加上一个变化量，而跳过其间的奇数层，这就是蛙跳式这个名称的来历。
在这里插入图片描述
2.2对圣维南方程方程差分
上述格式用于圣维南方程，便得到Q和Z的显式表达式。
无侧向入流时完整偏微分方程为：
在这里插入图片描述(1-2)
按式(1-1)对(1-2)进行离散可得：
在这里插入图片描述(1-3)
上式就是由连续方程得到的在i点和j时间层的水位显式表达式。
从动量方程得：
在这里插入图片描述（1-4）
于是，对于在i点和j时间层得流量显式表达式为：
在这里插入图片描述（1-5）
阻力项也可以用下面差分形式：
在这里插入图片描述
这里Φ是权函数，取Φ=0.5则得：
在这里插入图片描述
从而得：
在这里插入图片描述
在这里插入图片描述
1.3稳定性
该格式的稳定性的必要条件：
在这里插入图片描述
考虑到阻力项的影响，∆t还必须满足：在这里插入图片描述
在这里插入图片描述
其中 在这里插入图片描述 为特征流速， 在这里插入图片描述 为重力波速， 在这里插入图片描述 为摩阻比降。

Fortran版

!b=The width of the rectangular cross-section,unit:m
!length=Total length of the stream under consideration,unit:m
!g=Gravitational acceleration,unit:m/s^2
!deltat=Time step,unit:s
!dx=Spatial step
!rn=Roughness coefficient
!t=Total computation time
!h0=Initial water depth,unit:m
!qq0=Initial discharge
!areao=Cross Section Area, Unit: m^2
!cko=Hydromodulus

program leapfrog
      use ifport
      implicit none
      integer::n,i,j,itert,t,start(2), finish(2), ierr
      real b,a,g,deltat,dx,length,lambda,&
           qq0,bb,rn,h0,Total_Time
      real,dimension(:),allocatable::z0,z,z1,zh,&
                     q0,q,q1,qh,cko,bwo,areao
      real area,ck,random1
      

      open(1,file='dataf.txt',status='old')
      open(2,file='outf.txt',status='old')
      read(1,*)b,h0,length,qq0
      read(1,*)g,deltat,dx,rn,t

      call gettimeofday(start, ierr)

      bb=length/dx
      itert=ifix(bb)
      n=itert
      write(*,*)'n=',n    
   
      allocate(z(0:n+1),q(0:n+1),z0(0:n+1),q0(0:n+1))
      !do i=1,n
       !     z(i)=h0
       !     q(i)=qq0

      !end do
      call random_seed()
      do i=1,n
           z(i)=random1(h0,h0+1.0)
           q(i)=random1(qq0,qq0+1.0)
           !write(*,*) z(i)
      end do    

      allocate(cko(n),bwo(n),areao(0:n+1))
      do i=1,n
            cko(i)=ck(b,z(i),rn)
            bwo(i)=b
            areao(i)=area(b,z(i))
      end do

      z(0) = z(n)
      z(n+1) = z(1)
      q(0) = q(n)
      q(n+1) = q(1)
      areao(0)=areao(n)     
      areao(n+1)=areao(1)

      do i=1,n+1
          z0(i) = z(i)
          q0(i) = q(i)
      end do

      lambda=deltat/dx
     

      allocate(zh(n),qh(n),z1(0:n+1),q1(0:n+1))
      !First step is FTBS
      do i=1,n
            zh(i)=z(i)-lambda*(1./bwo(i))*(q(i)-q(i-1))
      end do

      do i=1,n
            qh(i)=(1./(g*deltat*(areao(i)*abs(q(i))/cko(i)**2)+1.))&
                  *(q(i)*(1.-(g*deltat*((areao(i)*abs(q(i)))/cko(i)**2)))-2.&
                  *lambda*(q(i)/areao(i))*(q(i)-q(i-1))+lambda&
                  *((q(i)**2*bwo(i)-g*areao(i)**3)/areao(i)**2)*&
                  (z(i)-z(i-1))+lambda*(q(i)**2/areao(i)**2)*z(i)*(areao(i)-areao(i-1)))
      end do

      do i=1,n
            z1(i)=zh(i)
            q1(i)=qh(i)
      end do

      z1(0) = z1(n)
      z1(n+1) = z1(1)
      q1(0) = q1(n)
      q1(n+1) = q1(1)
      
      do i=1,n+1
            areao(i)=area(b,z1(i))
      end do

      do j=1,t/deltat
            do i=1,n
                  zh(i)=z(i)-lambda*(1./bwo(i))*(q1(i+1)-q1(i-1))
            end do

            do i=1,n
                  qh(i)=(1./(g*deltat*(areao(i)*abs(q1(i))/cko(i)**2)+1.))&
                       *(q(i)*(1.-(g*deltat*((areao(i)*abs(q1(i)))/cko(i)**2)))-&
                       2.*lambda*(q1(i)/areao(i))*(q1(i+1)-q1(i-1))+lambda&
                       *((q1(i)**2*bwo(i)-g*areao(i)**3)/areao(i)**2)&
                       *(z1(i+1)-z1(i-1))+lambda*(q1(i)**2/areao(i)**2)*z1(i)*(areao(i+1)-areao(i-1)))
            end do

            do i=0,n+1
                  z(i)=z1(i)
                  q(i)=q1(i)
            end do
            do i=1,n
                  z1(i)=zh(i)
                  q1(i)=qh(i)
            end do

            z1(0) = z1(n)
            z1(n+1) = z1(1)
            q1(0) = q1(n)
            q1(n+1) = q1(1)
      end do
      
      call gettimeofday(finish, ierr)
      Total_Time = finish(1)-start(1)+(finish(2)-start(2))/1000000.0
      write(*,*) 'Total time is',Total_time,'s'

      do i=1,n
            write(*,*) i,z0(i),z1(i),abs(z1(i)-z0(i)),q0(i),q1(i),abs(q1(i)-q0(i))
      end do
end program

function area(b,y)
      implicit none
      real b,y,area
      area=b*y
end function

function ck(b,y,rn)
      implicit none
      real ck,b,y,rn
      ck=(b*y/(2*y+b))**0.667*b*y/rn
      return
end function

function random1(lb,ub)
      implicit none
      integer i
      real lb,ub
      real random1
      real t,len

      len=ub-lb
      call random_number(t)
      random1=lb+len*t
      return
end function

C语言版

#include<stdio.h>
#include<stdlib.h>
#include<math.h>


double area(double b, double y)
{
	double a;
	a = b * y;
	return a;
}

double ck(double b, double y, double rn)
{
	double c,a,d;
	a = b * y / (2 * y + b);
	d = pow(a, 0.667);
	c = d*b*y / rn;
	return c;
}
void main()
  {
	 int i, j, itert, n, k ;
	 int m = 0;
	 double v[7];
	 double *q0, *q, *q1, *qh, *z, *z0, *z1, *zh, *cko, *areao;
	 double lambda ;
	 double delta_x, delta_t;
	 double b , length, qq0, bb, tt, h0,  t,  g = 9.81, rn=0.023;
	 FILE *fp;
	 fp = fopen("data.txt", "r");
	 if (fp == NULL)
	 	 return 0;
	 while (fscanf(fp, "%lf", &v[m]) != EOF)
		 m++;
	 fclose(fp);
	 b = v[0];
	 h0 = v[1];
	 qq0 = v[2];
	 length = v[3];
	 delta_x = v[4];
	 delta_t = v[5];
	 t = v[6];
	 bb = length / delta_x;
	 tt = t / delta_t;
	 itert = floor(bb);
	 k = round(tt);
	 n = itert;
	 printf("n=%d,k=%d\n", n, k);
	 printf("b=%f,h0=%f,qq0=%f\n", b, h0, qq0);
	 printf("length=%f,delta_x=%f,delta_t=%f,t=%f\n", length, delta_x, delta_t, t);
	 q0 = (double *)malloc((n + 1) * sizeof(double));
	 q = (double *)malloc((n + 1) * sizeof(double));
	 q1 = (double *)malloc((n + 1) * sizeof(double));
	 qh = (double *)malloc((n + 1) * sizeof(double));
	 z0 = (double *)malloc((n + 1) * sizeof(double));
	 z = (double *)malloc((n + 1) * sizeof(double));
	 z1=(double *)malloc((n + 1) * sizeof(double));
	 zh = (double *)malloc((n + 1) * sizeof(double));
	 areao = (double *)malloc((n + 1) * sizeof(double));
	 cko = (double *)malloc((n + 1) * sizeof(double));


	 for (i = 1; i <= n; i++)
	 {
		 z[i] = 1;
		 q[i] = 1;
	 }
	 for (i = 1; i <= n; i++)
	 {
		 cko[i] = ck(b, z[i], rn);
		 areao[i] = area(b, z[i]);
	 }

	 q[0] = qq0;
	q[n+1] = q[1];
	z[0] = h0;
	z[n + 1] = z[1];
	areao[0] = b * z[0];
	areao[n + 1]= areao[1];

	for (i = 1; i <= n+1; i++)
	{
		q0[i] = q[i];
		z0[i] = z[i];
	}
	
	lambda = delta_t / delta_x;
	
	//First time step is FTBS

		for (i = 1; i <= n; i++)
		{
			zh[i] = z[i] - lambda * (1 / b)*(q[i] - q[i - 1]);
			qh[i] = (1 / (g*delta_t*((areao[i] * fabs(q[i])) / (cko[i] * cko[i])) + 1))*(q[i] * (1 - (g*delta_t*((areao[i] * fabs(q[i])) / (cko[i] * cko[i])))) - 2 * lambda*(q[i] / areao[i])*(q[i] - q[i - 1]) + lambda * ((q[i] * q[i]*b - g * areao[i] * areao[i] * areao[i]) / (areao[i] * areao[i])*(z[i] - z[i - 1])) + lambda * ((q[i] * q[i]) / (areao[i] * areao[i]))*z[i] * (areao[i] - areao[i - 1]));
		}
		for (i = 1; i <= n; i++)
		{
			printf("%f,%f\n", zh[i], qh[i]);
		}

		for (i = 1; i <= n; i++)
		{
			z1[i] = zh[i];
			q1[i] = qh[i];
			
		}
		for (i = 1; i <= n; i++)
		{
			areao[i] = area(b, z1[i]);
		}
		q1[0] = q1[n];
		q1[n + 1] = q1[n];
		z1[0] = z1[n];
		z1[n + 1] = z1[n];
		areao[0]= b * z1[0];
		areao[n+1] = areao[n];

		for (j = 0; j < k; j++)
		{
			for (i = 1; i <= n; i++)
			{
				zh[i] = z[i] - lambda * (1 / b)*(q1[i + 1] - q1[i - 1]);
				qh[i] = (1/(g*delta_t*((areao[i] * fabs(q1[i])) / (cko[i] * cko[i])) + 1))*(q[i] * (1 - (g*delta_t*((areao[i] * fabs(q1[i])) / (cko[i] * cko[i])))) - 2 * lambda*(q1[i] / areao[i])*(q1[i+1] - q1[i - 1]) + lambda * ((q1[i] * q1[i] * b - g * areao[i] * areao[i] * areao[i]) / (areao[i] * areao[i])*(z1[i+1]-z1[i-1])) + lambda * ((q1[i] * q1[i]) / (areao[i] * areao[i]))*z1[i] * (areao[i+1] - areao[i - 1]));
               
			}

			for (i = 0; i <= n + 1; i++)
			{
				z[i] = z1[i];
				q[i] = q1[i];
			}
			for (i = 1; i <= n; i++)
			{
				z1[i] = zh[i];
				q1[i] = qh[i];
			}

			q1[0] = q1[n];
			q1[n + 1] = q1[n];
			z1[0] = z1[n];
			z1[n + 1] = z1[n];

			for (i = 1; i <= n; i++)
			{
				printf("%d  %f,%f,%f,%f\n", i, z0[i], z1[i], q0[i], q1[i]);
			}
		}
}

结果
在这里插入图片描述

lennonq
关注
  • 3
    点赞
  • 踩
  • 17
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
[计算流体力学][Matlab] 使用 A,B,C 格式与蛙跳格式求解二维对流问题
PriceCheap的博客
05-05 2706
1. 题目 2. 转述 原题目要求可以简化为： 对于二维对流方程： ∂u/∂t+∂u/∂x+∂u/∂y=0 u(x,y,0)={█(1,when-4≤x≤4,-4≤y≤4@0,other)┤ 取计算范围为 -16≤x≤16,-16≤y≤16，Δx=Δy=1，Δt/Δx=Δt/Δy=0.5,1,2，t>0，使用 A，B，C格式与蛙跳格式计算90个时间步长，制作 u-x-y 图随时间的动画，依此分析四种格式对于二维对流方程的稳定性。 3. 分析 之前分析一维对流方程的代码 https://blog.c
对流方程蛙跃格式matlab,对流方程各种格式代码matlab
weixin_39943220的博客
03-19 429
对流方程——偏微分方程的数值解法用迎风格式解对流方程function u = peYF(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);if a>0for j=1:(n+M)u0(j) = IniU(minx+(j-M-1)*h);endelsefor j=1:(n+M)u0(j) = IniU(minx+(j-1)*h);enden...
对称蛙跳格式的稳定性分析 (1995年)
06-12
利用待定系数法,建立了解色散方程ut=auxxx的若干蛙跳型对称差分格式.其中稳定性最好的格式稳定性条件是｜R｜≤3.2470.此外,得到一个高精度蛙跳型对称显格式的截断误差为O(τ2+h6),但其稳定性条件仅为｜R｜≤0.16208.
几种简单差分格式
热门推荐
lk040384的博客
05-05 1万+
以对流方程为例： 一阶迎风 （时间向前空间向前） （时间向前空间向后） 二阶迎风 三阶迎风 二阶中心 （无条件不稳定） 将上式的取成平均值，就得到了Lax-Friedrichs格式 在时间和空间上都用二阶中心差分离散，可以得到蛙跳（leap-frog）格式 将在进行Taylor级数展开 根据对流方程有 代入Taylor展开式得 用中...
[轻笔记]蛙跳积分法
昔风不起，唯有努力生存！
04-02 1455
知乎上JimKarrey关于蛙跳积分法的介绍潜显易懂。 现简单记录如下： 已知：积分的起点是ttt时刻的坐标x(t)x(t)x(t)，t−12δtt-\frac{1}{2}\delta tt−21​δt时刻的速度v(t−12δt)v(t-\frac{1}{2}\delta t)v(t−21​δt)，t−δtt-\delta tt−δt时刻的加速度a(t−δt)a(t-\delta t)a(t−δt)。 step1：计算出x(t)x(t)x(t)处的受力，得到ttt时刻的加速度a(t)a(t)a(t)，完
蛙跳
xusworld的专栏
04-27 439
题目9:一只青蛙一次可以跳上一级台阶，也可以跳上两级。求该青蛙跳上一个n级台阶总共有多少种情况。 #include using namespace std; //N 代表目标数字 #define N 10 class Solution { private: int * mNumArray = new int[N + 1]; int mSteps[2]; int mArrayIn
基本蛙跳程序 matlab
12-09
基本蛙跳程序 matlab/C混合蛙跳算法(SFLA)是一种基于群智能的亚启发式进 化算法，并应用于离散组合优化问题的求解。该算法是又一 种受自然界生物现象启示而产生的基于群体的协同搜索方 法。
python编程作业--盐度对流方程的差分格式设计与讨论
【更新于公众号 | 简谱学记】
03-16 2709
python盐度对流差分格式记录
一阶电路误差分析_PDE有限差分方法(12)——对流方程数值格式的分析方法
weixin_39900180的博客
10-27 803
今天直接把张强书上的第6章大部分内容完成掉. 主要讲对流方程数值格式的相容性、稳定性以及数值振荡的分析方法。参考书：(1) J.W. Thomas - Numerical Partial Differential Equations_ Finite Difference Methods (1995, Springer)(2) 张强《偏微分方程的有限差分方法》科学出版社，2019年1月版。(3) K...
波动方程HLF格式的稳定区域 (1992年)
05-23
波动方程HLF格式是指对方程ωtt=ωxx的Hamilton形式构造的显式蛙跳．(Ieap-Frog)格式，其精度为O(△t2r+△x21)，P,q为正整数。我们研究格式的强稳定性与时间方向精度的关系．结论是：对1≤p≤6格式稳定区域的变化如波形，大小交替，对p≥7，则趋于同一个区域。这种现象是扩散方程的显格式从未见到过的．
遗传算法之蛙跳算法
12-21
实现蛙跳算法，对一系列的高维复杂函数，利用蛙跳算法求其最优值
基本蛙跳算法（Frog leaping algorithm，FLA）
m0_73970831的博客
12-02 3660
蛙跳算法是基于种群进化的元启发式算法之一，通过模拟自然界中青蛙觅食过程中种群所体现出的交流与合作行为，以实现对问题的求解。
差分方程_PDE有限差分方法(11)——常系数对流方程中的常见格式
weixin_36035585的博客
01-04 1826
从今天开始进入双曲型方程的差分方法构造. 主要介绍几个常见的常系数对流方程的构造方法以及CFL方法。下一次我们来复习常系数对流方程差分格式的其他研究手段。与抛物型方程相比, 双曲线方程缺乏耗散机制, 相应数值困难更为突出.参考书：(1) J.W. Thomas - Numerical Partial Differential Equations_ Finite Difference Method...
随手笔记
weixin_51755752的博客
12-15 172
方法的重写（override/overwrite） 1.重写：子类继承父类以后，可以对父类中同名同参的方法，进行覆盖操作。 2.应用：重写以后，当创建子类对象以后，通过子类对象调用子父类中的同名同参的方法时， 实际上执行的是子类重写父类的方法。 3.重写的规定： 方法的声明：权限修饰符 返回值类型 方法名（形参列表）{ //方法体 } 约定俗成：子类中的叫重写的方法，父类中叫被重写的方法 ①子类重写的方法的方法名和形参列表与父类被重写的方法的方法名 和形参列表相同。 ②子类重写的方法的权限修饰符不小于父类被
【数值积分（欧拉积分、Velet积分与蛙跳积分）】
qq_40999248的博客
03-28 1046
数值积分（欧拉积分、Velet积分与蛙跳积分）数值积分级数展开与积分Velet积分蛙跳积分 数值积分 我们使用的各种数值积分方法都是对积分做近似，近似的本质类似于泰勒级数展开式。在实际模拟中，不可能使用泰勒级数展开中的无穷项，所以积分方法只用级数的截断类似。更健壮的积分方法是使用级数中的更多项。积分方法对模拟的精度有非常大的影响。最简单的积分有欧拉积分，其次有Velet积分（类似蛙跳积分）、蛙跳积分。理解这些概念首先需了解级数展开与积分。 级数展开与积分 泰勒级数展开 s(t+h)=s(t)+s˙(t)h+
java 蛙跳台阶问题解决算法
diu_brother的博客
03-20 1475
题目：一只青蛙可以一次跳上一级台阶，也可以一次跳上二级，求该青蛙跳上一个n级台阶总共有多少种跳法解题思路：我们把n级台阶的跳法看成是n的函数，记为f(n),当n>2时，第一次调的时候就有两种不同的选择，一是第一次跳一级，则此时的跳法数目等于后面剩下的n-1级台阶的跳法数目，即为f(n-1)；另一种选择是一次跳2级，此时跳法数目等于后面剩下的n-2级台阶跳法数目，即为f(n-2).因此n级台阶不同跳法
蛙跳算法原理
qq_43102141的博客
03-22 4460
这里写自定义目录标题蛙跳算法原理来源 蛙跳算法原理 蛙跳算法原理 来源 链接: link.
蛙跳格式matlab解双曲型方程
最新发布
06-12
蛙跳格式是一种数值解法，可以用于求解双曲型偏微分方程。以下是使用蛙跳格式求解双曲型方程的 MATLAB 代码示例： 假设要求解的方程为： ∂u/∂t + a∂u/∂x = 0 其中 a 为常数，初始条件为 u(x,0) = f(x)，边界条件为 u(0,t) = g(t)。 首先定义一些参数： ``` a = 1; % 常数 a h = 0.1; % 空间步长 k = 0.02; % 时间步长 x = 0:h:1; % 空间网格 t = 0:k:1; % 时间网格 n = length(x); % 空间网格数 m = length(t); % 时间网格数 ``` 然后初始化 u 数组： ``` u = zeros(n,m); % 初始化 u 数组 u(:,1) = f(x); % 初始条件 u(1,:) = g(t); % 边界条件 ``` 接着使用蛙跳格式求解： ``` for j = 2 : m for i = 2 : n-1 u(i,j) = u(i,j-1) - a*k/(2*h) * (u(i+1,j-1) - u(i-1,j-1)) + ... (a*k/h)^2/2 * (u(i+1,j-1) - 2*u(i,j-1) + u(i-1,j-1)); end end ``` 最后绘制结果图像： ``` [X,T] = meshgrid(t,x); surf(X,T,u'); ``` 完整代码如下： ``` a = 1; % 常数 a h = 0.1; % 空间步长 k = 0.02; % 时间步长 x = 0:h:1; % 空间网格 t = 0:k:1; % 时间网格 n = length(x); % 空间网格数 m = length(t); % 时间网格数 f = @(x) exp(-100*(x-0.5).^2); % 初始条件 g = @(t) sin(2*pi*t); % 边界条件 u = zeros(n,m); % 初始化 u 数组 u(:,1) = f(x); % 初始条件 u(1,:) = g(t); % 边界条件 for j = 2 : m for i = 2 : n-1 u(i,j) = u(i,j-1) - a*k/(2*h) * (u(i+1,j-1) - u(i-1,j-1)) + ... (a*k/h)^2/2 * (u(i+1,j-1) - 2*u(i,j-1) + u(i-1,j-1)); end end [X,T] = meshgrid(t,x); surf(X,T,u'); ``` 注意：此处仅提供了一个简单的示例，实际应用中需要根据具体问题进行调整和优化。
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值