CN111695309A - 基于统计动力学的高性能大规模流固耦合流体仿真方法 - Google Patents
基于统计动力学的高性能大规模流固耦合流体仿真方法 Download PDFInfo
- Publication number
- CN111695309A CN111695309A CN202010488622.9A CN202010488622A CN111695309A CN 111695309 A CN111695309 A CN 111695309A CN 202010488622 A CN202010488622 A CN 202010488622A CN 111695309 A CN111695309 A CN 111695309A
- Authority
- CN
- China
- Prior art keywords
- fluid
- space
- scale
- simulation
- solid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000012530 fluid Substances 0.000 title claims abstract description 110
- 239000007787 solid Substances 0.000 title claims abstract description 104
- 238000000034 method Methods 0.000 title claims abstract description 94
- 238000004088 simulation Methods 0.000 title claims abstract description 92
- 230000008878 coupling Effects 0.000 title claims abstract description 62
- 238000010168 coupling process Methods 0.000 title claims abstract description 62
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 62
- 238000013507 mapping Methods 0.000 claims abstract description 41
- 238000005457 optimization Methods 0.000 claims abstract description 12
- 238000005070 sampling Methods 0.000 claims abstract description 10
- 238000005315 distribution function Methods 0.000 claims description 12
- 238000009826 distribution Methods 0.000 claims description 9
- 238000007654 immersion Methods 0.000 claims description 8
- 238000012417 linear regression Methods 0.000 claims description 7
- 230000003993 interaction Effects 0.000 claims description 6
- 238000005259 measurement Methods 0.000 claims description 6
- 238000003672 processing method Methods 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 claims description 5
- 238000012549 training Methods 0.000 claims description 5
- 230000002457 bidirectional effect Effects 0.000 claims description 4
- 238000005094 computer simulation Methods 0.000 claims description 2
- 238000010606 normalization Methods 0.000 claims description 2
- 230000005540 biological transmission Effects 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 15
- 238000010586 diagram Methods 0.000 description 13
- 230000000694 effects Effects 0.000 description 12
- 230000008569 process Effects 0.000 description 11
- 239000006185 dispersion Substances 0.000 description 10
- 238000009877 rendering Methods 0.000 description 9
- 230000008859 change Effects 0.000 description 7
- 239000011159 matrix material Substances 0.000 description 5
- 239000002245 particle Substances 0.000 description 5
- 238000013461 design Methods 0.000 description 4
- 238000013508 migration Methods 0.000 description 4
- 230000005012 migration Effects 0.000 description 4
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 3
- 150000001875 compounds Chemical class 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241000283690 Bos taurus Species 0.000 description 1
- 208000012322 Raynaud phenomenon Diseases 0.000 description 1
- 238000007664 blowing Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010924 continuous production Methods 0.000 description 1
- 230000001808 coupling effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 239000000779 smoke Substances 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 230000017105 transposition Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Computing Systems (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Processing Or Creating Images (AREA)
Abstract
本发明提出了一种基于统计动力学的高性能大规模流固耦合流体仿真方法。该方法基于统计动力学，采用高阶非正交中心矩弛豫格子玻尔兹曼模型，首次提出基于优化的自适应高阶松弛参数控制，遵循伽利略不变性，可以在高雷诺数下仿真流固耦合的强湍流现象。为了提升仿真稳定性，本发明构造了统计动力学介观模型和宏观流体物理模型的通用映射对应关系，能够实现在时间尺度上的自适应流固耦合仿真，在空间尺度上对求解区域进行自适应尺度的采样和连续尺度展开，大大提升了计算效率。同时，本发明提出的方法具有良好的并行性和可扩展性。在高分辨率下，通过将其扩展到多节点多GPU的大型服务器设备上，可以实现高性能大规模场景的流固耦合仿真。
Description
技术领域
本发明涉及一种基于统计动力学方法的高性能大规模流固耦合的流体仿真方法，是一种全新的稳定的基于统计动力学的介观模型仿真方法。
背景技术
流体仿真和流固耦合在计算流体力学(Computational Fluid Dynamics，简称CFD)和计算机图形学(Computer Graphics，简称CG)两个领域都有着悠久的历史。
在CFD领域，由于目前的低阶求解格式的稳定性和可靠性，被广泛应用于实际工程计算中。但是低阶格式存在大的数值耗散和色散误差，对于复杂问题，比如湍流，流固耦合等问题，必须采用低耗散和色散的高阶形式。近几年发展的高阶形式有：有限差分高阶格式、间断Galerkin有限元法、ENO/WENO有限体积法、有限谱差分法等。高质量的计算网格是CFD求解的前提条件，也是影响CFD求解精度最重要的因素之一。而生成网格需要大量的工作量，导致CFD求解效率低下。流体仿真在计算机图形学中的发展超过二十几年。在计算机图形学领域，流体仿真和流固耦合的发展，使得流体仿真在电影特效制作以及游戏中有着广泛的应用。一个稳定的精确的流固耦合的流体仿真方法对于研究流体及其应用是非常必要的。如何在有限的计算资源和高效率的情况下做到数值稳定，精确和灵活的流体仿真方法受到众多学者的青睐。高性能大规模流体仿真平台对于工业设计有着至关重要的作用。但是目前的流体仿真技术几乎很难满足上面的要求。
目前常用的CFD方法有有限差分方法(Finite Difference Method(FDM))，有限体积方法(Finite Volume Method(FVM))、粒子法(Particle Projection Method)、格子玻尔兹曼方法(Lattice Boltzmann Method(LBM))等方法。有限差分方法基于网格求解不可压缩N-S方程，存在较大的数值耗散，在准确度上没有达到很高的精度，后期发展中有学者提出加入湍流模型的方法，虽然可以增加更多的流体细节，但是大部分都不是遵循物理规律。并且，现有的直接方法在复杂的边界上的湍流仿真中存在较大缺陷，需要求解一个全局的线性方程组，效率比较低。有限体积法是将流体的欧拉控制方程在单元控制体内进行积分后离散求解。这种方法需要对求解空间进行网格化，需要很多的计算资源，计算效率低下，不易于广泛应用。粒子法是一种基于拉格朗日的近似方法，无需网格，需要对流体和固体的本身进行离散。光滑粒子流体动力学(Smoothed-particle hydrodynamics(SPH))方法虽然快，但是准确度难以保证，适合于追求视觉效果的场景。
LBM方法属于介观尺度的方法。该方法区别于其他CFD方法，没有直接求解流体的纳维-斯托克斯(N-S)方程，而是通过计算介观粒子间的迁移(streaming)和碰撞(collision)两个过程，从而模拟流体的运动。早期基于统计物理学发展出Bhatnagar–Gross–Krook(BGK)模型，逐渐发展成一种数值计算方法，其中有多尺度弛豫(MRT)模型和级联格子玻尔兹曼(Cascaded LBM)或者非正交中心矩驰豫格子玻尔兹曼模型(NO-CMLBM)模型。相比于N-S方程，LBM的主要特点是不需要求解全局方程，具有很好的局部特性，有利于并行实现，具体很高的计算效率。现有的BGK模型和MRT模型中，较低的粘度情况下，数值求解的稳定性很难保证。Cascaded LBM和NO-CMLBM模型由于模型本身提供的高阶参数不能达到高精度，也会导致数值扩散。
传统的LBM方法通常采用的是相同尺度的网格来求解，这样很难捕捉到空间和时间上流体连续变化的物理细节，除非用非常高分辨率的计算网格，但是又会导致需要庞大的计算资源。因此一些基于多重网格的方法提出来，用于解决这个问题，但是在多重网格的方法中，不同尺度的网格之间的变化是整数倍的，所以也同样很难描述流体在空间中的连续变化过程，不够灵活，同时基于多重网格的流体模拟，在湍流中会产生不连续的结构。
流固耦合是近年来计算机图形学和CFD领域的研究热点。在计算机图形学中，边界处理尤其是边界上满足不可压缩性的压力条件(压力需要提供正确的作用到固体运动上的力)做了很多研究工作。体素化边界和全欧拉方法在强湍流下的流固耦合和小时间步长会产生强的色散误差。基于高精度的网格优化的边界一致性的网格方法在CFD和计算机图形学领域应用广泛，但是其动态优化网格需要大量的计算代价；剪切单元格方法基于FVM思想有着较高的精度但是需要重构新鲜网格细胞(仿真中上一时刻在固体内部的网格，当前时刻在流体网格中)和死亡网格细胞(仿真中上一时刻在流体网格中，当前时刻在固体网格中)，尤其是在运动变化过大的流固耦合情况下，其重构方法的精度远远不够。基于SPH方法的流固耦合中依然存在固体边界上压力的边界条件准确性弱的问题。CFD领域中的浸入式边界方法，因其不需要解决新鲜网格问题和其强的局部性使其有着强并行性而被广泛被采用。由于介观流体和宏观固体运动目前没有一种通用的形式映射，所以浸入式边界方法往往很难做到真实的流固耦合的流体仿真。
发明内容
本发明的目的是：提出一种创新性的高阶参数优化模型，通过建立介观统计力学和宏观流体力学的通用映射，可以做到复杂湍流和动态边界的自适应空间和时间尺度上的多分辨率耦合仿真并且易于扩展到多节点多GPU的高性能大规模流固耦合仿真。
为了达到上述目的，本发明的技术方案是提供了一种基于统计动力学的高性能大规模流固耦合流体仿真方法，其特征在于，包括以下步骤：
步骤(1)流体仿真建模：在非正交中心矩驰豫格子玻尔兹曼模型(Non-OrthogonalCentral Moment LBM，以下简称NO-CMLBM)的基础上，基于六阶埃尔米特多项式近似N-S方程的NO-CMLBM提出目标度量方程
xk表示空间中k点，t表示LBM空间的时间，∈(xk，t)表示t时刻空间中k点的度量方程值，ρ表示LBM空间的密度，u表示LBM空间的速度，Π表示LBM空间二阶矩，δt(·)＝(·)t+Δt-(·)t，Δt表示真实物理空间的时间步长，k表示迭代次数的倍数，在多种不同仿真环境下优化目标度量方程得到足够多的训练数据，然后基于训练数据用线性回归的方法去自适应控制NO-CMLBM的高阶参数，本发明提出了高阶参数的局部表达式
其中，
表示弛豫时间，θ表示是线性回归得到的参数向量，T表示矩阵转置，sp表示LBM空间p点的状态向量，sp＝(ρp，||ρpup||，||Πp||，1)，ρp表示LBM空间p点的密度，up表示LBM空间p点的速度，Πp表示LBM空间p点的二阶矩，基于训练数据用线性回归拟合，得到NO-CMLBM的θ的具体数值,其高阶参数的控制方法可以应用于其他各种场景的仿真。
采用上述步骤，本发明大大减少了求解方法的耗散和色散误差，提高求解方法的准确性和稳定性。
在本步骤中，基于六阶埃尔米特多项式近似N-S方程的NO-CMLBM模型，本发明首次提出目标度量方程，遵循伽利略不变性，能够稳定的模拟高雷诺数下的强湍流。本发明设计高阶弛豫参数ri(i表示第i个高阶参数，i＝{9，...，26})和优化参数τi的关系：ri＝1/(τi+0.5)，提出NO-CMLBM的优化参数τi和求解方程中的耗散和色散误差存在一个最优的高阶参数满足最小的耗散和色散的数值误差的关系。然后优化目标度量方程，生成足够多数据，基于线性回归得到高阶参数的自适应计算方法。
步骤(2)通用映射构造：根据流体的性质，雷诺数在物理空间下保持不变，基于宏观基础物理量归一化到介观格子玻尔兹曼方程中，建立宏观流体物理量和介观流体物理量的通用映射关系，其中，宏观流体和介观流体的通用映射物理量包括空间步长、时间步长、求解空间尺度大小、密度、速度、单位体积力、粘度；
步骤(3)时间空间连续尺度自适应仿真：基于步骤(2)的通用映射关系，在时间和空间连续尺度上进行自适应仿真，构造连续尺度之间的映射关系；
步骤(4)高性能流固耦合仿真：使用浸入式边界处理方法，将流体方程求解的欧拉网格和固体运动求解的拉格朗日网格耦合，满足滑动和非滑动的边界条件，利用最小插值内核构建低耗散的流固边界，其中，耦合力方法基于步骤(2)的通用映射关系，计算介观流体和宏观刚体之间的耦合力的映射关系并且提供单向流固耦合和双向流固耦合流体仿真；
步骤(5)建立低分变率实时仿真和高性能大规模流固交互流体仿真平台：在NVIDIA RTX 2080Ti显卡上可以做到100×200×100分辨率的实时仿真，并且其易扩展可以在多节点多GPU上高效并行化实现。
优选地，步骤(2)中，所述通用映射关系为：
式中，t表示真实物理空间的时间，
表示LBM空间的参照速度，uref表示物理空间的参照速度，Δx表示真实物理空间的空间步长，
表示LBM空间的时间，Δt表示真实物理空间的时间步长，
表示LBM空间的时间步长，x表示真实物理空间的求解空间尺度大小，
表示LBM空间的求解空间尺度大小，
表示LBM空间的空间步长，ρ表示真实物理空间的密度，ρ0表示参照的物理空间的空气密度，
表示LBM空间的密度，u表示真实物理空间的速度，
表示LBM空间的速度，F表示真实物理空间的单位体积力，
表示LBM空间的单位体积力，v表示真实物理空间的粘度，
表示LBM空间的粘度。
优选地，流体在真实世界中，空间和时间上是连续变化的，之前方法在不同区域采用整数倍尺度过渡进行流体仿真，湍流情况下会损失流体的湍流结构，带来流体结构上的不连续。所以需要保证流体模拟仿真在空间上的尺度变化是一个连续的过程，因此，步骤(3)中本发明构造了空间上连续变化的尺度。由于流体结构的不连续现象一般是发生在不同尺度的边界部分，所以在不同尺度的边界部分，增加它们边界的耦合部分，使得流体能够在连续尺度的边界上过渡的更加自然。同时保证空间连续尺度之间的雷诺数一致，根据不同尺度的空间步长的倍数关系，通过通用映射关系中时间步长关系
推导不同尺度上迭代次数的倍数k，然后基于倍数k在不同尺度上完成时间同步一致，追踪固体的运动轨迹，在固体周围重建精细尺度网格，精细尺度网格和粗糙尺度网格部分区域耦合在一起形成完整的仿真空间。在精细尺度网格内部进行流固耦合计算，然后在其他粗糙尺度网格区域仅仿真流体运动。
时间尺度上基于每个时刻求出在介观尺度下的速度最大值
保持
在[0.18,0.22]之间，通过通用映射关系中速度映射关系
已知LBM空间的速度
和
以及真实物理空间的速度u求出uref，然后更新介观LBM空间的
分布，根据新的介观物理量
重构介观分布函数fi，继续进行流体仿真。
优选地，步骤(4)中，所述浸入式边界处理方法包括以下步骤：
首先优化固体网格，基于每个三角形网格采样一点或者泊松采样构造固体拉格朗日网格点集分布，随后将欧拉流体求解网格和固体拉格朗日网格耦合，满足滑动和非滑动的边界条件，利用最小插值内核
求解低耗散的流固边界力的分布，其中，r表示插值点到插值中心的距离。
其中非滑动的边界条件是通过速度的耦合边界条件推导出流固耦合中固体对流体运动作用力，计算方法为：
式中，Xp表示空间中p点的位置，Fs→f(Xp)表示空间p点处固体s对流体f的作用力，ub(Xp)表示表示固体边界b上p点的速度u，u(Xp)表示流体在p点的的速度。然后基于作用力与反作用力得到流体运动对固体运动的作用力，求解方法为：Ff→s(Xp)＝-Fs→f(Xp)
滑动的边界条件推导出流固耦合中固体运动对流体运动作用力，计算方法为：
式中，n表示固体边界b上p点的的法向量。流体运动对固体运动作用力的求解方法为：Ff→s(Xp)＝-Fs→f(Xp)。在固体运动(刚体运动)中所受流体作用的合力Fs以及合转矩τs是通过Fs＝∑pFf→s(Xp)Δsp，τs＝∑p(Xp-B)×Ff→s(Xp)Δsp计算得到，式中，B表示固体质心，Δsp表示离散积分面元。
优选地，步骤(5)中，利用格子玻尔兹曼方程的并行性，通过并行优化可以在NVIDIA RTX 2080Ti显卡上做到100×200×100分辨率的实时仿真；基于统一计算架构(CUDA)和消息传递接口(MPI)并行技术，由于本发明的可扩展性，优化并行实现的数据结构，可以在多节点多GPU的服务器上实现高性能大规模流固耦合的流体仿真。
本发明提出的高精度的高阶松弛参数控制的高阶非正交中心矩松弛模型，是一种全新的稳定的基于统计动力学的介观模型仿真方法，大大降低仿真过程的数值耗散和色散误差。基于介观和宏观物理量的通用映射，通过浸入式边界方法，可以在强湍流、大速度变化下做到复杂的流固耦合仿真，从而提供了一种高效、准确和鲁棒的流固耦合仿真方法。其高效性和扩展性使得本发明通过GPU并行优化可以在NVIDIA RTX 2080Ti显卡上做到100×200×100分辨率的实时仿真。并且本发明的方法自身有很强的可扩展性，基于统一计算架构(CUDA)和消息传递接口(MPI)并行技术，优化数据结构后在多节点多GPU的服务器上实现流固耦合的高性能大规模仿真。
提出的高性能大规模流固耦合的流体仿真方法，是一种全新的稳定的基于统计动力学的介观模型仿真方法。基于高精度的高阶松弛参数的高阶非正交中心矩松弛模型，构造介观与宏观物理量的通用映射，可以在强湍流，大速度变化下进行流固耦合模拟仿真，并且通过并行优化可以得到实时的较低分辨率仿真，同时基于易扩展性，可以在多节点多GPU的服务器上实现高性能大规模流固耦合的流体模拟仿真。由于兼具精度、稳定性和效率，本发明提供的方法可以广泛应用于流体仿真的设计领域，如航空，汽车，高铁，轮船，潜水艇，建筑，环境等，同时也可以应用于基于计算机物理图形动画的电影，广告特效等领域。本发明与现有技术相比的优点在于：
1.本发明的方法基于NO-CMLBM模型，提出的高精度的高阶松弛参数自适应控制的高阶非正交中心矩松弛模型，可以保证求解器流体低耗散和低色散，提高了传统格子玻尔兹曼模型的稳定性和精确度。相比于现有的流体仿真方法物理更加真实，模型更加稳定、灵活、精确，算法更加高效。
2.本发明的方法在低粘度、湍流很强、大速度变化的情况下可以进行准确，稳定的流体仿真。
3.本发明可以在动态固体边界的条件下，支持固体和流体相互作用，仿真出真实的单向和双向的流固耦合的流体仿真。
4.本发明提出的高性能流体仿真平台可以做到高效准确稳定的大规模流体仿真，为工业产品设计及影视特效提供重要的流固耦合仿真平台。
5.相比现有的流体仿真方法，本发明的方法可以在很复杂的边界情况下进行稳定准确的流体仿真。
6.相比于现有的流体仿真方法，在满足物理真实的同时，模型具有更高的并行度和计算效率。
7.相比现有的流体仿真方法，本发明可以做到较低分辨率下的实时仿真。
附图说明
图1示出本发明中采用的3D格子点结构示意图；
图2示出本发明中高阶弛豫参数和度量方程之间的关系示意图；
图3示出本发明中浸入式方法的示意图；
图4(a)及图4(b)示出本发明中固体边界的采样点算法示意图；
图5示出本发明中低分辨率实时仿真的渲染效果示意图；
图6示出本发明中动态下落平板与烟雾流固单向交互作用的渲染效果示意图；
图7示出本发明中喷射流吹一堆木箱的流固双向耦合作用的渲染效果示意图；
图8示出本发明中流固双向和单向交互作用同时存在场景的渲染效果示意图；
图9示出本发明中静态大规模虚拟高铁风洞试验的渲染效果示意图；
图10示出本发明中动态大规模龙卷风与汽车，牛耦合作用的渲染效果示意图；
图11示出本发明中火箭推进器仿真的渲染效果示意图；
图12示出本发明中与真实绕圆柱实验对比的渲染效果示意图；
图13示出本发明中与真实斜板下落实验对比的速度场渲染效果示意图；
图14示出本发明流程图。
具体实施方式
下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
本发明提供的一种基于统计动力学的高性能大规模流固耦合流体仿真方法所采用的高精度高阶松弛参数自适应控制的高阶非正交中心矩松弛模型，可以保证求解流体低耗散和低色散误差，大大提高了传统格子玻尔兹曼模型的稳定性和精确度。本发明还构造了介观格子玻尔兹曼仿真与宏观纳维-斯托克斯(Navier-Stokes)仿真之间的通用映射关系，基于此映射进行时间和空间连续尺度的自适应计算。本发明建立了强湍流和大速度变化的高性能流固耦合的流体仿真方法。本发明最终能够实现在较低分辨率下实时仿真和高性能大规模流固耦合的流体仿真。
具体而言，本发明包括以下内容：
1.高精度高阶松弛参数的高阶非正交中心矩松弛模型：
在NO-CMLBM的模型基础上，本发明提出的高精度高阶松弛参数的高阶非正交中心矩松弛模型，采用DnQm(n维离散空间和m维速度空间)的离散波尔兹曼格子结构如图1所示。传统的离散格子玻尔兹曼方程如下所示：
fi(x+ciΔt，t+Δt)-fi(x，t)＝Ωi+Fi
式中，fi(x，t)表示时刻t空间x的第i个方向概率分布函数值，ci表示离散格子结构第i个方向的网格速度，Ωi第i个方向的LBM的碰撞操作，Fi力在第i个方向的投影。在求解过程中可以分为迁移(streaming)过程：fi *(x，t+Δt)＝fi(x-ciΔt，t)和碰撞(collision)过程：
以及边界处理三个步骤，式中，τ表示弛豫时间。在碰撞过程中，首先将每个格子的概率分布函数映射到矩空间，其次在矩空间做弛豫操作，最后将矩空间映射到概率分布函数空间：
其中Ω表示碰撞操作，M表示概率分布函数映射到矩空间的矩阵，R是弛豫参数对角矩阵，包含低阶参数和高阶参数，f表示概率分布函数，m表示矩空间的矩。其中mj＝∑i(ci，a-ua)α(ci，b-ub)β(ci，c-uc)γfi，ci，a、ci，b、ci，c表示第i个方向的网格速度的三个分量，ua、ub、uc表示当前位置速度u的三个分量，α、β、γ表示指数,可以取[0,1,2]，fi表示当前位置的第i个方向概率分布函数值。通过6阶埃尔米特多项式展开逼近NS方程，在非正交中心矩松弛模型中，矩空间的平衡态表达为：
力Fi在矩空间的弛豫过程是：
式中，
表示力Fi在弛豫之后的分布函数，R表示弛豫参数对角矩阵。
其中力Fi在矩空间的弛豫空间的表达：
式中，
表示空间i点的力F在矩空间的弛豫空间对应第1，2，3，10，11，12，23，24，25网格方向的值，Fx、Fy、Fz表示力F在三个维度x,y,z方向的分量。
弛豫参数对角矩阵R中ri＝(3ν+0.5)-1是与实际物理粘度v相关。首次提出目标度量方程，ri表示第i个对角元素，i表示第i阶数，i＝{4，...，8}，设计高阶弛豫参数ri(i表示第i阶数，i＝{9，...，26})和优化参数τi的关系：ri＝(3τi+1/2)-1。首先提出NO-CMLBM的优化参数τi和求解耗散和色散误差存在关系，如图2所示。然后根据度量方程：
其中，xk表示空间中k点，t表示真实物理空间的时间，∈(xk，t)表示t时刻空间中k点的度量方程值，ρ表示真实物理空间的密度，u表示真实物理空间的速度，Π表示二阶矩，δt(·)＝(·)t+Δt-(·)t，Δt表示真实物理空间的时间步长，k表示迭代次数的倍数。优化度量方程，得到足够多数据(由密度场，速度场，高阶参数ri组成的时间序列组成)，其中优化过程是枚举所有的τi，得到最小的度量方程值对应的τi，然后求出当前弛豫过程的高阶参数ri。最后基于数据线性回归得到高阶参数的自适应求解方法。具体表达方式是：
其中，sp是LBM空间p点的状态向量，p表示空间p点，sp＝(ρp，||ρpup||，||Пp||，1)，ρp表示LBM空间p点密度，up表示LBM空间p点速度，Πp表示LBM空间p点二阶矩；θ是线性回归得到的参数向量；
表示LBM空间p点的弛豫时间。
2.介观和宏观物理量的通用映射方法：
基于流体的特性，雷诺数在物理空间保持不变的原则，通过宏观基础物理量(速度、密度、空间步长、时间步长等)归一化到介观格子玻尔兹曼方程中，建立起宏观流体物理量和介观流体物理量的通用映射关系。宏观流体和介观流体的通用映射物理量主要包括：空间步长Δx，时间步长Δt，求解空间尺寸大小x，密度ρ，速度u，单位体积力F，粘度v。关系如下：
t表示真实物理空间的时间，
表示LBM空间的参照速度uref，uref表示物理空间的参照速度，Δx表示真实物理空间的空间步长，
表示LBM空间的时间，Δt表示真实物理空间的时间步长，
表示LBM空间的时间步长，x表示真实物理空间的求解空间尺度大小，
表示LBM空间的求解空间尺度大小，
表示LBM空间的空间步长，ρ表示真实物理空间的密度，ρ0表示参照的物理空间的空气密度，
表示LBM空间的密度，u表示真实物理空间的速度，
表示LBM空间的速度，F表示真实物理空间的单位体积力，
表示LBM空间的单位体积力，v表示真实物理空间的粘度，
表示LBM空间的粘度。
3.基于通用映射的时间和空间尺度的自适应仿真计算方法：
基于步骤2的通用映射构造，在时间和空间上构造连续尺度的映射关系，在时间和空间尺度上进行自适应仿真。
时间尺度上基于每个时刻求出在介观尺度下的速度最大值
保持
在[0.18，0.22]之间，通过通用映射中速度映射关系
已知LBM空间速度
和
以及物理速度u求出uref，然后更新介观LBM空间的速度场
分布，根据新的介观观物理量
重构新的介观分布函数fi，继续进行流体仿真。
由于流体的特性，需要满足不同空间尺度之间具有相同数值的雷诺数，根据不同尺度的空间步长的倍数关系，通过通用映射中时间步长关系
推导不同尺度上迭代次数的倍数k，然后基于倍数k在不同尺度上完成时间的同步一致性。追踪固体的运动轨迹，在固体周围重建精细尺度网格，由于流体结构的不连续现象一般会是发生在不同尺度的边界部分，所以在不同尺度的边界部分，增加它们边界的耦合部分，使得流体能够在连续尺度的边界上过渡的更加自然。精细尺度网格和粗糙尺度网格部分区域耦合在一起形成完整的仿真空间。
首先在粗糙尺度网格上仿真，粗糙尺度网格边界的介观分布函数fi部分方向在做迁移操作：fi *(x，t+Δt)＝fi(x-ciΔt，t)的时候需要通过在精细尺度网格内部采用高精度插值得到，插值的精度满足质量和动量守恒。然后在精细尺度网格内部进行仿真，根据迭代倍数k可以得到精细尺度网格需要迭代的次数，如果k不是整数，迭代次数等于k向上取整的整数。由于固体存在精细尺度网格内部，精细尺度网格在每次仿真的时候，需要进行流固耦合仿真计算；精细尺度网格边界的部分fi在迁移操作求解过程为：通过在粗糙尺度网格内部采用高精度插值后得到介观分布函数fi 1；由于精细尺度网格迭代一次相当于粗糙尺度网格的
次，由当前时刻fi 0，下一时刻fi 1，在时间上线性插值得到
最后由于向上取整取迭代次数，在精细尺度网格仿真k次后，需要将精细尺度网格所有点的时间采用线性插值同步到与粗糙尺度网格的同一时刻。
4.基于浸入式边界处理方法的流固耦合仿真方法：
本发明使用的浸入式边界处理方法，如图3所示，首先优化固体网格，固体拉格朗日网格点采样是基于每个三角形网格采样一点或者泊松采样，得到采样点分布，如图4(a)及图4(b)所示。将欧拉流体求解网格和固体拉格朗日网格耦合，满足滑动和非滑动的边界条件，利用最小插值内核
求解低耗散的流固边界力的分布，r表示插值点到插值中心的距离。
非滑动的边界条件是通过速度的耦合边界条件，推导出流固耦合中固体对流体运动作用力的计算方法为：
式中，Xp表示空间中p点的位置，Fs→f(Xp)表示空间p点处固体s对流体f的作用力，ub(Xp)表示固体边界b上p点的速度，u(Xp)表示流体在p点的的速度。然后基于作用力与反作用力得到流体运动对固体运动作用力的求解方法为：Ff→s(Xp)＝-Fs→f(Xp)。
滑动的边界条件推导出流固耦合中固体运动对流体运动作用力的计算方法为：
式中，n表示固体边界b上p点的的法向量。同样的流体运动对固体运动作用力的求解方法为：Ff→s(Xp)＝-Fs→f(Xp)。在固体运动(刚体运动)中所受流体作用的合力以及合转矩是通过Fs＝∑pFf→s(Xp)Δsp，τs＝∑p(Xp-B)×Ff→s(Xp)Δsp计算得到，式中，B表示固体边界b的质心位置，Δsp表示离散积分面元。
5.建立低分变率实时仿真和高性能大规模流固交互流体仿真平台：
本发明利用格子玻尔兹曼方程的并行性，基于统一计算架构(CUDA)和消息传递接口(MPI)并行技术，通过并行优化可以在低分变率实时仿真，如图5所示。在NVIDIA RTX2080Ti显卡上能够做到100×200×100分辨率的实时仿真。由于本发明的可扩展性，优化并行实现的数据结构，可以在多节点多GPU的服务器上实现高性能大规模流固耦合的流体模拟仿真。
Claims (5)
1.一种基于统计动力学的高性能大规模流固耦合流体仿真方法，其特征在于，包括以下步骤：
步骤(1)流体仿真建模：在非正交中心矩驰豫格子玻尔兹曼模型的基础上，基于六阶埃尔米特多项式近似N-S方程的非正交中心矩驰豫格子玻尔兹曼模型提出目标度量方程
xk表示空间中k点，t表示LBM空间的时间，∈(xk，t)表示t时刻空间中k点的度量方程值，ρ表示LBM空间的密度，u表示LBM空间的速度，Π表示LBM空间的二阶矩，δt(·)＝(·)t+Δt-(·)t，Δt表示真实物理空间的时间步长，k表示迭代次数的倍数，在多种不同仿真环境下优化目标度量方程得到足够多的训练数据，然后基于训练数据用线性回归的方法去自适应控制非正交中心矩驰豫格子玻尔兹曼模型的高阶参数；
步骤(2)通用映射构造：根据流体的性质，雷诺数在物理空间下保持不变，基于宏观基础物理量归一化到介观格子玻尔兹曼方程中，建立宏观流体物理量和介观流体物理量的通用映射关系，其中，宏观流体和介观流体的通用映射物理量包括空间步长、时间步长、求解空间尺度大小、密度、速度、单位体积力、粘度；
步骤(3)时间空间连续尺度自适应仿真：基于步骤(2)的通用映射关系，在时间和空间连续尺度上进行自适应仿真，构造连续尺度之间的映射关系；
步骤(4)高性能流固耦合仿真：使用浸入式边界处理方法，将流体方程求解的欧拉网格和固体运动求解的拉格朗日网格耦合，满足滑动和非滑动的边界条件，利用最小插值内核构建低耗散的流固边界，其中，耦合力方法基于步骤(2)的通用映射关系，计算介观流体和宏观刚体之间的耦合力的映射关系并且提供单向流固耦合和双向流固耦合流体仿真；
步骤(5)建立低分变率实时仿真和高性能大规模流固交互流体仿真平台：在NVIDIARTX 2080Ti显卡上可以做到100×200×100分辨率的实时仿真，并且其易扩展可以在多节点多GPU上高效并行化实现。
2.如权利要求1所述的一种基于统计动力学的高性能大规模流固耦合流体仿真方法，其特征在于，步骤(2)中，所述通用映射关系为：
式中，t表示真实物理空间的时间，
表示LBM空间的参照速度，uref表示物理空间的参照速度，Δx表示真实物理空间的空间步长，
表示LBM空间的时间，Δt表示真实物理空间的时间步长，
表示LBM空间的时间步长，x表示真实物理空间的求解空间尺度大小，
表示LBM空间的求解空间尺度大小，
表示LBM空间的空间步长，ρ表示真实物理空间的密度，ρ0表示参照的物理空间的空气密度，
表示LBM空间的密度，u表示真实物理空间的速度，
表示LBM空间的速度，F表示真实物理空间的单位体积力，
表示LBM空间的单位体积力，ν表示真实物理空间的粘度，
表示LBM空间的粘度。
3.如权利要求2所述的一种基于统计动力学的高性能大规模流固耦合流体仿真方法，其特征在于，步骤(3)中，在不同尺度的流体结构的边界部分，增加边界的耦合部分，使得流体能够在连续尺度的边界上过渡地更加自然，同时保证空间连续尺度之间的雷诺数一致，根据不同尺度的空间步长的倍数关系，通过通用映射关系中时间步长关系
推导不同尺度上迭代次数的倍数k，然后基于倍数k在不同尺度上完成时间同步一致；
时间尺度上基于每个时刻求出在介观尺度下的速度最大值
保持
在[0.18,0.22]之间，通过通用映射关系中速度映射关系
已知LBM空间的速度
和
以及真实物理空间的速度u求出uref，然后更新介观LBM空间的
分布，根据新的介观物理量
重构介观分布函数fi，继续进行流体仿真。
4.如权利要求1所述的一种基于统计动力学的高性能大规模流固耦合流体仿真方法，其特征在于，步骤(4)中，所述浸入式边界处理方法包括以下步骤：
首先优化固体网格，基于每个三角形网格采样一点或者泊松采样构造固体拉格朗日网格点集分布，随后将欧拉流体求解网格和固体拉格朗日网格耦合，满足滑动和非滑动的边界条件，利用最小插值内核
求解低耗散的流固边界力的分布，其中，r表示插值点到插值中心的距离。
5.如权利要求1所述的一种基于统计动力学的高性能大规模流固耦合流体仿真方法，其特征在于，步骤(5)中，利用玻尔兹曼方程的并行性，通过并行优化可以在NVIDIA RTX2080Ti显卡上做到100×200×100分辨率的实时仿真；基于统一计算架构和消息传递接口并行技术，优化并行实现的数据结构，可以在多节点多GPU的服务器上实现高性能大规模流固耦合的流体仿真。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010488622.9A CN111695309B (zh) | 2020-06-02 | 2020-06-02 | 基于统计动力学的高性能大规模流固耦合流体仿真方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010488622.9A CN111695309B (zh) | 2020-06-02 | 2020-06-02 | 基于统计动力学的高性能大规模流固耦合流体仿真方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111695309A true CN111695309A (zh) | 2020-09-22 |
| CN111695309B CN111695309B (zh) | 2023-03-21 |
Family
ID=72479202
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010488622.9A Active CN111695309B (zh) | 2020-06-02 | 2020-06-02 | 基于统计动力学的高性能大规模流固耦合流体仿真方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111695309B (zh) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112861374A (zh) * | 2021-03-05 | 2021-05-28 | 深圳泽森软件技术有限责任公司 | 基于预控制器的多物理耦合仿真处理方法、装置和设备 |
| CN113420364A (zh) * | 2021-01-25 | 2021-09-21 | 中国第一汽车股份有限公司 | 一种基于流固耦合的电泳过程车身结构变形仿真方法 |
| CN113505518A (zh) * | 2021-06-30 | 2021-10-15 | 同济大学 | 用于质子交换膜燃料电池催化剂浆料制备过程的模拟方法 |
| CN113688487A (zh) * | 2021-08-16 | 2021-11-23 | 苏州同元软控信息技术有限公司 | 乘用车二维空调管路流体动力学的仿真方法及装置 |
| CN113935169A (zh) * | 2021-10-14 | 2022-01-14 | 深圳泽森软件技术有限责任公司 | 物理仿真方法、装置、计算机设备和存储介质 |
| CN115526091A (zh) * | 2022-11-22 | 2022-12-27 | 中国人民解放军国防科技大学 | 一种面向多物理场应用的分离式耦合数值模拟方法和装置 |
| CN115630559A (zh) * | 2022-12-08 | 2023-01-20 | 中国人民解放军国防科技大学 | 一种基于粒子网格适配算法的流固耦合方法以及装置 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104268943A (zh) * | 2014-09-28 | 2015-01-07 | 北京航空航天大学 | 一种基于欧拉-拉格朗日耦合方法的流体仿真方法 |
| WO2017084106A1 (zh) * | 2015-11-20 | 2017-05-26 | 田川 | 一种数值模拟飞行器流场的系统及方法 |
| CN109376463A (zh) * | 2018-11-16 | 2019-02-22 | 重庆科技学院 | 一种尾流下弹性支撑圆柱驰振流固耦合分析方法 |
| CN110110802A (zh) * | 2019-05-14 | 2019-08-09 | 南京林业大学 | 基于高阶条件随机场的机载激光点云分类方法 |
-
2020
- 2020-06-02 CN CN202010488622.9A patent/CN111695309B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104268943A (zh) * | 2014-09-28 | 2015-01-07 | 北京航空航天大学 | 一种基于欧拉-拉格朗日耦合方法的流体仿真方法 |
| WO2017084106A1 (zh) * | 2015-11-20 | 2017-05-26 | 田川 | 一种数值模拟飞行器流场的系统及方法 |
| CN109376463A (zh) * | 2018-11-16 | 2019-02-22 | 重庆科技学院 | 一种尾流下弹性支撑圆柱驰振流固耦合分析方法 |
| CN110110802A (zh) * | 2019-05-14 | 2019-08-09 | 南京林业大学 | 基于高阶条件随机场的机载激光点云分类方法 |
Non-Patent Citations (3)
| Title |
|---|
| FARSHAD GHARIBI等: "Simulation of high-viscosity-ratio multicomponent fluid flow using a pseudopotential model based on the nonorthogonal central-moments lattice Boltzmann method", 《PHYSICAL REVIEW E》 * |
| WEI LI等: "Continuous-Scale Kinetic Fluid Simulation", 《IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS》 * |
| 刘克同等: "基于MRT-LBM大涡模拟的桥梁气动参数数值仿真", 《四川大学学报(工程科学版)》 * |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113420364A (zh) * | 2021-01-25 | 2021-09-21 | 中国第一汽车股份有限公司 | 一种基于流固耦合的电泳过程车身结构变形仿真方法 |
| CN113420364B (zh) * | 2021-01-25 | 2022-07-22 | 中国第一汽车股份有限公司 | 一种基于流固耦合的电泳过程车身结构变形仿真方法 |
| CN112861374A (zh) * | 2021-03-05 | 2021-05-28 | 深圳泽森软件技术有限责任公司 | 基于预控制器的多物理耦合仿真处理方法、装置和设备 |
| CN113505518A (zh) * | 2021-06-30 | 2021-10-15 | 同济大学 | 用于质子交换膜燃料电池催化剂浆料制备过程的模拟方法 |
| CN113688487A (zh) * | 2021-08-16 | 2021-11-23 | 苏州同元软控信息技术有限公司 | 乘用车二维空调管路流体动力学的仿真方法及装置 |
| CN113688487B (zh) * | 2021-08-16 | 2024-01-02 | 苏州同元软控信息技术有限公司 | 乘用车二维空调管路流体动力学的仿真方法及装置 |
| CN113935169A (zh) * | 2021-10-14 | 2022-01-14 | 深圳泽森软件技术有限责任公司 | 物理仿真方法、装置、计算机设备和存储介质 |
| CN115526091A (zh) * | 2022-11-22 | 2022-12-27 | 中国人民解放军国防科技大学 | 一种面向多物理场应用的分离式耦合数值模拟方法和装置 |
| CN115630559A (zh) * | 2022-12-08 | 2023-01-20 | 中国人民解放军国防科技大学 | 一种基于粒子网格适配算法的流固耦合方法以及装置 |
| CN115630559B (zh) * | 2022-12-08 | 2023-03-10 | 中国人民解放军国防科技大学 | 一种基于粒子网格适配算法的流固耦合方法以及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111695309B (zh) | 2023-03-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111695309B (zh) | 基于统计动力学的高性能大规模流固耦合流体仿真方法 | |
| Stam | Real-time fluid dynamics for games | |
| CN104268943A (zh) | 一种基于欧拉-拉格朗日耦合方法的流体仿真方法 | |
| Wandel et al. | Teaching the incompressible Navier–Stokes equations to fast neural surrogate models in three dimensions | |
| US20070239414A1 (en) | Method of simulating detailed movements of fluids using derivative particles | |
| US9984489B2 (en) | Fluid dynamics framework for animated special effects | |
| CN110717269B (zh) | 一种基于网格和粒子耦合的流体表面细节保护方法 | |
| Quigley et al. | Real-time interactive tree animation | |
| CN105278346A (zh) | 一种基于离散格子Boltzmann双分布模型的热流体仿真方法 | |
| Weißmann et al. | Underwater rigid body dynamics | |
| Wandel et al. | Unsupervised deep learning of incompressible fluid dynamics | |
| CN103426196A (zh) | 一种流体环境下的关节动画建模技术 | |
| Zong et al. | Neural stress fields for reduced-order elastoplasticity and fracture | |
| Gao et al. | An efficient FLIP and shape matching coupled method for fluid–solid and two-phase fluid simulations | |
| Gao et al. | Accelerating liquid simulation with an improved data‐driven method | |
| CN104463934B (zh) | 一种“质点‑弹簧”系统驱动的点集模型动画自动生成方法 | |
| Song et al. | Derivative particles for simulating detailed movements of fluids | |
| CN110941894A (zh) | 一种基于超级计算机的大规模流体仿真的原位可视化方法 | |
| CN116933674A (zh) | 基于lbm的快速流体模拟方法 | |
| Movania et al. | A Novel GPU‐Based Deformation Pipeline | |
| Wu et al. | Improved divergence‐free smoothed particle hydrodynamics via priority of divergence‐free solver and SOR | |
| Vyatkin et al. | Optimized finite element method using free-form volume patches for deformation of three-dimensional objects | |
| Li et al. | Real-time physically plausible simulation of forest | |
| Liu et al. | Ellipsoidal-blob approximation of 3D models and its applications | |
| Sewall et al. | Fast Fluid Simulation Using Residual Distribution Schemes. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |