流函数-涡量法的二维方腔流数值模拟(matlab编程)
流函数和势函数公式
流函数和势函数公式流函数与势函数是描述流体运动的两个重要概念，在流体力学中被广泛应用。
本文将介绍流函数和势函数的基本概念、性质以及求解方法。
1.流函数的概念和性质流函数是描述在二维定常流动中，各个流线上速度矢量的旋转情况的函数。
对于二维流动，假设流体流动的速度场为V(x,y)，则流函数Ψ(x,y)定义为：V=∇Ψ=(∂Ψ/∂x,∂Ψ/∂y)其中，∇Ψ是流函数Ψ的梯度向量。
流函数的性质如下：1)斜率定理：沿着流线的方向，流函数的局部斜率等于流体的速度分量。
2)流线定理：流线上的流函数值保持不变，即Ψ为常数。
3)流函数的连续性：在空间中的流函数是连续的，除非在相应的流体内有边界。
4)流函数的耗散性：流函数对时间是线性的，即流函数在时间方向上是耗散的。
2.势函数的概念和性质势函数是描述流体在无旋力场中流动时所具备的性质的函数。
无旋力场是指速度场的旋度等于零。
对于二维流动，假设流体流动的速度场为V(x,y)，则势函数φ(x,y)定义为：V=∇φ=(∂φ/∂x,∂φ/∂y)其中，∇φ是势函数φ的梯度向量。
势函数的性质如下：1)势函数的梯度向量是速度向量。
2)势流是不可压缩的，即∇·V=0。
3)势函数满足拉普拉斯方程，即∇²φ=0。
4)由于速度场的旋度等于零，势函数是无旋的。
3.流函数和势函数的关系在二维流动中，流函数和势函数之间存在一种特殊的关系，称为流函数-势函数耦合关系。
根据流函数和势函数的定义，可以得到流函数和势函数的关系：Ψ = ∫(∂φ/∂y)dx + f(y)φ = ∫(∂Ψ/∂x)dy + g(x)其中，f(y)和g(x)是任意常数函数。
根据流函数-势函数耦合关系可以求解流体的速度场，并且满足连续性方程和运动方程。
4.求解流函数和势函数的方法求解流函数和势函数的方法有多种，常用的方法有分离变量法、解析法和数值法。
4.1分离变量法分离变量法是将流函数和势函数分解为各自的变量函数，并通过解偏微分方程的边值问题来确定这些变量函数。
浓度场数值模拟matlab
浓度场数值模拟matlab1. 介绍浓度场数值模拟是一种通过数值方法来计算和预测物质浓度分布的技术。
这种模拟方法通常基于流体力学和传质过程的基本方程，通过数值求解这些方程来得到浓度场的分布。
在很多领域中，如环境科学、化学工程、材料科学等，浓度场数值模拟都是一种重要的工具。
在本文中，我们将介绍如何使用MATLAB进行浓度场数值模拟。
我们将首先介绍数值模拟的基本原理和方法，然后详细讨论如何在MATLAB中实现这些方法，并给出一些实际应用的示例。
2. 数值模拟方法数值模拟浓度场的方法有很多种，其中常用的方法包括有限差分法（finite difference method）、有限元法（finite element method）和有限体积法（finite volume method）等。
这些方法都是将浓度场分割成网格，然后通过求解离散化的方程来得到各个网格点上的浓度值。
在MATLAB中，我们可以使用这些方法的内置函数或者自己编写代码来实现浓度场的数值模拟。
下面我们将以有限差分法为例，介绍如何在MATLAB中实现浓度场的数值模拟。
有限差分法是一种将微分方程离散化的方法，它将连续的空间域划分为离散的网格点，并在每个网格点上计算浓度的近似值。
有限差分法的基本思想是使用差分近似来代替微分运算，从而将微分方程转化为代数方程。
在浓度场数值模拟中，有限差分法通常用于求解扩散方程。
3. 在MATLAB中实现数值模拟在MATLAB中，我们可以使用内置函数pdepe来求解偏微分方程。
这个函数可以用于求解一维、二维和三维的偏微分方程，并支持不同的边界条件和初值条件。
下面是一个使用pdepe函数求解一维扩散方程的示例：function [c,f,s] = diffusion_eqn(x,t,u,DuDx)c = 1;f = DuDx;s = 0;function diffusion_simulation()x = linspace(0,1,100);t = linspace(0,1,100);m = 0;sol = pdepe(m,@diffusion_eqn,@diffusion_ic,@diffusion_bc,x,t);u = sol(:,:,1);surf(x,t,u)在这个示例中，我们定义了一个名为diffusion_eqn的函数来描述扩散方程。
二维扩散方程的matlab程序
二维扩散方程是描述在二维空间中物质或能量传播的数学模型。
在科学和工程领域中，二维扩散方程被广泛应用于描述热传导、物质扩散等现象。
Matlab作为一种强大的科学计算软件，在求解和可视化二维扩散方程方面具有很大的优势。
本文将介绍如何使用Matlab编写二维扩散方程的求解程序，并通过实例演示其应用。
一、二维扩散方程模型二维扩散方程可以用以下偏微分方程表示：∂u/∂t = D(∂^2u/∂x^2 + ∂^2u/∂y^2)其中，u(x, y, t)是描述扩散物质浓度或能量分布的函数，D是扩散系数，x和y分别是空间坐标，t是时间。
上式描述了u随时间和空间坐标的变化规律，求解这个偏微分方程即可得到扩散过程中u的分布情况。
二、二维扩散方程的差分格式为了在计算机上求解二维扩散方程，我们需要将其离散化。
常用的方法是采用有限差分法，即将空间和时间分割成若干个小区间，然后在每个小区间上使用近似的差分格式来表示偏微分方程。
对于二维扩散方程，我们可以使用以下的差分格式：u(i, j, t+Δt) = u(i, j, t) + DΔt/Δx^2 * (u(i+1, j, t) - 2u(i, j, t) + u(i-1, j, t)) + DΔt/Δy^2 * (u(i, j+1, t) - 2u(i, j, t) + u(i, j-1, t))这个差分格式将时间t+Δt的u(i, j)用t时刻的值和邻近点的值表示出来，通过迭代求解，即可得到u随时间和空间的变化。
三、 Matlab程序设计在Matlab中，我们可以很方便地编写二维扩散方程的求解程序。
我们需要定义计算区域的空间网格和时间步长：```matlabLx = 1; 区域长度Ly = 1; 区域宽度Nx = 100; 空间网格数Ny = 100;dx = Lx/Nx; 空间步长dy = Ly/Ny;D = 0.1; 扩散系数dt = 0.01; 时间步长T = 1; 总的模拟时间```接下来，我们初始化u在空间上的分布，并使用差分格式进行迭代计算：```matlabu_init = zeros(Nx, Ny); 初始化uu_init(Nx/2, Ny/2) = 1; 在中心点加入扩散物质u = u_init;for t = 0 : dt : Tu_new = u;for i = 2 : Nx-1for j = 2 : Ny-1u_new(i, j) = u(i, j) + D * dt/dx^2 * (u(i+1, j) - 2*u(i, j) + u(i-1, j)) + D * dt/dy^2 * (u(i, j+1) - 2*u(i, j) + u(i, j-1));endendu = u_new;end```我们可以使用Matlab的绘图功能将扩散物质的分布进行可视化：```matlab[X, Y] = meshgrid(1:Nx, 1:Ny);surf(X, Y, u);xlabel('x');ylabel('y');zlabel('u');```四、实例演示接下来，我们通过一个具体的例子来演示上述程序的应用。
《计算流体力学》3-涡量流函数形式
§6. 平面不可压粘性流动的差分方法－涡量/流函数形式对于平面二维不可压粘性流，可用涡量 ζ 和流函数 y 来描述，称为控制方程的涡量/流函数形式22222222u v t x y x y x y ζζζζζνζì骣ï抖抖 ÷ïç÷++=+çï÷çï÷ç抖 抖桫ïïïíïï抖ï+=ïï抖ïïîy y 其中，流函数 y 和涡量 ζ 的定义分别为v x y ¶=-¶ ， u y y ¶=¶ ， v ux y抖=-抖z 方程组中的第一个方程是涡量的非线性二维对流/扩散方程，可以将以前介绍过的，对模型方程的差分近似、对非线性问题的冻结系数法、二维问题的ADI 方法等组合起来，为这一方程构造差分格式111222111222,,1,1,,1,1,,1,,1,,1,,122122222n n n n n n j k j kj k j kj k j k nnj kj kn n n n n n j k j k j k j k j k j k u v xyt x y ζζζζζζζζζζζζν++++-+-++++-+----++D D D 骣÷-+-+ç÷ç÷=+ç÷ç÷D D ç÷ç桫111222111222111,,1,1,,1,1,,1111,,1,,1,,122122222n n n n n n j k j kj k j kj k j k nnj kj kn n n n n n j k j k j k j k j k j k u v xyt x y ζζζζζζζζζζζζν+++++++-+-+++++++-+----++D D D 骣÷-+-+ç÷ç÷=+ç÷ç÷D D ç÷ç桫式中采用了ADI 方法中的P-R 格式，冻结系数 ,n j k u 和 ,nj kv 的计算需根据流函数的定义做中心差分近似,1,1,2n n j k j k n j ku y+--=D y y ， 1,1,,2n n j k j kn j kv x+--=-D y y方程组中的第二个方程是流函数的Poisson 方程，其差分格式为1111111,,1,,1,,11,2222n n n n n n j k j k j kj k j k j k n j k x y ζ+++++++-+-+-+-++=D D y y y y y y利用已经求出的 1,n j k ζ+ ，就可以从中解出 1,n j k +y （用迭代法和时间相关法）。
matlab算λci准则 涡识别准则
MATLAB算λCI准则涡识别准则一、概述在流体力学领域，涡识别是一项重要的研究内容。
涡是流体中旋转运动的局部区域，它在空气动力学、海洋工程、土木工程等许多领域都有着重要的应用。
而MATLAB算λCI准则是一种常用的涡识别方法，本文将对其进行介绍和分析。
二、MATLAB算λCI准则的原理在MATLAB中，算λCI准则采用的是基于特征值分解的方法。
该方法首先通过数值模拟或实验测量得到流体速度场的矢量场数据，然后利用特征值分解得到速度梯度张量，并将其代入λCI准则的计算公式中。
λCI准则是一种用来判断涡旋的强度和方向的指标，其计算公式如下：λCI = 0.5 * (λ2 + λ1) / (λ2 - λ1)其中，λ1和λ2分别为速度梯度张量的两个特征值，根据λCI的数值大小可以判断出该点是否存在涡旋以及其强度和方向。
三、MATLAB算λCI准则的实现在MATLAB中，算λCI准则的实现通常需要经过以下步骤：1. 读取流场数据：首先需要从数值模拟或实验测量中得到流场的速度场数据，通常以二维或三维矢量场的形式存储。
2. 特征值分解：利用MATLAB提供的特征值分解函数，将速度梯度张量转化为特征值和特征向量。
3. 计算λCI值：根据上文提到的λCI计算公式，将特征值代入计算得到每个点的λCI值。
4. 涡识别：根据λCI值的大小和正负判断出涡旋的存在和方向，以此将涡旋的位置标出或者进行进一步的分析。
四、MATLAB算λCI准则的优缺点分析1. 优点：（1）算法简单：λCI准则的计算方法相对简单直观，便于理解和实现。
（2）对流场特征敏感：λCI能够有效地识别出流场中的涡旋结构，对边界层等复杂流动情况也有较好的适应性。
2. 缺点：（1）计算结果受网格分辨率影响：当流场数据的网格分辨率较低时，λCI准则的计算结果可能会受到较大的误差影响。
（2）无法区分不同类型的涡旋：λCI准则只能判断出流场中是否存在涡旋和其强度方向，但无法明确区分出不同类型的涡旋结构。
以SIMPLE算法求解瞬态二维顶板驱动方腔流动问题
P旳
(lb)
也+也 =0
dx dy
(lc)
式(la) ~ (lc)中,t一时间,s；p—流体压力,Pa；
u、v一 ％,y方向上的速度,m/s；a—运动粘性系数,
m2/s.
2问题描述
如图1所示的方腔内充满粘性不可压缩流体, 运动粘性系数U在方腔内为常数，其左、右、下壁固 定,上壁以U = S运动.求在这样一种驱动力作用 下，方腔内的流体动力特性•在图1中，方腔的高度
1基本控制方程
对于二维顶板驱动方腔不可压均质粘性瞬态
流动问题,其流场的流体控制方程由式(la) ~(lc)
给出•
=£(u +il du + d( uu) +。(他)
dt dx
du\
丄型
~dx]
i dy) p dx
(la)
dv dt
uv) dx
1
d(仞)
=— d / V dx\
dx]
+斗 v^\ dy\ dyl
(9)
5方程求解
对于本算例的瞬态问题，在每一时间步采用
SIMPLE算法来求解时，需要计算在t + △力时刻的:
1）速度分量叫八％的初算值；2）皱八 %•的修正值此广、巧：阳；3）压力刃J的修正值
刃广•在交错网格的编码系统中，节点GJ）和（/J）
点处的速度分量初算值＜/+1、吩小和压力修正方
程的隐式格式离散方程可表达为：
\-VI,N + 2 -^VI,N+1 ~VI,N
(1 = 1,...⑻
4)压力和压力修正量物理边界条件的插值
为配合采用SIMPLE算法求解，需要将压力P
和压力修正值P'扩充到相应的虚拟网格节点中去. 对于本算例，由于所有垂宜于固壁的速度分量均为
matlab的streamslice用法
matlab的streamslice用法MATLAB的streamslice用法简介streamslice是MATLAB中用于绘制流线图的函数。
流线图是一种显示矢量场方向和强度的图形表示方法，常用于可视化流体动力学和电磁场分析等领域。
以下是streamslice函数的一些常见用法。
输入参数streamslice函数通常需要以下参数：•x和y：定义网格的坐标点矩阵。
•u和v：二维矢量场的速度分量。
•startx、starty和startz：流线起点的坐标。
•spacing：流线间隔的距离。
绘制基本流线图要绘制基本的流线图，可以使用以下代码：[x, y] = meshgrid(linspace(-2, 2, 20), linspace(-2, 2, 20));u = -y;v = x;streamslice(x, y, u, v);axis equal;上述代码首先定义了一个网格，然后计算了矢量场的速度分量。
最后通过调用streamslice函数绘制了流线图，并使用axis equal 设置了坐标轴的比例。
自定义流线起点和间距通过指定startx、starty和spacing参数，可以自定义流线的起点和间距。
以下示例代码展示了如何使用这些参数：[x, y] = meshgrid(linspace(-2, 2, 20), linspace(-2, 2, 20));u = -y;v = x;startx = linspace(-2, 2, 10);starty = linspace(-2, 2, 10);spacing = [, , , ];streamslice(x, y, u, v, startx, starty, spacing);axis equal;上述代码中，startx和starty定义了流线起点的坐标，spacing定义了流线的间距。
控制流线的可见性和密度streamslice函数提供了一些参数，可以控制流线的可见性和密度：•arrows: 指定是否在流线上绘制箭头，默认为true。
Matlab函数总结
Matlab函数大全.txt7温暖是飘飘洒洒的春雨；温暖是写在脸上的笑影；温暖是义无反顾的响应；温暖是一丝不苟的配合。
8尊重是一缕春风，一泓清泉，一颗给人温暖的舒心丸，一剂催人奋进的强心剂Matlab函数大全信源函数randerr 产生比特误差样本randint 产生均匀分布的随机整数矩阵randsrc 根据给定的数字表产生随机矩阵wgn 产生高斯白噪声信号分析函数biterr 计算比特误差数和比特误差率eyediagram 绘制眼图scatterplot 绘制分布图symerr 计算符号误差数和符号误差率信源编码compand mu律/A律压缩/扩张dpcmdeco DPCM（差分脉冲编码调制）解码dpcmenco DPCM编码dpcmopt 优化DPCM参数lloyds Lloyd法则优化量化器参数quantiz 给出量化后的级和输出值误差控制编码bchpoly 给出二进制BCH码的性能参数和产生多项式convenc 产生卷积码cyclgen 产生循环码的奇偶校验阵和生成矩阵cyclpoly 产生循环码的生成多项式decode 分组码解码器encode 分组码编码器gen2par 将奇偶校验阵和生成矩阵互相转换gfweight 计算线性分组码的最小距离hammgen 产生汉明码的奇偶校验阵和生成矩阵rsdecof 对Reed-Solomon编码的ASCII文件解码rsencof 用Reed-Solomon码对ASCII文件编码rspoly 给出Reed-Solomon码的生成多项式syndtable 产生伴随解码表vitdec 用Viterbi法则解卷积码（误差控制编码的低级函数）bchdeco BCH解码器bchenco BCH编码器rsdeco Reed-Solomon解码器rsdecode 用指数形式进行Reed-Solomon解码rsenco Reed-Solomon编码器rsencode 用指数形式进行Reed-Solomon编码调制与解调ademod 模拟通带解调器ademodce 模拟基带解调器amod 模拟通带调制器amodce 模拟基带调制器apkconst 绘制圆形的复合ASK-PSK星座图ddemod 数字通带解调器ddemodce 数字基带解调器demodmap 解调后的模拟信号星座图反映射到数字信号dmod 数字通带调制器dmodce 数字基带调制器modmap 把数字信号映射到模拟信号星座图（以供调制）qaskdeco 从方形的QASK星座图反映射到数字信号qaskenco 把数字信号映射到方形的QASK星座图专用滤波器hank2sys 把一个Hankel矩阵转换成一个线性系统模型hilbiir 设计一个希尔伯特变换IIR滤波器rcosflt 升余弦滤波器rcosine 设计一个升余弦滤波器（专用滤波器的低级函数）rcosfir 设计一个升余弦FIR滤波器rcosiir 设计一个升余弦IIR滤波器信道函数awgn 添加高斯白噪声伽罗域计算gfadd 伽罗域上的多项式加法gfconv 伽罗域上的多项式乘法gfcosets 生成伽罗域的分圆陪集gfdeconv 伽罗域上的多项式除法gfdiv 伽罗域上的元素除法gffilter 在质伽罗域上用多项式过滤数据gflineq 在至伽罗域上求Ax=b的一个特解gfminpol 求伽罗域上元素的最小多项式gfmul 伽罗域上的元素乘法gfplus GF（2^m）上的元素加法gfpretty 以通常方式显示多项式gfprimck 检测多项式是否是基本多项式gfprimdf 给出伽罗域的MATLAB默认的基本多项式gfprimfd 给出伽罗域的基本多项式gfrank 伽罗域上矩阵求秩gfrepcov GF（2）上多项式的表达方式转换gfroots 质伽罗域上的多项式求根gfsub 伽罗域上的多项式减法gftrunc 使多项式的表达最简化gftuple 简化或转换伽罗域上元素的形式工具函数bi2de 把二进制向量转换成十进制数de2bi 把十进制数转换成二进制向量erf 误差函数erfc 余误差函数istrellis 检测输入是否MATLAB的trellis结构（structure）marcumq 通用Marcum Q 函数oct2dec 八进制数转十进制数poly2trellis 把卷积码多项式转换成MATLAB的trellis描述vec2mat 把向量转换成矩阵——————————————————————————————————————————————————A aabs 绝对值、模、字符的ASCII码值acos 反余弦acosh 反双曲余弦acot 反余切acoth 反双曲余切acsc 反余割acsch 反双曲余割align 启动图形对象几何位置排列工具all 所有元素非零为真angle 相角ans 表达式计算结果的缺省变量名any 所有元素非全零为真area 面域图argnames 函数M文件宗量名asec 反正割asech 反双曲正割asin 反正弦asinh 反双曲正弦assignin 向变量赋值atan 反正切atan2 四象限反正切atanh 反双曲正切autumn 红黄调秋色图阵axes 创建轴对象的低层指令axis 控制轴刻度和风格的高层指令B bbar 二维直方图bar3 三维直方图bar3h 三维水平直方图barh 二维水平直方图base2dec X进制转换为十进制bin2dec 二进制转换为十进制blanks 创建空格串bone 蓝色调黑白色图阵box 框状坐标轴break while 或for 环中断指令brighten 亮度控制C ccapture （3版以前）捕获当前图形cart2pol 直角坐标变为极或柱坐标cart2sph 直角坐标变为球坐标cat 串接成高维数组caxis 色标尺刻度cd 指定当前目录cdedit 启动用户菜单、控件回调函数设计工具cdf2rdf 复数特征值对角阵转为实数块对角阵ceil 向正无穷取整cell 创建元胞数组cell2struct 元胞数组转换为构架数组celldisp 显示元胞数组内容cellplot 元胞数组内部结构图示char 把数值、符号、内联类转换为字符对象chi2cdf 分布累计概率函数chi2inv 分布逆累计概率函数chi2pdf 分布概率密度函数chi2rnd 分布随机数发生器chol Cholesky分解clabel 等位线标识cla 清除当前轴class 获知对象类别或创建对象clc 清除指令窗clear 清除内存变量和函数clf 清除图对象clock 时钟colorcube 三浓淡多彩交叉色图矩阵colordef 设置色彩缺省值colormap 色图colspace 列空间的基close 关闭指定窗口colperm 列排序置换向量comet 彗星状轨迹图comet3 三维彗星轨迹图compass 射线图compose 求复合函数cond （逆）条件数condeig 计算特征值、特征向量同时给出条件数condest 范 -1条件数估计conj 复数共轭contour 等位线contourf 填色等位线contour3 三维等位线contourslice 四维切片等位线图conv 多项式乘、卷积cool 青紫调冷色图copper 古铜调色图cos 余弦cosh 双曲余弦cot 余切coth 双曲余切cplxpair 复数共轭成对排列csc 余割csch 双曲余割cumsum 元素累计和cumtrapz 累计梯形积分cylinder 创建圆柱D ddblquad 二重数值积分deal 分配宗量deblank 删去串尾部的空格符dec2base 十进制转换为X进制dec2bin 十进制转换为二进制dec2hex 十进制转换为十六进制deconv 多项式除、解卷delaunay Delaunay 三角剖分del2 离散Laplacian差分demo Matlab演示det 行列式diag 矩阵对角元素提取、创建对角阵diary Matlab指令窗文本内容记录diff 数值差分、符号微分digits 符号计算中设置符号数值的精度dir 目录列表disp 显示数组display 显示对象内容的重载函数dlinmod 离散系统的线性化模型dmperm 矩阵Dulmage-Mendelsohn 分解dos 执行DOS 指令并返回结果double 把其他类型对象转换为双精度数值drawnow 更新事件队列强迫Matlab刷新屏幕dsolve 符号计算解微分方程E eecho M文件被执行指令的显示edit 启动M文件编辑器eig 求特征值和特征向量eigs 求指定的几个特征值end 控制流FOR等结构体的结尾元素下标eps 浮点相对精度error 显示出错信息并中断执行errortrap 错误发生后程序是否继续执行的控制erf 误差函数erfc 误差补函数erfcx 刻度误差补函数erfinv 逆误差函数errorbar 带误差限的曲线图etreeplot 画消去树eval 串演算指令evalin 跨空间串演算指令exist 检查变量或函数是否已定义exit 退出Matlab环境exp 指数函数expand 符号计算中的展开操作expint 指数积分函数expm 常用矩阵指数函数expm1 Pade法求矩阵指数expm2 Taylor法求矩阵指数expm3 特征值分解法求矩阵指数eye 单位阵ezcontour 画等位线的简捷指令ezcontourf 画填色等位线的简捷指令ezgraph3 画表面图的通用简捷指令ezmesh 画网线图的简捷指令ezmeshc 画带等位线的网线图的简捷指令ezplot 画二维曲线的简捷指令ezplot3 画三维曲线的简捷指令ezpolar 画极坐标图的简捷指令ezsurf 画表面图的简捷指令ezsurfc 画带等位线的表面图的简捷指令F ffactor 符号计算的因式分解feather 羽毛图feedback 反馈连接feval 执行由串指定的函数fft 离散Fourier变换fft2 二维离散Fourier变换fftn 高维离散Fourier变换fftshift 直流分量对中的谱fieldnames 构架域名figure 创建图形窗fill3 三维多边形填色图find 寻找非零元素下标findobj 寻找具有指定属性的对象图柄findstr 寻找短串的起始字符下标findsym 机器确定内存中的符号变量finverse 符号计算中求反函数fix 向零取整flag 红白蓝黑交错色图阵fliplr 矩阵的左右翻转flipud 矩阵的上下翻转flipdim 矩阵沿指定维翻转floor 向负无穷取整flops 浮点运算次数flow Matlab提供的演示数据fmin 求单变量非线性函数极小值点（旧版）fminbnd 求单变量非线性函数极小值点fmins 单纯形法求多变量函数极小值点（旧版）fminunc 拟牛顿法求多变量函数极小值点fminsearch 单纯形法求多变量函数极小值点fnder 对样条函数求导fnint 利用样条函数求积分fnval 计算样条函数区间内任意一点的值fnplt 绘制样条函数图形fopen 打开外部文件for 构成for环用format 设置输出格式fourier Fourier 变换fplot 返函绘图指令fprintf 设置显示格式fread 从文件读二进制数据fsolve 求多元函数的零点full 把稀疏矩阵转换为非稀疏阵funm 计算一般矩阵函数funtool 函数计算器图形用户界面fzero 求单变量非线性函数的零点G ggamma 函数gammainc 不完全函数gammaln 函数的对数gca 获得当前轴句柄gcbo 获得正执行"回调"的对象句柄gcf 获得当前图对象句柄gco 获得当前对象句柄geomean 几何平均值get 获知对象属性getfield 获知构架数组的域getframe 获取影片的帧画面ginput 从图形窗获取数据global 定义全局变量gplot 依图论法则画图gradient 近似梯度gray 黑白灰度grid 画分格线griddata 规则化数据和曲面拟合gtext 由鼠标放置注释文字guide 启动图形用户界面交互设计工具H hharmmean 调和平均值help 在线帮助helpwin 交互式在线帮助helpdesk 打开超文本形式用户指南hex2dec 十六进制转换为十进制hex2num 十六进制转换为浮点数hidden 透视和消隐开关hilb Hilbert矩阵hist 频数计算或频数直方图histc 端点定位频数直方图histfit 带正态拟合的频数直方图hold 当前图上重画的切换开关horner 分解成嵌套形式hot 黑红黄白色图hsv 饱和色图I iif-else-elseif 条件分支结构ifft 离散Fourier反变换ifft2 二维离散Fourier反变换ifftn 高维离散Fourier反变换ifftshift 直流分量对中的谱的反操作ifourier Fourier反变换i, j 缺省的"虚单元"变量ilaplace Laplace反变换imag 复数虚部image 显示图象imagesc 显示亮度图象imfinfo 获取图形文件信息imread 从文件读取图象imwrite 把imwrite 把图象写成文件ind2sub 单下标转变为多下标inf 无穷大info MathWorks公司网点地址inline 构造内联函数对象inmem 列出内存中的函数名input 提示用户输入inputname 输入宗量名int 符号积分int2str 把整数数组转换为串数组interp1 一维插值interp2 二维插值interp3 三维插值interpn N维插值interpft 利用FFT插值intro Matlab自带的入门引导inv 求矩阵逆invhilb Hilbert矩阵的准确逆ipermute 广义反转置isa 检测是否给定类的对象ischar 若是字符串则为真isequal 若两数组相同则为真isempty 若是空阵则为真isfinite 若全部元素都有限则为真isfield 若是构架域则为真isglobal 若是全局变量则为真ishandle 若是图形句柄则为真ishold 若当前图形处于保留状态则为真isieee 若计算机执行IEEE规则则为真isinf 若是无穷数据则为真isletter 若是英文字母则为真islogical 若是逻辑数组则为真ismember 检查是否属于指定集isnan 若是非数则为真isnumeric 若是数值数组则为真isobject 若是对象则为真isprime 若是质数则为真isreal 若是实数则为真isspace 若是空格则为真issparse 若是稀疏矩阵则为真isstruct 若是构架则为真isstudent 若是Matlab学生版则为真iztrans 符号计算Z反变换J j , K kjacobian 符号计算中求Jacobian 矩阵jet 蓝头红尾饱和色jordan 符号计算中获得 Jordan标准型keyboard 键盘获得控制权kron Kronecker乘法规则产生的数组L llaplace Laplace变换lasterr 显示最新出错信息lastwarn 显示最新警告信息leastsq 解非线性最小二乘问题（旧版）legend 图形图例lighting 照明模式line 创建线对象lines 采用plot 画线色linmod 获连续系统的线性化模型linmod2 获连续系统的线性化精良模型linspace 线性等分向量ln 矩阵自然对数load 从MAT文件读取变量log 自然对数log10 常用对数log2 底为2的对数loglog 双对数刻度图形logm 矩阵对数logspace 对数分度向量lookfor 按关键字搜索M文件lower 转换为小写字母lsqnonlin 解非线性最小二乘问题lu LU分解M mmad 平均绝对值偏差magic 魔方阵maple &nb, sp。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
流函数- 涡量法的二维方腔流数值模拟基本方程：① 在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S 方程为221Re u v tx y ξξξξψξ∂∂∂⎧++=∇⎪∂∂∂⎨⎪∇=-⎩其中,u v y xψψ∂∂==-∂∂ Re 为雷诺数差分格式：采用FTCS 格式有：1,,1,1,,1,11,,1,,1,,1,,2222122Re ()()n n n nn nn n n n n ni j i ji j i ji j i j i j i j i j i j i j i j nni ji juvtxyx y ξξξξξξξξξξξξ++-+-+-+-⎛⎫----+-+=--++ ⎪ ⎪∆∆∆∆∆⎝⎭1,,1,,1,,1,2222=0()()i j i j i j i j i j i j i j x y ψψψψψψξ+-+--+-+++∆∆,1,11,1,,,,22n n n n i j i j i j i jn i ji j u v yxψψψψ+-+---==-∆∆对于本问题，将方腔四边同时分为n 等分，则有x y h ∆=∆=∆ 故()()1,1,,1,1,1,,,1,1,,,1,124+()2Re n n n n ni j i j i j i j i j n n n n n n n ni ji ji j i j i j i j i j i j th u v h ξξξξξξξξξξξ+-+-++-+-⎡⎤+++-∆∆⎡⎤=--+-+⎢⎥⎣⎦∆⎢⎥⎣⎦21,1,,1,1,1,,,()4n n n n n i j i j i j i j i j n n n i ji ji j h ψψψψξψψψ+-+-+⎡⎤++++⨯∆=+-⎢⎥⎢⎥⎣⎦,1,11,1,,,,22n n n n i j i j i j i jni ji j uv hhψψψψ+-+---==-∆∆② 在直角坐标系下，不可压定常流体所满足的流函数涡量形式的N-S 方程为221Re =0u v x y ξξξψξ∂∂⎧+=∇⎪∂∂⎨⎪∇+⎩其中,u v y xψψ∂∂==-∂∂ Re 为雷诺数差分格式：采用FTCS 格式有：1,1,,1,11,,1,,1,,1,22221+=22Re ()()i j i ji j i j i j i j i j i j i j i j n i ji ju vxyx y ξξξξξξξξξξ+-+-+-+-+⎛⎫---+-++ ⎪ ⎪∆∆∆∆⎝⎭1,,1,,1,,1,2222=0()()i j i j i j i j i j i j i j x y ψψψψψψξ+-+--+-+++∆∆,1,11,1,,,,22i j i j i j i ji j i j u v yxψψψψ+-+---==-∆∆对于本问题，将方腔四边同时分为n 等分，则有x y h ∆=∆=∆，则有即()()1,1,,1,11,,,1,1,,,1,1,Re =+48n n n n i j i j i j i j n n n n n n n n n i j i j i j i j i j i j i j i j i j h u v ξξξξξξξξξξξ+-+-++-+-⎧⎫+++⨯∆⎪⎪⎡⎤+----⎨⎬⎣⎦⎪⎪⎩⎭21,1,,1,1,1,,,()4n n n n n i j i j i j i j i j n n n i ji ji j h ψψψψξψψψ+-+-+⎡⎤++++⨯∆=+-⎢⎥⎢⎥⎣⎦,1,11,1,,,,22n n n n i j i j i j i jni ji j uv hhψψψψ+-+---==-∆∆边界条件：在腔体的两侧和顶边，*22()()s s s s h ψψψξ=-=-∆ （第二式由泰勒级数展开得到）在底边*22()()s s s s h h ψψψξ=-+∆=-∆ （第二式由泰勒级数展开得到）其中s 代表边界，*s 代表与边界相邻的节点。
而,,x s h y s ∆⎧∆=⎨∆⎩当代表左右两边当代表上下两边x y h ∆=∆=∆即2,1,1,22(+)()j j j h h ψψξ-∆=-∆,2,1,122()()i i i h ψψξ-=-∆,1,1,22()()I j I j I j h ψψξ++-=-∆,,1,122()()i J i J i J h ψψξ++-=-∆Matlab 程序为：① 不可压非定常流体clear;%参数设置Re=10; %雷诺数取10,100，500，1000 L=1; %空穴几何尺寸 n=100;dh=L/n;%delta hdt=1e-4; %时间步长 psi=zeros(n+1,n+1); xi=zeros(n+1,n+1); rho=1;for k=1:1000000 err=0;%边界条件 for i=2:nxi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2; endfor j=2:nxi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2; xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2; end%控制方程 for i=2:n for j=2:nu(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh); v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));err1=(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+xi(i,j)*dh^2)/4-psi(i,j);psi(i,j)=psi(i,j)+rho*err1;err2=dt*(-dh/2*(u(i,j)*(xi(i+1,j)-xi(i-1,j)) ...+v(i,j)*(xi(i,j+1)-xi(i,j-1))) ...+(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1)-4*xi(i,j))/Re)/dh^2; xi(i,j)=xi(i,j)+rho*err2;temp=max(abs(err1),abs(err2));if err<temperr=temp;endendendif (mod(k,1000)==0) %每千步显示结果kerrcontour(psi,100);%contour求迹线pause(0.5)endif err<1e-6break;endendkerrrhodtcontour(psi,100);Re=10时，k=9216，err=9.9957e-07，rho=1，dt=1.0000e-04；Re=100时，k=10043，err=9.9973e-07，rho=1，dt=1.0000e-03；Re=500时，k=11275，err=9.9948e-07，rho=1，dt=0.0100;Re=1000时，k=16458，err=9.9983e-07，rho=1，dt=0.0100；②不可压定常流体clear;Re=10; %雷诺数取100，500，1000L=1; %空穴几何尺寸n=100;dh=L/n;%delta hpsi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1.0;for k=1:100000err=0;for i=2:nxi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2;endfor j=2:nxi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2;xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2;endfor i=2:nfor j=2:nu(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh);v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));err1=(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+xi(i,j)*dh^2)/4 -psi(i,j);psi(i,j)=psi(i,j)+rho*err1;err2=(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1))/4 ...-Re*dh*(u(i,j)*(xi(i+1,j)-xi(i-1,j))+v(i,j)*(xi(i,j+1)-xi(i,j-1)) )/8-xi(i,j);xi(i,j)=xi(i,j)+rho*err2;temp=max(abs(err1),abs(err2));if err<temperr=temp;endendendif err<1e-6endendkerrrho%psicontour(psi,100);Re=10时，k=6445，err=9.9978e-07，rho=1.0; Re=100时，k =7533，err=9.9953e-07，rho=1.0;Re=500时,k =10707，err =9.9973e-07，rho=1.0；Re=1000时,在不调节松弛因子时，其发散了，通过减小其松弛因子得到k =27600，err=9.9970e-07，rho=0.6000；最后使用时间相关法再次求解该问题，此时在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S 方程为221Re u v tx y t ξξξξξψξ∂∂∂⎧++=∇⎪∂∂∂⎪⎨∂⎪=∇+⎪∂⎩其中,u v y xψψ∂∂==-∂∂ Re 为雷诺数差分格式：采用FTCS 格式有：1,,1,1,,1,11,,1,,1,,1,,2222122Re ()()n n n nn nn n n n n ni j i ji j i ji j i j i j i j i j i j i j i j nni ji juvtxyx y ξξξξξξξξξξξξ++-+-+-+-⎛⎫----+-+=--++ ⎪ ⎪∆∆∆∆∆⎝⎭1,,1,,1,,1,,1,2222()()n nn n n n n n i j i ji j i j i j i j i j i j n i jtx y ψψψψψψψψξ++-+---+-+=++∆∆∆,1,11,1,,,,22n n n n i j i j i j i jn i j i j u v yxψψψψ+-+---==-∆∆对于本问题，将方腔四边同时分为n 等分，则有x y h ∆=∆=∆ 故()()1,1,,1,1,1,,,1,1,,,1,124+()2Re n n n n ni j i j i j i j i j n n n n n n n ni ji ji j i j i j i j i j i j th u v h ξξξξξξξξξξξ+-+-++-+-⎡⎤+++-∆∆⎡⎤=--+-+⎢⎥⎣⎦∆⎢⎥⎣⎦1,1,,1,1,1,,,24()n n n n ni j i j i j i j i j n n ni j i j i j t h ψψψψψψψξ+-+-+⎛⎫+++-=+∆+ ⎪ ⎪∆⎝⎭,1,11,1,,,,22n n n n i j i j i j i jni ji j uv hhψψψψ+-+---==-∆∆Matlab 程序为：clear;%参数设置Re=1000; %雷诺数取100，500，1000 L=1; %空穴几何尺寸 n=100;dh=L/n;%delta hdt=5e-5; %时间步长 psi=zeros(n+1,n+1); xi=zeros(n+1,n+1); rho=1.0;for k=1:1000000 err=0; for i=2:nxi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2; endfor j=2:nxi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2; xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2; endfor i=2:nfor j=2:nu(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh);v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));err1=dt*((psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)-4*psi(i,j)) /dh^2+xi(i,j));psi(i,j)=psi(i,j)+rho*err1;err2=dt/dh^2*(-dh/2*(u(i,j)*(xi(i+1,j)-xi(i-1,j)) ...+v(i,j)*(xi(i,j+1)-xi(i,j-1))) ...+(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1)-4*xi(i,j))/Re);xi(i,j)=xi(i,j)+rho*err2;temp=max(abs(err1),abs(err2));if err<temperr=temp;endendendif (mod(k,1000)==0)kerrcontour(psi,100);pause(0.5)endif err<1e-6break;endendkerrcontour(psi,100);对于该方法，在此处难以收敛，可能需继续调试，这里就不再讨论了。