流场粒子追踪精度数值实验

在计算流线,拉格朗日拟序结构等流场后处理时,我们常常需要计算无质量的粒子在流场中迁移时的轨迹,无质量意味着粒子的速度为流场当地的速度。此时,求解粒子的位移这个问题是一个非常简单的常微分方程问题。
假设流场中存在 i 个粒子,那么对于每个粒子有 v i = x i 初始条件为: x i ( t 0 ) 假设流场中存在i个粒子,那么对于每个粒子有 \\v_i = x_i \\初始条件为:x_i(t_0) 假设流场中存在i个粒子,那么对于每个粒子有vi=xi初始条件为:xi(t0)
1、欧拉向前差分
欧拉向前差分是最简单的方法,其主要步骤为:
x i t + d t = x i + v i t ∗ d t v i t 为 t 时刻在 x i 处的流场速度 \\x_{i}^{t+dt} = x_{i}+v_{i}^{t}*dt \\ v_{i}^{t}为t时刻在x_i处的流场速度 xit+dt=xi+vitdtvitt时刻在xi处的流场速度
显然,该方法的精度取决于时间步长,时间步长越大,误差越大。
该方法需要注意几个事项:
(1)流场插值处理
我们实现知道的流场数据是时间和空间离散的,对于任意一个时刻,我们只知道空间离散点处的流场速度,加入粒子的位置恰好不在这些点上,我们就需要通过插值得到其速度。
matlab提供了这样的插值函数,interp2()。
假设我们现在已知流场的空间离散点矩阵x_f;y_f以及速度U,V;
我们想通过插值获得,ftle_x,ftle_y处粒子的速度ux,uy

ux = interp2(x_f,y_f,U,ftle_x,ftle_y,'linear');
uy = interp2(x_f,y_f,V,ftle_x,ftle_y,'linear');
%%上面的代码中ftle_x,ftle_y可以是矩阵形式,这样得到的速度ux,uy也是矩阵    

那么接下来粒子迁移就十分简单了

for t = f_time:dt:(T+f_time)
    %获取t时刻的速度数据
    U = reshape(vel(1:Nx*Ny,t),Ny,Nx);
    V = reshape(vel(Nx*Ny+1:end,t),Ny,Nx);
    %%计算粒子迁移t-->t+dt
    ux = interp2(x_f,y_f,U,ftle_x,ftle_y,'linear');
    uy = interp2(x_f,y_f,V,ftle_x,ftle_y,'linear');

    X = ftle_x + dt*tstep*ux;
    Y = ftle_y + dt*tstep*uy;
    ftle_x=X;
    ftle_y=Y;
end

另外,好像对于这个问题不能用梯型方法?梯型方法的步骤如下:
y n + 1 = y n + h 2 ( f ( x n , y n ) + f ( x n + 1 , y n + 1 ) y_{n+1}=y_n+\frac{h}{2}\bigl(f\bigl(x_n,y_n\bigr)+f\bigl(x_{n+1},y_{n+1}) yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+1)
主要问题在于我没法知道下一步粒子在哪,即这里的x(n+1),y(n+1),因为我就是在求这个东西?
暂时没想明白怎么用梯型方法。
2、龙格-库塔法
这个就比较玄学了,这个我也没弄清原理,只是知道matlab里有这个函数可以用,并且还不错。
由于龙格-库塔法需要两个时间步中间的时刻,我们不仅需要空间插值,还需要时间插值。
插值的方法也很简单,matlab提供了相应的函数
(1)首先我们需要获得粒子的时间,空间离散结果

for k = 1:size(vel,2)
    XX(:,:,k) = X0;
    YY(:,:,k) = Y0;%%这里是给定了每个时刻的离散空间位置,我们采用欧拉法求解流场,所以离散点是固定的
    u(:,:,k) = reshape(vel(1:ny*nx,k),[ny nx]);
    v(:,:,k) = reshape(vel(ny*nx+1:end,k),[ny nx]);
    TT(:,:,k) = t(k)*ones(ny,nx);%%这里存储着每个粒子的时间信息,因为全场内每个粒子的时间是相同的,所以用了ones()函数
end

(2)对其进行转置(该步骤不是必要的)

%% 该部分的操作相当于对一个三维矩阵的前两维进行转置
P = [2 1 3];    
Xprime = permute(XX,P);
Yprime = permute(YY,P);
Tprime = permute(TT,P);
Uprime = permute(u,P);
Vprime = permute(v,P);

(3)时空间插值

u_interp = griddedInterpolant(Xprime,Yprime,Tprime,Uprime,'linear');
v_interp = griddedInterpolant(Xprime,Yprime,Tprime,Vprime,'linear'); 

(4)ode45()函数求解

 [t, X_out] = ode45(@velocity_interp_func,intspan(tspan(k),T),[X0(i,j),Y0(i,j)]);
 %% 这里ode45函数需要输入三个参数
 @velocity_interp_func这是句柄,包含着方程的信息
 intspan(tspan(k),T) 这个是从某个时刻计算到某个时刻
 这里具体含义是从tspan(k)时刻推进到tspan(k)+T时刻
 [X0(i,j),Y0(i,j)]是初始场

根据这个函数就可以计算处初始位置在[X0(i,j),Y0(i,j)],初始时刻为tspan(k),经过T时间的运动后,最终的位置X_out。X_out包含了每个时刻对应的粒子位置,t即是这些时刻的时间信息。
但是我们只需要最后一个时刻的信息,因此

X_adv(i,j) = X_out(end,1);
Y_adv(i,j) = X_out(end,2);  

需要注意的是,只是计算单独一个粒子的轨迹,要想计算全场的粒子,需要对其进行遍历,然鹅经过实操发现,这个过程还是很耗时的。
这个方法的具体运行原理我还没有搞明白,因为最简单的梯型法现在还没想清楚,所以仅放在这里。

注意:上面两种方法都是利用离散数据得到粒子迁移的,这里我们通过解析解获得时间间隔为0.1下的0:16s时刻共161个流场数据,流场网格划分为均匀划分,nx =201;ny = 101;
流场设置为double_gyre流场
我们计算的问题是:
流场初始0时刻布置201*101个粒子,计算10s时刻时这些粒子的位置。
3、对照组
为了方便对比上面两种方法的精度,我设置了一组对照组,由于流场设置为double_gyre流场,这个流场是存在解析解的,因此我们可以获得任意时刻,任意位置处的精确速度信息。
通过设置时间步长非常小,并且采用欧拉向前推进,可以将这个方法得到的结果精度提升到很高。
我计算了两种精度的结果1、dt=0.01s,2、dt=0.001s

然后将上面两种方法得到的结果与该方法进行对比,从而得到两种方法的误差范围。
4、实验结果

figure(1)
dx = sigma_x - X_ana;
dy = sigma_y - Y_ana;
e_ds = sqrt(dx.^2+dy.^2);
pcolor(e_ds)
shading interp

figure(2)
dx = X_adv - X_ana;
dy = Y_adv - Y_ana;
o_ds = sqrt(dx.^2+dy.^2);
pcolor(o_ds)
shading interp

欧拉法与dt=0.01解析解对比,后者的时间精度是前者的十倍
结果如下:
在这里插入图片描述

ode45方法结果与dt=0.01解析方法对比
结果如下:
在这里插入图片描述
可以看出,欧拉方法的误差要大于ode45方法,但是其实差距不是很大。

我们对误差进行全场平均:
欧拉:0.71
ode45:0.67

误差最大值:
欧拉:2.18
ode45:1.97
这里其实误差还是蛮大的,因为流场范围才是2*1.
还有一个比较神奇的点是,这个误差云图很像ftle图,这是因为ftle图反应了粒子拉伸扭曲的程度,而恰好在这些位置处,误差非常敏感,比较有意思。

最后我们对比了dt=0.01和dt =0.001时解析方法之间的对比结果
其误差云图如下:
在这里插入图片描述
可以看出,dt=0.01时时间精度就已经很高了,与dt=0.001的结果很接近。

补充:总觉得上面的对比还是不够明显,于是又进行了几组测试
dt=0.02;dt=0.05
测试结果,这里仅对比误差平均值:
dt=0.05
欧拉:0.082
ode45:0.091
欧拉方法更接近dt=0.05
dt=0.02
欧拉:0.13
ode45:0.038
ode45方法更接近dt=0.02

可以看出ode45方法其实精度超过了dt=0.05时的解析方法。
对比到此结束,仍然需要指出。虽然平均误差看起来比较大,实际上这些是主要由所谓的FTLE场的脊位置处的粒子贡献的。大多数位置处的粒子误差水平还是比较小的。
下图左为欧拉,右为ode45,蓝色区域的误差还是比较小的。
在这里插入图片描述
欧拉方法的误差矩阵:
在这里插入图片描述

ode45方法的误差矩阵:
在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/30510.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

Java版本+企业电子招投标系统源代码之电子招投标系统建设的重点和未来趋势

计算机与网络技术的不断发展,推动了社会各行业信息化的步伐。时至今日,电子政务、电子商务已经非常普及,云计算、大数据、工业4.0、“互联网”等发展理念也逐步深入人心,如何将传统行业与互联网科技有效结合起来,产生1…

Vue实现元素沿着坐标数组移动,超出窗口视图时页面跟随元素滚动

一、实现元素沿着坐标数组移动 现在想要实现船沿着下图中的每个河岸移动。 实现思路: 1、将所有河岸的位置以 [{x: 1, y: 2}, {x: 4, y: 4}, …] 的形式保存在数组中。 data() {return {coordinateArr: [{ x: 54, y: 16 }, { x: 15, y: 31 }, { x: 51, y: 69 }…

升级Nginx

目录 前言 一、升级Nginx 1)首先在官网下载一个新版本的Nginx 2)首先将下载的压缩包进行解包 3)进入已解包的目录中 4)配置安装路径 5)make 6)备份原来Nginx的资源 7)重启Nginx服务 8&#…

【2023最全教程】Web自动化测试怎么做?Web自动化测试的详细流程和步骤

一、什么是web自动化测试 自动化(Automation)是指机器设备、系统或过程(生产、管理过程)在没有人或较少人的直接参与下,按照人的要求,经过自动检测、信息处理、分析判断、操纵控制,实现预期的目…

毕业季Android开发面试,有哪些常见的题?

前言 对于计算机行业早已烂大街,随之而来的毕业季。还会有大批的程序员涌进来,而我们想要继续进入Android开发岗位的人员,最先考虑的是面试。面试题是我们决定踏进工作的重要环节。 对于刚毕业的实习生来说,如何在应聘中脱颖而出…

LightningChart .NET 10.5.1 Crack LightningChart 2023

LightningChart .NET v.10.5.1 已经发布! DataCursor 和 3D TransparencyRenderMode 现在可用。 为所有 3D、Polar 和 Smith 系列启用 DataCursor 在早期阶段,LightningChart 提供了不同的工具,需要用户编写额外的代码才能启用数据跟踪功能。…

控制您的数据:Web3私有链为数据主权带来的突破性变革

在数字化时代,数据已经成为企业和个人最宝贵的资产之一。然而,随着大规模数据泄露和滥用事件的频发,数据主权和隐私保护成为了备受关注的问题。在这个背景下,Web3私有链的出现为数据主权带来了一场突破性的变革。 首先&#xff0c…

风景类Midjourney prompt提示词

稳定输出优美风景壁纸的Midjourney prompt提示词。 1\在夏夜,有淡蓝色的星空,海边,流星,烟花,海滩上全是蓝色的玫瑰和绿色的植物,由Ivan Aivazovsky和Dan Mumford,趋势在cgsociety,…

windows2022证书配置.docx

Windows证书的配置 要求两台主机,一台作为域,一台进入域 按要求来选择角色服务 确认之后安装 安装完以后配置证书服务 选择服务 按要求配置 注:此处不用域用户登陆无法使用企业CA 按要求来 创建新的私钥 这几处检查无误后默认即可 有效期…

AJAX概述

1.1什么是AJAX. Ajax即AsynchronousJavascript And XML:异步数据回调。 使用Ajax技术网页应用能够快速地将更新呈现在用户界面上,不需要重载(刷新)整个页面【只刷新局部】,这使得程序能够更快地回应用户的操作。、 1…

2023年5月青少年机器人技术等级考试理论综合试卷(四级)

青少年机器人技术等级考试理论综合试卷(四级)2023.6 分数: 100 题数: 30 一、 单选题(共 20 题, 共 80 分) 1.Arduino C 语言, 部分程序如下, 串口监视器输出结果是“D”时, 变量 i …

【集群】Haproxy搭建Web群集

文章目录 一、Haproxy 相关概念1. Haproxy 的概述2. Haproxy 的主要特性3. 常见的 Web 集群调度器4. 常见的应用分析4.1 LVS 应用4.2 Haproxy 应用4.3 LVS、Nginx、Haproxy的区别 5. Haproxy 调度算法原理5.1 roundrobin5.2 static-rr5.3 leastconn5.4 source5.5 uri5.6 url_pa…

SpringBoot + Vue前后端分离项目实战 || 二:Spring Boot后端与数据库连接

系列文章: SpringBoot Vue前后端分离项目实战 || 一:Vue前端设计 文章目录 新建Spring后台项目添加依赖 新建数据库IDEA 连接数据库IDEA 自动创建类实体定义数据传递至前端的格式 B站视频讲解:2023全网最简单但实用的SpringBootVue前后端分离…

RTC

文章目录 前言驱动应用程序运行 前言 RTC(Real Time Clock,实时时钟)是个常用的外设,通过 RTC 我们可以知道日期和时间信息,因此在需要记录时间的场合就需要实时时钟。 可以使用专用的实时时钟芯片来完成此功能&#…

扫雷小游戏【C语言】

目录 前言 一、基本实现逻辑 二、实现步骤 1. 我们希望在进入游戏时有一个菜单让我们选择 2. 我们希望可以重复的玩(一把玩完了还可以接着玩) 3. 采用多文件形式编程 4.要扫雷先得有棋盘(创建棋盘R*N) 5.初始化棋盘 6.打…

【网络安全】深入解析 PHP 代码审计技术与实战

前言 登录某个网站并浏览其页面时,注意到了一些看起来不太对劲的地方。这些迹象可能是该网站存在漏洞或被黑客入侵的标志。为了确保这个网站的安全性,需要进行代码审计,这是一项专门针对软件代码进行检查和分析的技术。在本文中,…

一、Docker介绍

学习参考:尚硅谷Docker实战教程、Docker官网、其他优秀博客(参考过的在文章最后列出) 目录 前言一、Docker是什么?二、Docker能干撒?三、容器虚拟化技术 和 虚拟机有啥区别?1.虚拟机2.容器虚拟化技术3.对比4.Docker为啥比VM虚拟机…

献给蓝初小白系列(二)——Liunx应急响应

1、Linux被入侵的症状​​ ​​https://blog.csdn.net/weixin_52351575/article/details/131221720​​ 2、Linux应急措施 顺序是:隔离主机--->阻断通信--->清除病毒--->可疑用户--->启动项和服务--->文件与后门--->杀毒、重装系统、恢复数据 …

AAC ADTS格式分析

标题 1.AAC简介2. AAC ADTS格式分析2.1 adts_fixed_header详细介绍2.2 adts_variable_header详细介绍 1.AAC简介 AAC音频格式:Advanced Audio Coding(⾼级⾳频解码),是⼀种由MPEG-4标准定义的有损⾳频压缩格式,由Fraunhofer发展,Dolby, Sony…

vue3 + TS + elementplus + pinia实现后台管理系统左侧菜单联动实现 tab根据路由切换联动内容

效果图&#xff1a; home.vue页面代码 <template><el-container><el-aside width"collapse ? 200px : 70px"><el-button color"#626aef" click"collapseToggle()"><el-icon><Expand v-if"collapse"…