CN111859529A - A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft - Google Patents

A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft Download PDF

Info

Publication number
CN111859529A
CN111859529A CN202010528125.7A CN202010528125A CN111859529A CN 111859529 A CN111859529 A CN 111859529A CN 202010528125 A CN202010528125 A CN 202010528125A CN 111859529 A CN111859529 A CN 111859529A
Authority
CN
China
Prior art keywords
grid
domain
boundary
dynamic calculation
calculation domain
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
Application number
CN202010528125.7A
Other languages
Chinese (zh)
Other versions
CN111859529B (en
Inventor
蒋崇文
胡姝瑶
高振勋
李椿萱
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202010528125.7A priority Critical patent/CN111859529B/en
Publication of CN111859529A publication Critical patent/CN111859529A/en
Application granted granted Critical
Publication of CN111859529B publication Critical patent/CN111859529B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • G06F30/15Vehicle, aircraft or watercraft design
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法,针对现有多重网格方法在模拟飞行器绕流时存在大量无效计算的问题,通过在细网格和粗网格上采用对流、粘性两类动态计算域,使得流动控制方程在每个迭代步中均在仅包含未收敛受扰单元的对流动态计算域中求解,而粘性效应也仅在只包含粘性效应主导单元的粘性动态计算域内考虑,从而实现了仅对未收敛受扰单元求解、仅在局部区域考虑粘性效应的求解思路。通过有效减少单个迭代步的计算量和达到收敛态所需的迭代步数,本发明可显著提高多重网格方法模拟飞行器绕流的计算效率。

Figure 202010528125

The invention discloses a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft. Aiming at the problem that a large number of invalid calculations exist in the existing multi-grid method when simulating the flow around an aircraft, by adjusting the fine grid and coarse grid Convective and viscous dynamic computational domains are adopted in the above, so that the flow governing equations are solved in the convective dynamic computational domain containing only the unconverged disturbed elements in each iteration step, and the viscous effect is only contained in the viscous effect dominant element. The viscous dynamic calculation domain of , realizes the solution idea of only solving the unconverged disturbed element and only considering the viscous effect in the local area. By effectively reducing the calculation amount of a single iteration step and the number of iteration steps required to reach a convergent state, the invention can significantly improve the calculation efficiency of the multi-grid method for simulating the flow around the aircraft.

Figure 202010528125

Description

一种飞行器绕流数值模拟的多重网格扰动域更新加速方法A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft

技术领域technical field

本发明涉及计算流体力学技术领域,尤其涉及一种飞行器绕流数值模拟的多重网格扰动域更新加速方法。The invention relates to the technical field of computational fluid dynamics, in particular to a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft.

背景技术Background technique

随着计算机技术的飞速发展,计算流体力学数值模拟已成为当代飞行器气动设计的重要技术支撑之一。尽管数值模拟能够精准预测飞行器的静、动态气动特性,但对于飞行器气动设计而言,仍存在数值模拟计算效率难以满足日益增长的实际应用需求的问题。因此,提出适用于飞行器绕流数值模拟的高效方法对提高飞行器设计迭代效率、缩短研制周期具有重要意义和工程应用价值。With the rapid development of computer technology, computational fluid dynamics numerical simulation has become one of the important technical supports for the aerodynamic design of contemporary aircraft. Although numerical simulation can accurately predict the static and dynamic aerodynamic characteristics of aircraft, for the aerodynamic design of aircraft, there is still a problem that the computational efficiency of numerical simulation cannot meet the increasing demand for practical applications. Therefore, it is of great significance and engineering application value to propose an efficient method for numerical simulation of flow around aircraft to improve the efficiency of aircraft design iteration and shorten the development cycle.

多重网格方法是飞行器绕流数值模拟的主流加速技术之一。该方法既可用于飞行器静态气动特性预测中定常绕流流场的迭代求解,也适用于飞行器动态气动特性预测中非定常绕流流场迭代推进中的内迭代求解。多重网格方法是一种在粗、细网格间循环求解流动控制方程的加速收敛技术,其加速原理有两方面:其一,粗网格能以更大的迭代步长推进求解,从而可以带动细网格更快达到收敛态;其二,粗、细网格间的循环求解有助于衰减阻碍细网格收敛的高频误差,从而可以使细网格的整体误差能更快地衰减。Multigrid method is one of the mainstream acceleration techniques for numerical simulation of flow around aircraft. The method can be used not only for the iterative solution of the steady surrounding flow field in the prediction of the static aerodynamic characteristics of the aircraft, but also for the inner iterative solution in the iterative propulsion of the unsteady surrounding flow field in the prediction of the dynamic aerodynamic characteristics of the aircraft. The multi-grid method is an accelerated convergence technique that cyclically solves the flow control equations between coarse and fine grids. The acceleration principle has two aspects: First, the coarse grid can advance the solution with a larger iterative step size, so that it can be solved. Drive the fine mesh to reach the convergence state faster; secondly, the cyclic solution between the coarse and fine meshes helps to attenuate the high-frequency error that hinders the convergence of the fine mesh, so that the overall error of the fine mesh can be attenuated faster .

目前,多重网格方法针对飞行器绕流数值模拟时仍存在一个弊端,即每个迭代步均必须对飞行器流场网格预设计算域中的所有网格单元更新求解。这种全局更新的求解方式忽略了求解过程中飞行器绕流流场演化的特点,从而会造成不同程度的无效计算。对于亚声速飞行器,流动的数学性质为椭圆型,因此,为避免数值模拟的有限边界所导致的非物理响应,飞行器流场网格的预设计算域需要取数十倍于飞行器参考长度的范围进行模拟。然而,随着飞行器绕流流场逐渐趋近于定常态,不满足流动控制方程的区域也相应收缩,飞行器流场网格的预设计算域也将趋于冗余。对于超声速飞行器,流动的数学性质为双曲型,即流场中任一点仅受其上游依赖域内流动的影响。因此,在迭代过程中,位于飞行器绕流流场下游的区域不会先于其上游流场收敛,已达到定常态的上游流动也不会受到未收敛下游流动更新的影响。由此可知,在整个迭代求解过程中,始终对飞行器流场网格预设计算域内的所有网格单元全局更新求解将产生大量无效计算,从而显著降低飞行器绕流数值模拟的计算效率。因此,若能消除现有方法中的无效计算,必将显著提升飞行器绕流数值模拟的计算效率,从而有效解决飞行器气动设计中数值模拟计算效率难以满足日益增长的实际应用需求的问题。At present, the multi-grid method still has a drawback in the numerical simulation of the flow around the aircraft, that is, each iteration step must update and solve all the grid elements in the preset calculation domain of the aircraft flow field grid. This global update solution method ignores the characteristics of the evolution of the flow field around the aircraft during the solution process, resulting in different degrees of invalid calculation. For a subsonic aircraft, the mathematical properties of the flow are elliptical. Therefore, in order to avoid the unphysical response caused by the finite boundary of the numerical simulation, the preset computational domain of the aircraft flow field grid needs to be in the range of dozens of times the reference length of the aircraft. Do a simulation. However, as the flow field around the aircraft gradually approaches the steady state, the region that does not satisfy the flow control equations also shrinks accordingly, and the preset computational domain of the aircraft flow field grid will also tend to be redundant. For supersonic vehicles, the mathematical properties of the flow are hyperbolic, that is, any point in the flow field is only affected by the flow in its upstream dependent domain. Therefore, in the iterative process, the region located downstream of the flow field around the aircraft will not converge before its upstream flow field, and the upstream flow that has reached a steady state will not be affected by the update of the unconverged downstream flow. It can be seen that in the entire iterative solution process, the global update and solution of all grid elements in the preset calculation domain of the aircraft flow field grid will generate a large number of invalid calculations, which will significantly reduce the computational efficiency of the numerical simulation of the flow around the aircraft. Therefore, if the invalid calculation in the existing method can be eliminated, the computational efficiency of the numerical simulation of the flow around the aircraft will be significantly improved, thereby effectively solving the problem that the computational efficiency of numerical simulation in the aerodynamic design of the aircraft cannot meet the increasing practical application requirements.

发明内容SUMMARY OF THE INVENTION

有鉴于此,本发明提供了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法,用以解决现有多重网格方法在模拟飞行器绕流时存在大量无效计算的问题。In view of this, the present invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft, so as to solve the problem of a large number of ineffective calculations when simulating flow around an aircraft in the existing multi-grid method.

本发明提供的一种飞行器绕流数值模拟的多重网格扰动域更新加速方法,其特征在于,包括如下步骤:The invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft, which is characterized in that it includes the following steps:

S1:读入数据,包括飞行器流场的细网格、预设计算域、边界条件和计算设置;S1: Read in data, including the fine mesh of the aircraft flow field, the preset computational domain, boundary conditions and calculation settings;

S2:通过稀疏所述细网格,生成指定层数的粗网格;S2: Generate a coarse grid with a specified number of layers by thinning the fine grid;

S3:根据来流条件或根据给定流场,对所述粗网格和所述细网格上预设计算域中的流场进行初始化;S3: initialize the flow field in the preset computational domain on the coarse grid and the fine grid according to the incoming flow condition or according to a given flow field;

S4:根据流场初始化方式,建立对流动态计算域和粘性动态计算域;S4: According to the flow field initialization method, establish the convective dynamic calculation domain and the viscous dynamic calculation domain;

S5:将所有层网格求解所需存储的数据分为两类;第一类数据是单元的固有信息,采用静态数据结构存储所述预设计算域中所有网格单元的信息,包括网格坐标和流场变量;第二类数据是与求解更新相关的信息,采用动态数据结构仅存储所述对流动态计算域内网格单元的信息,包括守恒量更新量和当地迭代步长;S5: Divide the data to be stored for all layers of grid solution into two types; the first type of data is the inherent information of the unit, and the static data structure is used to store the information of all grid units in the preset computational domain, including the grid Coordinates and flow field variables; the second type of data is the information related to the solution update, and the dynamic data structure is used to store only the information of the grid cells in the convective dynamic calculation domain, including the update amount of the conservation quantity and the local iteration step size;

S6:根据边界条件类型,为当前层网格的边界虚网格的守恒量赋值;S6: According to the boundary condition type, assign value to the conservation quantity of the boundary virtual grid of the current layer grid;

S7:将流动控制方程的残差项分为无粘项与粘性项两类;其中,所述无粘项包括对流通量和粗网格的强迫函数,所述无粘项在对流动态计算域中计算;所述粘性项包括粘性通量和湍流模型方程源项,所述粘性项在粘性动态计算域中计算;S7: Divide the residual term of the flow control equation into two types: inviscid term and viscous term; wherein, the inviscid term includes the forcing function of convective flux and coarse grid, and the inviscid term is in the convective dynamic computational domain Calculated in ; the viscous term includes viscous flux and turbulence model equation source term, and the viscous term is calculated in the viscous dynamic computational domain;

S8:在所述对流动态计算域中,求解守恒量更新量,更新流场变量;S8: in the convective dynamic calculation domain, solve the update amount of the conserved quantity, and update the flow field variable;

S9:判断当前网格层是否为细网格;若是,则执行步骤S10;若否,则执行步骤S21;S9: determine whether the current grid layer is a fine grid; if yes, go to step S10; if not, go to step S21;

S10:判断细网格是否达到收敛条件;若否,则执行步骤S11;若是,则执行步骤S23;S10: Determine whether the fine mesh reaches the convergence condition; if not, execute step S11; if yes, execute step S23;

S11:对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动;若是,则执行步骤S12;若否,则返回步骤S11,对所述对流动态计算域的下一边界单元进行判断;遍历所述对流动态计算域的所有边界单元后,执行步骤S13;S11: Determine whether all the boundary elements of the convective dynamic calculation domain have been subjected to inviscid disturbance one by one; if so, go to step S12; Judging; after traversing all boundary elements of the convective dynamic computational domain, perform step S13;

S12:衡量无粘扰动的传播方向,并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域,返回步骤S11,对所述对流动态计算域的下一边界单元进行判断;S12: Measure the propagation direction of the inviscid perturbation, and include the unit located in the propagation direction in the adjacent units of the boundary unit into the convective dynamic calculation domain, and return to step S11 to perform the next boundary unit of the convective dynamic calculation domain. judge;

S13:对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件;若是,则执行步骤S14;若否,则返回步骤S13,对所述对流动态计算域的下一边界单元进行判断;遍历所述对流动态计算域的所有边界单元后,执行步骤S16;S13: Judging, one by one, whether all boundary elements of the convective dynamic calculation domain satisfy the four conditions of converged, located in compressible flow, located at the most upstream and no longer affected by other units in the convective dynamic calculation domain; , then go to step S14; if not, go back to step S13 to judge the next boundary unit of the convective dynamic computational domain; after traversing all the boundary units of the convective dynamic computational domain, go to step S16;

S14:将该边界单元从所述对流动态计算域中移除,并判断该边界单元是否存在于所述粘性动态计算域中;若是,则执行步骤S15;若否,则返回步骤S13,对所述对流动态计算域的下一边界单元进行判断;S14: Remove the boundary element from the convective dynamic calculation domain, and determine whether the boundary element exists in the viscous dynamic calculation domain; if so, execute step S15; The next boundary element of the convective dynamic computational domain is judged as described above;

S15:将该边界单元从所述粘性动态计算域中移除,返回步骤S13,对所述对流动态计算域的下一边界单元进行判断;S15: remove the boundary element from the viscous dynamic calculation domain, and return to step S13 to judge the next boundary element of the convective dynamic calculation domain;

S16:对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导;若是,则执行步骤S17;若否,则返回步骤S16,对所述粘性动态计算域的下一边界单元进行判断;遍历所述粘性动态计算域的所有边界单元后,执行步骤S18;S16: Determine whether all boundary elements of the viscous dynamic calculation domain are dominated by the viscous effect one by one; if so, go to step S17; if not, return to step S16 to judge the next boundary element of the viscous dynamic calculation domain ; After traversing all boundary elements of the viscous dynamic computational domain, execute step S18;

S17:将该边界单元的紧邻单元中位于所述对流动态计算域中的单元纳入所述粘性动态计算域,返回步骤S16,对所述粘性动态计算域的下一边界单元进行判断;S17: Include the unit located in the convective dynamic calculation domain in the adjacent units of the boundary unit into the viscous dynamic calculation domain, and return to step S16 to judge the next boundary element of the viscous dynamic calculation domain;

S18:根据细网格的对流动态计算域和粘性动态计算域,更新所有层粗网格的对流动态计算域和粘性动态计算域;S18: According to the convective dynamic calculation domain and the viscous dynamic calculation domain of the fine grid, update the convective dynamic calculation domain and the viscous dynamic calculation domain of the coarse grid of all layers;

S19:根据所述对流动态计算域,重新分配所有粗网格和细网格中第二类数据的存储空间;S19: According to the convective dynamic calculation domain, reallocate the storage space of the second type of data in all coarse grids and fine grids;

S20:将守恒量插值至更粗网格,返回步骤S6,进入更粗网格的计算;S20: Interpolate the conserved amount to a coarser grid, return to step S6, and enter the calculation of the coarser grid;

S21:判断是否需要计算更粗网格;若是,则返回步骤S20;若否,则执行步骤S22;S21: Determine whether a coarser grid needs to be calculated; if yes, go back to step S20; if not, go to step S22;

S22:将守恒量修正量插值回更细网格,并判断更细网格是否为细网格;若是,则返回步骤S6,进入下一迭代步的计算;若否,则返回步骤S6,进入更细网格的计算;S22: Interpolate the correction amount of the conservation amount back to a finer grid, and determine whether the finer grid is a fine grid; if so, return to step S6, and enter the calculation of the next iteration step; if not, return to step S6, enter Calculation of finer grids;

S23:输出结果。S23: output the result.

在一种可能的实现方式中,在本发明提供的上述多重网格扰动域更新加速方法中,步骤S3,根据来流条件或根据给定流场,对所述粗网格和所述细网格上预设计算域中的流场进行初始化,具体包括:In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, in step S3, according to the incoming flow conditions or according to a given flow field, the coarse grid and the fine grid are Initialize the flow field in the preset computational domain on the grid, including:

根据来流条件初始化,将所述粗网格和所述细网格上预设计算域中所有网格单元的守恒量赋为来流值;Initializing according to the incoming flow condition, assigning the conservation quantities of all grid cells in the preset computational domain on the coarse grid and the fine grid as the incoming flow value;

根据给定流场初始化,将所述细网格上预设计算域中所有网格单元的守恒量赋为给定流场值,所述粗网格上预设计算域中所有网格单元的守恒量根据细网格上所有网格单元的守恒量插值获得。According to the initialization of the given flow field, the conservation quantities of all grid cells in the preset calculation domain on the fine grid are assigned as the given flow field values, and the values of all grid cells in the preset calculation domain on the coarse grid The conserved quantities are obtained by interpolating the conserved quantities of all grid cells on the fine grid.

在一种可能的实现方式中,在本发明提供的上述多重网格扰动域更新加速方法中,步骤S4,根据流场初始化方式,建立对流动态计算域和粘性动态计算域,具体包括:In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, step S4, according to the flow field initialization method, establishes a convective dynamic calculation domain and a viscous dynamic calculation domain, which specifically includes:

根据来流条件初始化,对于对流动态计算域,最粗层网格取紧邻壁面的1层单元作为初始单元,细网格取壁面的L层相邻单元作为初始单元,L满足:According to the initialization of the inflow conditions, for the convective dynamic calculation domain, the coarsest grid layer takes the 1-layer element next to the wall as the initial element, and the fine grid takes the L-layer adjacent element of the wall as the initial element, and L satisfies:

L=2l (1)L=2 l (1)

其中,l表示粗网格层数,l为正整数,由步骤S1读入的计算设置给定;较粗层网格的初始单元根据细网格的初始单元生成;对于粘性动态计算域,所有网格层均取紧邻壁面的1层单元作为初始单元;Among them, l represents the number of coarse grid layers, and l is a positive integer, which is given by the calculation settings read in step S1; the initial elements of the coarser grids are generated according to the initial elements of the fine grid; for the viscous dynamic computational domain, all The grid layer takes the 1-layer unit next to the wall as the initial unit;

根据给定流场初始化,对于对流动态计算域,细网格的初始单元为给定流场中流动特性与来流条件不符的已受扰单元,满足:According to the initialization of the given flow field, for the convective dynamic computational domain, the initial cell of the fine grid is the disturbed cell whose flow characteristics in the given flow field are inconsistent with the incoming flow conditions, satisfying:

||W-W||/||ΔW(1)||max>εa,c (2)||WW ||/||ΔW (1) || max >ε a,c (2)

其中,W表示守恒量;W表示来流条件的守恒量;||ΔW(1)||max表示对流动Among them, W represents the conserved quantity; W represents the conserved quantity of the incoming flow condition; ||ΔW (1) || max represents the opposite flow

态计算域所有单元的守恒量更新量在第1迭代步的最大值;εa,c表示对流新增阈值,取10-4≤εa,c≤10-6;对于粘性动态计算域,细网格的初始单元为细网格对流动态计算域中的粘性效应主导单元,满足:is the maximum value of the update amount of the conservation quantity of all units in the viscous computational domain in the first iteration step; ε a,c represent the new threshold for convection, which is 10 -4 ≤ε a,c ≤10 -6 ; for the viscous dynamic computational domain, fine The initial element of the grid is the viscous effect dominant element in the fine-grid convection dynamic computational domain, which satisfies:

Figure BDA0002534371280000051
Figure BDA0002534371280000051

Figure BDA0002534371280000052
Figure BDA0002534371280000052

其中,Ψ表示粘性效应衡量参数,为粘性扰动与无粘扰动的质量流量之比;

Figure BDA0002534371280000053
表示第1迭代步的粘性效应衡量参数;I,J,K分别表示网格方向;
Figure BDA0002534371280000054
表示对流通量Jacobian矩阵沿i方向的谱半径,
Figure BDA0002534371280000055
表示粘性通量Jacobian矩阵沿i方向的谱半径;εa,v表示粘性新增阈值,取10-3≤εa,v≤10-2;所有粗网格层的对流动态计算域依据细网格的对流动态计算域生成,所有粗网格层的粘性动态计算域依据细网格的粘性动态计算域生成。Among them, Ψ represents the viscous effect measurement parameter, which is the ratio of viscous disturbance to mass flow without viscous disturbance;
Figure BDA0002534371280000053
Represents the viscous effect measurement parameter of the first iteration step; I, J, K represent the grid direction respectively;
Figure BDA0002534371280000054
represents the spectral radius of the Jacobian matrix of convective flux along the i direction,
Figure BDA0002534371280000055
Represents the spectral radius of the viscous flux Jacobian matrix along the i direction; ε a,v represents the new viscosity threshold, which is 10 -3 ≤ε a,v ≤10 -2 ; the convection dynamic computational domain of all coarse mesh layers is based on the fine mesh The convective dynamic computational domain of the grid is generated, and the viscous dynamic computational domain of all coarse grid layers is generated according to the viscous dynamic computational domain of the fine grid.

在一种可能的实现方式中,在本发明提供的上述多重网格扰动域更新加速方法中,步骤S7,将流动控制方程的残差项分为无粘项与粘性项两类;其中,所述无粘项包括对流通量和粗网格的强迫函数,所述无粘项在对流动态计算域中计算;所述粘性项包括粘性通量和湍流模型方程源项,所述粘性项在粘性动态计算域中计算,具体包括:In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, in step S7, the residual term of the flow control equation is divided into two types: inviscid term and viscous term; The inviscid term includes the convective flux and the forcing function of the coarse grid, and the inviscid term is calculated in the convective dynamic computational domain; the viscous term includes the viscous flux and the turbulence model equation source term, and the viscous term is in the viscous flow. Computation in the dynamic computational domain, including:

流动控制方程表示为:The flow control equation is expressed as:

Figure BDA0002534371280000061
Figure BDA0002534371280000061

其中,t表示时间,Fc表示对流通量,Fv表示粘性通量,QF表示粗网格的强迫函数,QT表示湍流模型方程源项,|Ω|表示网格单元的体积,Nf表示单元面的面数,ΔS表示单元面的面积;式(5)的右端项统称为流动控制方程的残差项,

Figure BDA0002534371280000062
为无粘项,
Figure BDA0002534371280000063
为粘性项;在对流动态计算域中求解残差项的无粘项,在粘性动态计算域中求解残差项的粘性项。where t is the time, F c is the convective flux, F v is the viscous flux, Q F is the forcing function of the coarse mesh, Q T is the source term of the turbulence model equation, |Ω| is the volume of the mesh element, N f represents the number of element faces, and ΔS represents the area of the element face; the right-hand term of equation (5) is collectively referred to as the residual term of the flow control equation,
Figure BDA0002534371280000062
for non-sticky terms,
Figure BDA0002534371280000063
is the viscous term; the inviscid term of the residual term is solved in the convective dynamic computational domain, and the viscous term of the residual term is solved in the viscous dynamic computational domain.

在一种可能的实现方式中,在本发明提供的上述多重网格扰动域更新加速方法中,步骤S11,对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动,具体包括:In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, step S11 is to judge whether all boundary elements of the convective dynamic calculation domain have been subjected to inviscid disturbance one by one, specifically including :

对所述对流动态计算域的所有边界单元逐个判断是否满足:Judge whether all boundary elements of the convective dynamic computational domain satisfy one by one:

||ΔW||/||ΔW(1)||max>εa,c (6)||ΔW||/||ΔW (1) || max >ε a,c (6)

其中,||ΔW||表示对流动态计算域的边界单元在当前迭代步中守恒量更新量的模值;Among them, ||ΔW|| represents the modulo value of the update quantity of the conservation quantity of the boundary element of the convective dynamic computational domain in the current iteration step;

步骤S12,衡量无粘扰动的传播方向,并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域,具体包括:In step S12, the propagation direction of the inviscid disturbance is measured, and the unit located in the propagation direction among the adjacent units of the boundary unit is included in the convective dynamic calculation domain, which specifically includes:

无粘扰动相对于流动以声速传播,沿无粘扰动的传播方向,无粘扰动的传播速度为正;令q表示单位方向向量,则无粘扰动沿q方向传播表示为:The inviscid disturbance propagates at the speed of sound relative to the flow, and along the propagation direction of the inviscid disturbance, the propagation velocity of the inviscid disturbance is positive; let q represent the unit direction vector, then the inviscid disturbance propagates along the q direction as:

u·q+a>0 (7)u·q+a>0 (7)

其中,u表示流动速度矢量,a表示声速;q取对流动态计算域边界单元的格心指向边界单元某一格点的单位向量,若q满足式(7),则将与对流动态计算域边界单元共享该格点的紧邻单元纳入对流动态计算域。Among them, u represents the flow velocity vector, a represents the speed of sound; q takes the unit vector of the grid center of the boundary element of the convective dynamic calculation domain pointing to a certain grid point of the boundary element, if q satisfies the formula (7), it will be combined with the boundary of the convective dynamic calculation domain. Cells that share this grid are included in the convective dynamic computational domain.

在一种可能的实现方式中,在本发明提供的上述多重网格扰动域更新加速方法中,步骤S13,对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件;若是,执行步骤S14;若否,则返回步骤S13,对所述对流动态计算域的下一边界单元进行判断;具体包括如下步骤:In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, in step S13, it is judged one by one for all the boundary elements of the convective dynamic calculation domain whether it satisfies the conditions of converged, compressible In the flow, at the most upstream and no longer affected by the other units in the convective dynamic calculation domain are four conditions; if so, go to step S14; The unit makes a judgment; it specifically includes the following steps:

S131:对所述对流动态计算域的所有边界单元逐个判断是否满足:S131: Determine whether all boundary elements of the convective dynamic computational domain satisfy one by one:

||ΔW||/||ΔW(1)||max<εd (8)||ΔW||/||ΔW (1) || maxd (8)

其中,εd表示删除阈值,取εd≤10-7;若是,则执行步骤S132;若否,则返回步骤S131,对所述对流动态计算域的下一边界单元进行判断;Wherein, ε d represents the deletion threshold, and ε d ≤ 10 −7 ; if yes, then go to step S132; if not, return to step S131 to judge the next boundary element of the convective dynamic computational domain;

S132:判断该边界单元的马赫数是否大于0.3;若是,则执行步骤S133;若否,则返回步骤S131,对所述对流动态计算域的下一边界单元进行判断;S132: judge whether the Mach number of the boundary unit is greater than 0.3; if so, execute step S133; if not, return to step S131 to judge the next boundary unit of the convective dynamic computational domain;

S133:判断该边界单元是否满足:S133: Determine whether the boundary unit satisfies:

Figure BDA0002534371280000071
Figure BDA0002534371280000071

其中,q取对流动态计算域边界单元的格心指向边界单元紧邻单元格心的单位向量,θd表示上游单元容差角,取5°≤θd≤10°;若是,则执行步骤S134;若否,则返回步骤S131,对所述对流动态计算域的下一边界单元进行判断;Wherein, q takes the unit vector of the cell center of the boundary cell of the convective dynamic calculation domain pointing to the cell center of the boundary cell immediately adjacent to the cell center, and θ d represents the tolerance angle of the upstream cell, which is 5° ≤θd≤10 °; if so, execute step S134; If not, return to step S131 to judge the next boundary element of the convective dynamic computational domain;

S134:对于超声速无粘流动,无需再判断,直接执行步骤S14;S134: For supersonic inviscid flow, no further judgment is required, and step S14 is directly executed;

对于亚声速流动和粘性流动,判断该边界单元是否满足:For subsonic flow and viscous flow, determine whether the boundary element satisfies:

||ΔW||+||Δ(ΔW)||<εd (10)||ΔW||+||Δ(ΔW)||<ε d (10)

Figure BDA0002534371280000081
Figure BDA0002534371280000081

其中,||Δ(ΔW)||表示边界单元守恒量更新量受边界单元对流动态计算域中紧邻单元影响的模值;Δt表示迭代步长,CCFL表示时间推进格式的CFL数,ΔRi表示对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响;Among them, ||Δ(ΔW)|| represents the modulo value of the update quantity of the boundary element conservation quantity affected by the adjacent elements in the convective dynamic calculation domain of the boundary element; Δt represents the iteration step size, C CFL represents the CFL number of the time-advance format, ΔR i Represents the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction;

对于亚声速无粘单元,对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响ΔRi表示为:For subsonic inviscid elements, the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction is expressed as:

Figure BDA0002534371280000082
Figure BDA0002534371280000082

其中,ΔFc表示对流通量的变化量;下标i+1表示对流动态计算域的边界单元沿正i方向的紧邻单元,下标i-1表示对流动态计算域的边界单元沿负i方向的紧邻单元,ΔSi+1/2表示对流动态计算域的边界单元沿正i方向单元面的面积;ΔSi-1/2表示对流动态计算域的边界单元沿负i方向单元面的面积;Among them, ΔF c represents the change of convective flux; the subscript i+1 represents the adjacent unit of the boundary element of the convective dynamic computing domain along the positive i direction, and the subscript i-1 represents the boundary element of the convective dynamic computing domain along the negative i direction ΔS i+1/2 represents the area of the boundary element of the convective dynamic computational domain along the positive i direction; ΔS i-1/2 represents the area of the boundary element of the convective dynamic computational domain along the negative i direction;

对于粘性单元,对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响ΔRi表示为:For a viscous element, the influence ΔR i of the adjacent element in the convective dynamic computational domain on the residual term of the boundary element along the i direction is expressed as:

Figure BDA0002534371280000083
Figure BDA0002534371280000083

若是,则执行步骤S14;若否,则返回步骤S131,对所述对流动态计算域的下一边界单元进行判断。If yes, go to step S14; if no, go back to step S131 to judge the next boundary unit of the convective dynamic computational domain.

在一种可能的实现方式中,在本发明提供的上述多重网格扰动域更新加速方法中,步骤S16,对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导,具体包括:In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, step S16 is to judge whether all boundary elements of the viscous dynamic calculation domain are dominated by the viscous effect one by one, which specifically includes:

对所述粘性动态计算域的所有边界单元逐个判断是否满足:Judge whether all boundary elements of the viscous dynamic computational domain satisfy one by one:

Figure BDA0002534371280000084
Figure BDA0002534371280000084

Figure BDA0002534371280000091
Figure BDA0002534371280000091

其中,Φ表示粘性新增阈值的缩比因子,||ΔW||max表示对流动态计算域所有单元的守恒量更新量在当前迭代步的最大值,(||ΔW||max)min表示第1迭代步至当前迭代步中单步守恒量更新量最大值的最小值。Among them, Φ represents the scaling factor of the newly added viscosity threshold, ||ΔW|| max represents the maximum value of the update amount of the conservation quantity of all elements in the convective dynamic computational domain in the current iteration step, (||ΔW|| max ) min represents the th 1 Iteration step to the minimum value of the maximum value of the single-step conservation amount update amount in the current iteration step.

本发明提供的上述飞行器绕流数值模拟的多重网格扰动域更新加速方法,残差估计仅在各层网格的对流动态计算域内计算残差的无粘项,仅在各层网格的粘性动态计算域内计算残差的粘性项,时间积分也仅在各层网格的对流动态计算域内执行,因此,本发明通过在各层网格采用对流、粘性两类动态计算域,可同时降低多重网格方法模拟飞行器绕流问题时各层网格的计算量,显著提升多重网格方法针对飞行器绕流数值模拟的计算效率,从而有效解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类;其中,对于第一类数据,以静态数据结构存储预设计算域中所有单元的信息;对于第二类数据,以动态数据结构仅存储各网格层对流动态计算域中网格单元的信息,并在每个迭代步及时根据各网格层的对流动态计算域调整第二类数据的存储空间,因此,相比于现有多重网格方法,本发明可有效降低飞行器绕流数值模拟的存储空间需求。针对飞行器流场的细网格,根据扰动传播特性仅将受扰动单元加入对流动态计算域中;通过衡量粘性效应的影响,仅将受粘性效应主导的网格单元加入细网格的粘性动态计算域中;利用求解收敛特征,在不影响收敛速率的前提下,将已收敛单元从对流、粘性两类动态计算域中移除;根据细网格的动态计算域,更新飞行器流场的所有粗网格的对流、粘性两类动态计算域;因此,本发明可始终保持在飞行器流场的所有网格层仅对求解未收敛的受扰单元求解,仅在粘性效应主导单元中考虑粘性,从而可有效避免现有多重网格方法在飞行器绕流数值模拟中的无效计算,解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。In the multi-grid disturbance domain update acceleration method for the numerical simulation of flow around an aircraft provided by the present invention, the residual estimation only calculates the inviscid term of the residual in the In the dynamic calculation domain, the viscous term of the residual is calculated, and the time integration is only performed in the convective dynamic calculation domain of each layer grid. Therefore, the present invention adopts two types of dynamic calculation domains, convective and viscous, in each layer grid. When the grid method simulates the flow around the aircraft, the calculation amount of each layer of grids significantly improves the calculation efficiency of the multi-grid method for the numerical simulation of the flow around the aircraft, thus effectively solving the problem that the calculation efficiency of numerical simulation is difficult to meet the increasing practical application of aircraft aerodynamic design. question of needs. The data to be stored is divided into two categories: the inherent information of all units and the information related to solving and updating; wherein, for the first type of data, the information of all units in the preset computational domain is stored in a static data structure; for the second type of data data, only the information of grid cells in the convective dynamic computing domain of each grid layer is stored in a dynamic data structure, and the storage space of the second type of data is adjusted according to the convective dynamic computing domain of each grid layer in time at each iteration step. Therefore, , compared with the existing multi-grid method, the invention can effectively reduce the storage space requirement of the numerical simulation of the flow around the aircraft. For the fine grid of the aircraft flow field, only the disturbed elements are added to the convective dynamic calculation domain according to the disturbance propagation characteristics; by weighing the influence of the viscous effect, only the grid elements dominated by the viscous effect are added to the viscous dynamic calculation of the fine grid In the domain; using the solution convergence feature, without affecting the convergence rate, the converged elements are removed from the two types of dynamic calculation domains of convection and viscosity; Convective and viscous two types of dynamic computational domains of the grid; therefore, the present invention can always keep all grid layers in the flow field of the aircraft to solve only the unconverged disturbed element, and only consider the viscosity in the viscous effect dominant element, thus The ineffective calculation of the existing multi-grid method in the numerical simulation of the flow around the aircraft can be effectively avoided, and the problem that the numerical simulation calculation efficiency is difficult to meet the increasing practical application requirements of the aerodynamic design of the aircraft can be solved.

附图说明Description of drawings

图1为本发明实施例1提供的多重网格扰动域更新加速方法的流程图;1 is a flowchart of a multi-grid disturbance domain update acceleration method provided in Embodiment 1 of the present invention;

图2为本发明实施例1中多重网格根据细网格生成粗网格的示意图;2 is a schematic diagram of generating a coarse grid from a multi-grid according to a fine grid in Embodiment 1 of the present invention;

图3为本发明实施例1中多重网格循环示意图;3 is a schematic diagram of a multi-grid cycle in Embodiment 1 of the present invention;

图4为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速NACA0012翼型无粘绕流的流场与对流动态计算域的演化过程图;Fig. 4 utilizes the multi-grid disturbance domain update acceleration method provided by the embodiment of the present invention to solve the flow field of the transonic NACA0012 airfoil inviscid flow and the evolution process diagram of the convective dynamic calculation domain;

图5a为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速NACA0012翼型无粘绕流的对流动态计算域网格量变化曲线;Fig. 5a is the grid quantity change curve of the convective dynamic computational domain for solving the inviscid flow around the transonic NACA0012 airfoil by utilizing the multi-grid disturbance domain update acceleration method provided by the embodiment 1 of the present invention;

图5b为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速NACA0012翼型无粘绕流的收敛曲线;Fig. 5b is the convergence curve of solving the inviscid flow around the transonic NACA0012 airfoil using the multi-grid disturbance domain update acceleration method provided by the embodiment 1 of the present invention;

图6为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速RAE2822翼型湍流绕流的流场与对流、粘性两类动态计算域的演化过程图;6 is a diagram showing the evolution process of the flow field and two types of dynamic calculation domains of convection and viscosity for solving the flow field of the transonic RAE2822 airfoil turbulent flow around the airfoil using the multi-grid disturbance domain update acceleration method provided in Embodiment 1 of the present invention;

图7a为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速RAE2822翼型湍流绕流的对流、粘性两类动态计算域网格量变化曲线;Fig. 7a is the grid quantity change curve of two types of dynamic computational domains, convective and viscous, for solving the turbulent flow around the transonic RAE2822 airfoil by using the multi-grid disturbance domain update acceleration method provided by Embodiment 1 of the present invention;

图7b为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速RAE2822翼型湍流绕流的收敛曲线。Fig. 7b is a convergence curve for solving the turbulent flow around a transonic RAE2822 airfoil by using the multi-grid disturbance domain update acceleration method provided in Embodiment 1 of the present invention.

具体实施方式Detailed ways

下面将结合本发明实施方式中的附图,对本发明实施方式中的技术方案进行清楚、完整的描述,显然,所描述的实施方式仅仅是作为例示,并非用于限制本发明。The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are merely illustrative and not intended to limit the present invention.

采用多重网格方法模拟飞行器绕流的求解过程具有如下特征:(1)流场的扰动始于控制方程无法满足之处;(2)携带着间断信息的扰动会随时间以有限速度逐渐传向周围流场,尚未受到扰动的区域仍保持初始状态;(3)可压缩流动中,上游流场的收敛不会晚于下游流场;(4)无论是否分离,粘性效应仅在有限区域起主导作用。为充分利用上述4点特征,避免现有多重网格方法中的无效计算,本发明提供了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法,通过在所有网格层仅对求解未收敛的收敛单元求解,仅在粘性效应主导单元中考虑粘性,有效减少单个迭代步的计算量和达到收敛态所需的迭代步数,从而显著提高多重网格方法模拟飞行器绕流的计算效率。The solution process of simulating the flow around an aircraft by using the multi-grid method has the following characteristics: (1) the disturbance of the flow field starts at the point where the governing equations cannot be satisfied; (2) the disturbance carrying discontinuous information will gradually propagate to the direction at a finite speed with time. (3) In compressible flow, the convergence of the upstream flow field will not be later than that of the downstream flow field; (4) Regardless of separation or not, the viscous effect only dominates in a limited region effect. In order to make full use of the above-mentioned four characteristics and avoid invalid calculation in the existing multi-grid method, the present invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft. The non-converged convergent element is solved, and the viscosity is only considered in the viscous effect-dominant element, which effectively reduces the calculation amount of a single iteration step and the number of iteration steps required to reach a convergent state, thereby significantly improving the calculation efficiency of the multi-grid method for simulating the flow around the aircraft .

本发明提供的一种飞行器绕流数值模拟的多重网格扰动域更新加速方法,包括如下步骤:The invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft, comprising the following steps:

S1:读入数据,包括飞行器流场的细网格、预设计算域、边界条件和计算设置;S1: Read in data, including the fine mesh of the aircraft flow field, the preset computational domain, boundary conditions and calculation settings;

S2:通过稀疏细网格,生成指定层数的粗网格;S2: Generate a coarse grid with a specified number of layers through a sparse fine grid;

S3:根据来流条件或根据给定流场,对粗网格和细网格上预设计算域中的流场进行初始化;S3: Initialize the flow field in the preset computational domain on the coarse grid and fine grid according to the incoming flow conditions or according to the given flow field;

S4:根据流场初始化方式,建立对流动态计算域和粘性动态计算域;S4: According to the flow field initialization method, establish the convective dynamic calculation domain and the viscous dynamic calculation domain;

S5:将所有层网格求解所需存储的数据分为两类;第一类数据是单元的固有信息,采用静态数据结构存储预设计算域中所有网格单元的信息,包括网格坐标和流场变量;第二类数据是与求解更新相关的信息,采用动态数据结构仅存储对流动态计算域内网格单元的信息,包括守恒量更新量和当地迭代步长;S5: Divide the data that needs to be stored for solving all layers of grids into two categories; the first type of data is the inherent information of the unit, and the static data structure is used to store the information of all grid units in the preset computational domain, including grid coordinates and Flow field variables; the second type of data is the information related to the solution update. The dynamic data structure is used to store only the information of the grid cells in the convective dynamic calculation domain, including the update amount of the conservation quantity and the local iteration step size;

S6:根据边界条件类型,为当前层网格的边界虚网格的守恒量赋值;S6: According to the boundary condition type, assign value to the conservation quantity of the boundary virtual grid of the current layer grid;

S7:将流动控制方程的残差项分为无粘项与粘性项两类;其中,无粘项包括对流通量和粗网格的强迫函数,无粘项在对流动态计算域中计算;粘性项包括粘性通量和湍流模型方程源项,粘性项在粘性动态计算域中计算;S7: The residual term of the flow control equation is divided into two types: inviscid term and viscous term; among them, the inviscid term includes the forcing function of convective flux and coarse grid, and the inviscid term is calculated in the convective dynamic calculation domain; the viscous term terms include viscous flux and turbulence model equation source terms, and the viscous term is calculated in the viscous dynamic computational domain;

S8:在对流动态计算域中,求解守恒量更新量,更新流场变量;S8: In the convective dynamic calculation domain, solve the update quantity of the conserved quantity and update the flow field variables;

S9:判断当前网格层是否为细网格;若是,则执行步骤S10;若否,则执行步骤S21;S9: determine whether the current grid layer is a fine grid; if yes, go to step S10; if not, go to step S21;

S10:判断细网格是否达到收敛条件;若否,则执行步骤S11;若是,则执行步骤S23;S10: Determine whether the fine mesh reaches the convergence condition; if not, execute step S11; if yes, execute step S23;

S11:对对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动;若是,则执行步骤S12;若否,则返回步骤S11,对对流动态计算域的下一边界单元进行判断;遍历对流动态计算域的所有边界单元后,执行步骤S13;S11: Judge whether all boundary elements of the convective dynamic calculation domain have been subjected to inviscidity disturbance one by one; if so, go to step S12; if not, return to step S11 to judge the next boundary element of the convective dynamic calculation domain; After all boundary elements of the domain are dynamically calculated, step S13 is performed;

S12:衡量无粘扰动的传播方向,并将该边界单元的紧邻单元中位于传播方向上的单元纳入对流动态计算域,返回步骤S11,对对流动态计算域的下一边界单元进行判断;S12: Measure the propagation direction of the inviscid disturbance, and include the unit located in the propagation direction in the adjacent units of the boundary unit into the convective dynamic calculation domain, and return to step S11 to judge the next boundary unit of the convective dynamic calculation domain;

S13:对对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受对流动态计算域中其他单元的影响四个条件;若是,则执行步骤S14;若否,则返回步骤S13,对对流动态计算域的下一边界单元进行判断;遍历对流动态计算域的所有边界单元后,执行步骤S16;S13: Judging, one by one, whether all the boundary elements in the convective dynamic computational domain satisfy the four conditions: converged, located in the compressible flow, located in the most upstream, and no longer affected by other elements in the convective dynamic computational domain; if so, go to step S13 S14; if not, return to step S13 to judge the next boundary element of the convective dynamic computational domain; after traversing all boundary elements of the convective dynamic computational domain, execute step S16;

S14:将该边界单元从对流动态计算域中移除,并判断该边界单元是否存在于粘性动态计算域中;若是,则执行步骤S15;若否,则返回步骤S13,对对流动态计算域的下一边界单元进行判断;S14: remove the boundary element from the convective dynamic calculation domain, and determine whether the boundary element exists in the viscous dynamic calculation domain; if so, go to step S15; The next boundary unit is judged;

S15:将该边界单元从粘性动态计算域中移除,返回步骤S13,对对流动态计算域的下一边界单元进行判断;S15: remove the boundary element from the viscous dynamic calculation domain, and return to step S13 to judge the next boundary element of the convective dynamic calculation domain;

S16:对粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导;若是,则执行步骤S17;若否,则返回步骤S16,对粘性动态计算域的下一边界单元进行判断;遍历粘性动态计算域的所有边界单元后,执行步骤S18;S16: Determine whether all the boundary elements of the viscous dynamic calculation domain are dominated by the viscous effect one by one; if so, go to step S17; if not, return to step S16 to judge the next boundary element of the viscous dynamic calculation domain; After calculating all boundary elements of the domain, perform step S18;

S17:将该边界单元的紧邻单元中位于对流动态计算域中的单元纳入粘性动态计算域,返回步骤S16,对粘性动态计算域的下一边界单元进行判断;S17: Include the unit located in the convective dynamic calculation domain in the adjacent units of the boundary unit into the viscous dynamic calculation domain, and return to step S16 to judge the next boundary element of the viscous dynamic calculation domain;

S18:根据细网格的对流动态计算域和粘性动态计算域,更新所有层粗网格的对流动态计算域和粘性动态计算域;S18: According to the convective dynamic calculation domain and the viscous dynamic calculation domain of the fine grid, update the convective dynamic calculation domain and the viscous dynamic calculation domain of the coarse grid of all layers;

S19:根据对流动态计算域,重新分配所有粗网格和细网格中第二类数据的存储空间;S19: Reallocate the storage space of the second type of data in all coarse grids and fine grids according to the convective dynamic computational domain;

S20:将守恒量插值至更粗网格,返回步骤S6,进入更粗网格的计算;S20: Interpolate the conserved amount to a coarser grid, return to step S6, and enter the calculation of the coarser grid;

S21:判断是否需要计算更粗网格;若是,则返回步骤S20;若否,则执行步骤S22;S21: Determine whether a coarser grid needs to be calculated; if yes, go back to step S20; if not, go to step S22;

S22:将守恒量修正量插值回更细网格,并判断更细网格是否为细网格;若是,则返回步骤S6,进入下一迭代步的计算;若否,则返回步骤S6,进入更细网格的计算;S22: Interpolate the correction amount of the conservation amount back to a finer grid, and determine whether the finer grid is a fine grid; if so, return to step S6, and enter the calculation of the next iteration step; if not, return to step S6, enter Calculation of finer grids;

S23:输出结果。S23: output the result.

下面通过一个具体的实施例对本发明提供的上述多重网格扰动域更新加速方法的具体实施进行详细说明。The specific implementation of the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention will be described in detail below through a specific embodiment.

实施例1:如图1所示,本发明提供的上述多重网格扰动域更新加速方法具体包括如下步骤:Embodiment 1: As shown in FIG. 1 , the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention specifically includes the following steps:

第一步:数据读入。The first step: data reading.

读入数据包括飞行器流场的细网格、预设计算域、边界条件和计算设置。The read-in data includes a fine mesh of the aircraft flow field, a preset computational domain, boundary conditions, and computational settings.

第二步:生成粗网格。Step 2: Generate a coarse mesh.

通过稀疏细网格,生成指定层数的粗网格。对于结构网格,常用的粗网格生成方法是将细网格沿所有网格方向稀疏一倍,图2以一维为例示意了该粗网格生成方法。图2中,h代表细网格,2h、4h分别代表第1层、第2层粗网格。即,第1层粗网格的每个单元包含2个细网格单元,第2层粗网格中的每个单元包含2个第1层粗网格单元。Generates a coarse mesh with a specified number of layers by sparse fine meshing. For structured grids, the commonly used method for generating coarse grids is to sparse the fine grids twice along all grid directions. Figure 2 illustrates the method for generating coarse grids by taking one dimension as an example. In Figure 2, h represents the fine grid, 2h and 4h represent the first layer and the second layer coarse grid, respectively. That is, each cell of the first layer coarse grid contains 2 fine grid cells, and each cell in the second layer coarse grid contains 2 first layer coarse grid cells.

第三步:流场初始化。The third step: flow field initialization.

流场初始化分为根据来流条件初始化和根据给定流场初始化两种方式。当根据来流条件初始化时,将粗网格和细网格上预设计算域中所有网格单元的守恒量赋为来流值。当根据给定流场初始化时,将细网格上预设计算域中所有网格单元的守恒量赋为给定流场值,粗网格上预设计算域中所有网格单元的守恒量根据细网格上所有网格单元的守恒量插值获得。以一维为例,粗网格单元的守恒量满足:The flow field initialization is divided into two ways: initialization according to incoming flow conditions and initialization according to a given flow field. When initialized according to inflow conditions, the conserved quantities of all grid cells in the preset computational domain on coarse and fine grids are assigned as inflow values. When initialized according to the given flow field, the conservation quantities of all grid cells in the preset computational domain on the fine grid are assigned to the given flow field value, and the conservation quantities of all grid cells in the pre-set computational domain on the coarse grid Interpolated from the conserved quantities of all grid cells on a fine grid. Taking one dimension as an example, the conservation of coarse grid cells satisfies:

Figure BDA0002534371280000141
Figure BDA0002534371280000141

其中,W2h表示粗网格网格单元的守恒量,Wh表示细网格网格单元的守恒量;|Ω|2h表示粗网格网格单元的体积,|Ω|h表示细网格网格单元的体积;2I-1表示沿网格方向第2I-1个单元,2I表示沿网格方向第2I个单元。Among them, W 2h represents the conserved quantity of the coarse mesh cell, W h represents the conservation quantity of the fine mesh cell; |Ω| 2h represents the volume of the coarse mesh cell, and |Ω| h represents the fine mesh The volume of grid cells; 2I-1 represents the 2I-1th cell along the grid direction, and 2I represents the 2Ith cell along the grid direction.

第四步:建立对流动态计算域和粘性动态计算域。The fourth step: establish the convective dynamic computational domain and the viscous dynamic computational domain.

(1)当根据来流条件初始化时,对流动态计算域和粘性动态计算域均根据壁面边界建立。携带着流场变化信息的扰动从物体表面产生,受粘性主导的单元也紧邻物体表面,因此,对流动态计算域取壁面的若干层相邻单元作为初始单元。为避免对流动态计算域的增长速度对多重网格法的收敛速率产生影响,最粗层网格取紧邻壁面的1层单元作为初始单元,则对于结构网格框架下的细网格,其对流动态计算域需取壁面的L层相邻单元作为初始单元,L满足:(1) When initialized according to the inflow conditions, both the convective dynamic calculation domain and the viscous dynamic calculation domain are established according to the wall boundary. The disturbance carrying flow field change information is generated from the surface of the object, and the units dominated by viscosity are also close to the surface of the object. Therefore, the convective dynamic computational domain takes several layers of adjacent elements on the wall as the initial units. In order to avoid the influence of the growth rate of the convective dynamic computational domain on the convergence rate of the multi-grid method, the first layer of elements next to the wall is taken as the initial element for the coarsest grid, and for the fine grid under the structural grid framework, the convective The dynamic calculation domain needs to take the L-layer adjacent element of the wall as the initial element, and L satisfies:

L=2l (2)L=2 l (2)

其中,l表示粗网格层数,l为正整数,由步骤S1读入的计算设置给定;较粗层网格的初始单元根据细网格的初始单元生成;对于粘性动态计算域,所有网格层均取紧邻壁面的1层单元作为初始单元,以避免粘性动态计算域增长过快而导致无效计算。Among them, l represents the number of coarse grid layers, and l is a positive integer, which is given by the calculation settings read in step S1; the initial elements of the coarser grids are generated according to the initial elements of the fine grid; for the viscous dynamic computational domain, all The grid layer is taken as the initial element of the 1-layer element next to the wall surface, so as to avoid the invalid calculation caused by the rapid growth of the viscous dynamic computational domain.

(2)当根据给定流场初始化时,对流动态计算域和粘性动态计算域根据给定流场的流动特性建立。对于对流动态计算域,细网格的初始单元为给定流场中的已受扰单元,这些单元的流动特性与来流条件不符,满足:(2) When initialized according to a given flow field, the convective dynamic computational domain and the viscous dynamic computational domain are established according to the flow characteristics of the given flow field. For the convective dynamic computational domain, the initial elements of the fine grid are the disturbed elements in the given flow field, and the flow characteristics of these elements do not match the incoming flow conditions and satisfy:

||W-W||/||ΔW(1)||max>εa,c (3)||WW ||/||ΔW (1) || max >ε a,c (3)

其中,W表示来流条件的守恒量;||ΔW(1)||max表示对流动态计算域所有单Among them, W represents the conserved quantity of the incoming flow condition; ||ΔW (1) || max represents all units in the convective dynamic computational domain

元的守恒量更新量在第1迭代步的最大值;εa,c表示对流新增阈值,取10-4≤εa,c≤10-6The maximum value of the update quantity of the conservation quantity of the element in the first iteration step; ε a,c represents the new threshold of convection, which is 10 -4 ≤ε a,c ≤10 -6 .

对于粘性动态计算域,细网格的初始单元为细网格对流动态计算域中的粘性效应主导单元。定义粘性效应衡量参数Ψ,令其表示粘性扰动与无粘扰动的质量流量之比。粘性效应主导单元的粘性效应衡量参数应处于较大量级,故粘性动态计算域中的单元应满足:For the viscous dynamic computational domain, the initial element of the fine grid is the viscous effect dominant element in the fine grid convective dynamic computational domain. The viscous effect measure parameter Ψ is defined so that it represents the ratio of the mass flow rate between viscous disturbances and inviscid disturbances. The viscous effect measurement parameter of the viscous effect-dominant element should be in a relatively large order, so the elements in the viscous dynamic computational domain should satisfy:

Figure BDA0002534371280000151
Figure BDA0002534371280000151

Figure BDA0002534371280000152
Figure BDA0002534371280000152

其中,Ψ表示粘性效应衡量参数,为粘性扰动与无粘扰动的质量流量之比;

Figure BDA0002534371280000153
表示第1迭代步的粘性效应衡量参数;I,J,K分别表示网格方向;
Figure BDA0002534371280000154
表示对流通量Jacobian矩阵沿i方向的谱半径,
Figure BDA0002534371280000155
表示粘性通量Jacobian矩阵沿i方向的谱半径;εa,v表示粘性新增阈值,取10-3≤εa,v≤10-2。Among them, Ψ represents the viscous effect measurement parameter, which is the ratio of viscous disturbance to mass flow without viscous disturbance;
Figure BDA0002534371280000153
Represents the viscous effect measurement parameter of the first iteration step; I, J, K represent the grid direction respectively;
Figure BDA0002534371280000154
represents the spectral radius of the Jacobian matrix of convective flux along the i direction,
Figure BDA0002534371280000155
Represents the spectral radius of the viscous flux Jacobian matrix along the i direction; ε a,v represents the new threshold of viscosity, which is 10 -3 ≤ε a,v ≤10 -2 .

所有粗网格层的对流动态计算域依据细网格的对流动态计算域生成,所有粗网格层的粘性动态计算域依据细网格的粘性动态计算域生成。The convective dynamic computational domains of all coarse grid layers are generated according to the convective dynamic computational domains of fine grids, and the viscous dynamic computational domains of all coarse grid layers are generated from the viscous dynamic computational domains of fine grids.

第五步:分配存储空间。Step 5: Allocate storage space.

将所有层网格求解所需存储的数据分为两类;第一类数据是单元的固有信息,采用静态数据结构存储预设计算域中所有网格单元的信息,包括网格坐标和流场变量;第二类数据是与求解更新相关的信息,采用动态数据结构仅存储对流动态计算域内网格单元的信息,包括守恒量更新量和当地迭代步长。The data needed to be stored in all layers of grid solution are divided into two categories; the first type of data is the inherent information of the unit, and the static data structure is used to store the information of all grid units in the preset computational domain, including grid coordinates and flow field. Variables; the second type of data is the information related to the solution update. The dynamic data structure is used to store only the information of the grid cells in the convective dynamic calculation domain, including the update amount of the conserved quantity and the local iterative step size.

对于第一类数据,采用静态数据结构存储预设计算域中所有网格单元的信息。例如,令I,J,K表示网格方向,Imax,Jmax,Kmax分别表示各网格方向的单元数,对于一块Imax×Jmax×Kmax的网格,若每个单元有10个量需要记录,则可采用一个实数型四维数组A1(10,Imax,Jmax,Kmax)进行存储。For the first type of data, a static data structure is used to store the information of all grid cells in the preset computational domain. For example, let I, J, K denote the grid direction, and I max , J max , and K max respectively denote the number of elements in each grid direction. For a grid of I max ×J max ×K max , if each element has If 10 quantities need to be recorded, a real four-dimensional array A 1 (10, I max , J max , K max ) can be used for storage.

对于第二类数据,采用动态数据结构仅存储对流动态计算域内网格单元的信息。例如,令I,J,K表示网格方向,具有相同J,K标号的单元统称为一行,定义一个特殊数据类型用于存储对流、粘性两类动态计算域中一行的信息,该特殊数据类型包括:a.一条单项链表——在两类动态计算域更新中记录当前迭代步该行单元的I标号范围;b.一个整型二维数组——在两类动态计算域更新时存储上一迭代步该行单元的I标号范围,在两类动态计算域更新后存储当前迭代步该行单元的I标号范围;c.一个实数型二维数组——存储该行位于对流动态计算域中单元的第二类数据,与整型二维数组的I标号范围相对应。对于一块Imax×Jmax×Kmax的网格,若每个单元有10个量需要记录,则采用一个特殊数据类型二维数组A2(Jmax,Kmax),A2中实数型二维数组为B(10,I1:I2),I1、I2分别表示行(J,K)中单元的I标号的最小值与最大值,1<I1,I2<ImaxFor the second type of data, the dynamic data structure is used to store only the information of grid cells in the convective dynamic computational domain. For example, let I, J, K denote the grid direction, the cells with the same J, K labels are collectively called a row, and a special data type is defined to store the information of a row in the two types of dynamic calculation domains of convection and viscosity. This special data type Including: a. a single-item list - record the I label range of the row unit in the current iteration step in the update of the two types of dynamic calculation fields; b. an integer two-dimensional array - store the last value in the update of the two types of dynamic calculation fields The range of the I label of the row unit in the iteration step, and the range of the I label of the row unit of the current iteration step is stored after the two types of dynamic calculation domains are updated; c. A real number two-dimensional array - to store the unit located in the convective dynamic calculation domain of the row The second type of data corresponds to the range of the I label of the integer two-dimensional array. For a grid of I max ×J max ×K max , if each unit has 10 quantities to be recorded, a special data type two-dimensional array A 2 (J max ,K max ) is used, and the real number type two in A 2 is used. The dimensional array is B(10, I 1 :I 2 ), I 1 and I 2 represent the minimum and maximum values of the I labels of the cells in the row (J, K), respectively, 1<I 1 , I 2 <I max .

现有多重网格方法中,所有网格层中所有网格单元的所有数据均存储在静态数据结构中。相比之下,本发明提供的多重网格扰动域更新加速方法,将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类;其中,对于第一类数据,以静态数据结构存储预设计算域中所有网格单元的信息;对于第二类数据,以动态数据结构仅存储各网格层对流动态计算域中网格单元的信息,并且,第十五步会在每个迭代步及时根据各网格层的对流动态计算域调整第二类数据的存储空间。因此,相比于现有多重网格方法,本发明可有效降低存储空间需求。In existing multigrid methods, all data for all grid cells in all grid layers are stored in static data structures. In contrast, the multi-grid disturbance domain update acceleration method provided by the present invention divides the data to be stored into two categories: inherent information of all units and information related to solution update; wherein, for the first type of data, the The static data structure stores the information of all grid cells in the preset computing domain; for the second type of data, the dynamic data structure only stores the information of the grid cells in the convective dynamic computing domain of each grid layer, and the fifteenth step will In each iteration step, the storage space of the second type of data is adjusted according to the convective dynamic computational domain of each grid layer in time. Therefore, compared with the existing multi-grid method, the present invention can effectively reduce the storage space requirement.

第六步:边界条件处理。Step 6: Boundary condition processing.

根据边界条件类型,为当前层网格的边界虚网格的守恒量赋值。According to the boundary condition type, assign the value of the conservation quantity of the boundary virtual grid of the current layer grid.

第七步:残差估计。Step 7: Residual estimation.

流动控制方程表示为:The flow control equation is expressed as:

Figure BDA0002534371280000161
Figure BDA0002534371280000161

其中,t表示时间,Fc表示对流通量,Fv表示粘性通量,QF表示粗网格的强迫函数,QT表示湍流模型方程源项,Nf表示单元面的面数,ΔS表示单元面的面积;式(6)的右端项统称为流动控制方程的残差项,

Figure BDA0002534371280000162
为无粘项,
Figure BDA0002534371280000171
为粘性项;在对流动态计算域中求解残差项的无粘项,在粘性动态计算域中求解残差项的粘性项。where t is the time, F c is the convective flux, F v is the viscous flux, Q F is the forcing function of the coarse grid, Q T is the source term of the turbulence model equation, N f is the number of element faces, and ΔS is the is the area of the element face; the right-hand term of equation (6) is collectively referred to as the residual term of the flow control equation,
Figure BDA0002534371280000162
for non-sticky terms,
Figure BDA0002534371280000171
is the viscous term; the inviscid term of the residual term is solved in the convective dynamic computational domain, and the viscous term of the residual term is solved in the viscous dynamic computational domain.

第八步:时间积分。Step 8: Time integration.

在对流动态计算域中,采用时间推进格式计算式(6)的左端项;所采用的时间积分方法与现有方法一致。In the convective dynamic computational domain, the left-hand term of Equation (6) is calculated using the time advance scheme; the time integration method used is consistent with the existing method.

现有多重网格方法中,最为耗时的残差估计、时间积分两步需对静态的预设计算域内各层网格的所有网格单元执行,约占99%的总求解时间。相比之下,本发明提供的多重网格扰动域更新加速方法,第七步的残差估计仅在各层网格的对流动态计算域内计算残差的无粘项,仅在各层网格的粘性动态计算域内计算残差的粘性项;第八步的时间积分也仅在各层网格的对流动态计算域内执行。因此,本发明可同时降低各层网格的计算量,显著提升多重网格方法的计算效率。Among the existing multi-grid methods, the two most time-consuming steps of residual estimation and time integration need to be performed on all grid elements of each grid in the static preset computational domain, accounting for about 99% of the total solution time. In contrast, in the multi-grid disturbance domain update acceleration method provided by the present invention, the residual estimation in the seventh step only calculates the inviscid term of the residual in the convective dynamic calculation domain of each layer of grids, and only calculates the inviscid term of the residual in each layer of grids. The viscous term of the residual is calculated in the viscous dynamic calculation domain of ; the time integration of the eighth step is also only performed in the convective dynamic calculation domain of each layer of grids. Therefore, the present invention can simultaneously reduce the calculation amount of each layer of grid, and significantly improve the calculation efficiency of the multi-grid method.

第九步:判断当前网格层是否为细网格;若是,则执行第十步;若否,则执行第十七步。Step 9: Determine whether the current grid layer is a fine grid; if so, perform Step 10; if not, perform Step 17.

第十步:判断细网格是否达到收敛条件;若否,则执行第十一步;若是,则执行第十九步。Step 10: judge whether the fine mesh meets the convergence condition; if not, execute step 11; if yes, execute step 19.

第十一步:增大对流动态计算域。Step 11: Increase the convective dynamic computational domain.

对对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动,即对对流动态计算域的所有边界单元逐个判断是否满足:Judging whether all the boundary elements of the convective dynamic calculation domain have been subjected to inviscid perturbation one by one, that is, to judge whether all the boundary elements of the convective dynamic calculation domain satisfy:

||ΔW||/||ΔW(1)||max>εa,c (7)||ΔW||/||ΔW (1) || max >ε a,c (7)

其中,||ΔW||表示对流动态计算域的边界单元在当前迭代步中守恒量更新量的模值;若对流动态计算域的边界单元满足式(7),则衡量无粘扰动的传播方向,并将该边界单元的紧邻单元中位于传播方向上的单元纳入对流动态计算域,返回第十一步,对对流动态计算域的下一边界单元进行判断;若对流动态计算域的边界单元不满足式(7),则返回第十一步,对对流动态计算域的下一边界单元进行判断;遍历对流动态计算域的所有边界单元后,执行第十二步。Among them, ||ΔW|| represents the modulo value of the update amount of the conservation quantity of the boundary element of the convective dynamic calculation domain in the current iteration step; if the boundary element of the convective dynamic calculation domain satisfies Equation (7), the propagation direction of the inviscid disturbance is measured. , and the unit located in the propagation direction in the adjacent units of the boundary unit is included in the convective dynamic calculation domain, and returns to the eleventh step to judge the next boundary unit of the convective dynamic calculation domain; if the boundary element of the convective dynamic calculation domain is not If formula (7) is satisfied, return to step 11 to judge the next boundary element of the convective dynamic computational domain; after traversing all boundary elements of the convective dynamic computational domain, perform the twelfth step.

无粘扰动相对于流动以声速传播,沿无粘扰动的传播方向,无粘扰动的传播速度为正;令q表示单位方向向量,则无粘扰动沿q方向传播表示为:The inviscid disturbance propagates at the speed of sound relative to the flow, and along the propagation direction of the inviscid disturbance, the propagation velocity of the inviscid disturbance is positive; let q represent the unit direction vector, then the inviscid disturbance propagates along the q direction as:

u·q+a>0 (8)u·q+a>0 (8)

其中,u表示流动速度矢量,a表示声速;q取对流动态计算域边界单元的格心指向边界单元某一格点的单位向量,若q满足式(8),则将与对流动态计算域边界单元共享该格点的紧邻单元纳入对流动态计算域。Among them, u is the flow velocity vector, a is the speed of sound; q is the unit vector of the grid center of the boundary element of the convective dynamic calculation domain pointing to a certain grid point of the boundary element. Cells that share this grid are included in the convective dynamic computational domain.

第十二步:缩小对流动态计算域和粘性动态计算域。The twelfth step: reduce the convective dynamic computational domain and the viscous dynamic computational domain.

为了在不影响收敛速率的前提下,尽可能地减少单个迭代步的计算量,可以缩小对流动态计算域的边界单元。对对流动态计算域的所有边界单元逐个判断以下条件实现对流动态计算域的缩小:In order to reduce the calculation amount of a single iteration step as much as possible without affecting the convergence rate, the boundary elements of the convective dynamic computational domain can be reduced. The following conditions are judged for all boundary elements of the convective dynamic computational domain one by one to reduce the convective dynamic computational domain:

(1)对对流动态计算域的所有边界单元逐个判断是否已收敛,即对对流动态计算域的所有边界单元逐个判断是否满足:(1) Judge whether all the boundary elements of the convective dynamic calculation domain have converged one by one, that is, judge whether all the boundary elements of the convective dynamic calculation domain satisfy:

||ΔW||/||ΔW(1)||max<εd (9)||ΔW||/||ΔW (1) || maxd (9)

其中,εd表示删除阈值,取εd≤10-7;若该边界单元满足式(9),则继续判断(2);若否,则对对流动态计算域的下一边界单元进行判断;Wherein, ε d represents the deletion threshold, and takes ε d ≤10 -7 ; if the boundary element satisfies the formula (9), continue to judge (2); if not, judge the next boundary element of the convective dynamic computational domain;

(2)判断该边界单元是否位于可压缩流动中,即判断该边界单元的马赫数是否大于0.3;若是,则继续判断(3);若否,则返回判断(1),对对流动态计算域的下一边界单元进行判断;(2) Determine whether the boundary element is in a compressible flow, that is, determine whether the Mach number of the boundary element is greater than 0.3; if so, continue to judge (3); The next boundary unit is judged;

(3)判断该边界单元是否位于最上游,即判断该边界单元是否满足:(3) Judging whether the boundary element is located at the most upstream, that is, judging whether the boundary element satisfies:

Figure BDA0002534371280000181
Figure BDA0002534371280000181

其中,q取对流动态计算域边界单元的格心指向边界单元紧邻单元格心的单位向量,θd表示上游单元容差角,取5°≤θd≤10°;若该边界单元满足式(10),则继续判断(4);若否,则返回判断(1),对对流动态计算域的下一边界单元进行判断;Among them, q is taken as the unit vector of the cell center of the boundary element of the convective dynamic calculation domain pointing to the cell center of the adjacent cell, and θd represents the tolerance angle of the upstream element, which is 5° ≤θd≤10 °; if the boundary element satisfies the formula ( 10), then continue to judge (4); if not, return to judge (1) to judge the next boundary element of the convective dynamic computational domain;

(4)判断该边界单元是否不再受对流动态计算域中其他单元的影响,对于超声速无粘流动,则无需判断、自动满足;对于亚声速流动和粘性流动,判断该边界单元是否满足:(4) Determine whether the boundary element is no longer affected by other elements in the convective dynamic calculation domain. For supersonic inviscid flow, no judgment is required and it is automatically satisfied; for subsonic flow and viscous flow, determine whether the boundary element satisfies:

||ΔW||+||Δ(ΔW)||<εd (11)||ΔW||+||Δ(ΔW)||<ε d (11)

Figure BDA0002534371280000191
Figure BDA0002534371280000191

其中,||Δ(ΔW)||表示边界单元守恒量更新量受边界单元对流动态计算域中紧邻单元影响的模值;Δt表示迭代步长,CCFL表示时间推进格式的CFL数,ΔRi表示对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响;对于亚声速无粘单元,对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响ΔRi表示为:Among them, ||Δ(ΔW)|| represents the modulo value of the update quantity of the boundary element conservation quantity affected by the adjacent elements in the convective dynamic calculation domain of the boundary element; Δt represents the iteration step size, C CFL represents the CFL number of the time-advance format, ΔR i Represents the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction; for subsonic inviscid elements, the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction ΔR i represents for:

Figure BDA0002534371280000192
Figure BDA0002534371280000192

其中,ΔFc表示对流通量的变化量;下标i+1表示对流动态计算域的边界单元沿正i方向的紧邻单元,下标i-1表示对流动态计算域的边界单元沿负i方向的紧邻单元,ΔSi+1/2表示对流动态计算域的边界单元沿正i方向单元面的面积;ΔSi-1/2表示对流动态计算域的边界单元沿负i方向单元面的面积;对于粘性单元,对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响ΔRi表示为:Among them, ΔF c represents the change of convective flux; the subscript i+1 represents the adjacent unit of the boundary element of the convective dynamic computing domain along the positive i direction, and the subscript i-1 represents the boundary element of the convective dynamic computing domain along the negative i direction ΔS i+1/2 represents the area of the boundary element of the convective dynamic computational domain along the positive i direction; ΔS i-1/2 represents the area of the boundary element of the convective dynamic computational domain along the negative i direction; For a viscous element, the influence ΔR i of the adjacent element in the convective dynamic computational domain on the residual term of the boundary element along the i direction is expressed as:

Figure BDA0002534371280000193
Figure BDA0002534371280000193

若该边界单元满足式(11)和式(12),则该边界单元已同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件,则将该边界单元从对流动态计算域中移除,若该边界单元同时存在于粘性动态计算域中,则将该边界单元从粘性动态计算域中移除,返回第十二步,对对流动态计算域的下一边界单元进行判断;若该边界单元不满足式(11)和式(12),则返回第十二步,对对流动态计算域的下一边界单元进行判断。遍历对流动态计算域的所有边界单元后,执行第十三步。If the boundary element satisfies equations (11) and (12), then the boundary element has been converged, located in the compressible flow, located in the most upstream and no longer affected by other elements in the convective dynamic computational domain. If the conditions are met, the boundary element is removed from the convective dynamic calculation domain. If the boundary element exists in the viscous dynamic calculation domain at the same time, the boundary element is removed from the viscous dynamic calculation domain and returns to step 12. Judge the next boundary element of the convective dynamic computational domain; if the boundary element does not satisfy equations (11) and (12), return to the twelfth step to judge the next boundary element of the convective dynamic computational domain. After traversing all boundary elements of the convective dynamic computational domain, perform the thirteenth step.

第十三步:增大粘性动态计算域。The thirteenth step: increase the viscous dynamic computational domain.

对粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导,即对粘性动态计算域的所有边界单元逐个判断是否满足:Determine whether all boundary elements of the viscous dynamic computational domain are dominated by the viscous effect one by one, that is, to judge whether all the boundary elements of the viscous dynamic computational domain satisfy:

Figure BDA0002534371280000201
Figure BDA0002534371280000201

Figure BDA0002534371280000202
Figure BDA0002534371280000202

其中,Φ表示粘性新增阈值的缩比因子,||ΔW||max表示对流动态计算域所有单元的守恒量更新量在当前迭代步的最大值,(||ΔW||max)min表示第1迭代步至当前迭代步中单步守恒量更新量最大值的最小值;若粘性动态计算域的边界单元满足式(15)和式(16),则将该边界单元的紧邻单元中位于对流动态计算域中的单元纳入粘性动态计算域,返回第十三步,对粘性动态计算域的下一边界单元进行判断;若粘性动态计算域的边界单元不满足式(15)和式(16),则返回第十三步,对粘性动态计算域的下一边界单元进行判断;遍历粘性动态计算域的所有边界单元后,执行第十四步。流动控制方程的粘性项呈椭圆型,因此,可以在不增大对流动态计算域的前提下,将已受粘性扰动单元的所有紧邻单元均加入粘性动态计算域。Among them, Φ represents the scaling factor of the newly added viscosity threshold, ||ΔW|| max represents the maximum value of the update amount of the conservation quantity of all elements in the convective dynamic computational domain in the current iteration step, (||ΔW|| max ) min represents the th 1 Iterative step to the minimum value of the maximum value of the update value of the single-step conservation quantity in the current iteration step; if the boundary element of the viscous dynamic computational domain satisfies the equations (15) and (16), the adjacent elements of the boundary element are located in the convection. The elements in the dynamic calculation domain are included in the viscous dynamic calculation domain, and return to step 13 to judge the next boundary element of the viscous dynamic calculation domain; if the boundary elements of the viscous dynamic calculation domain do not satisfy equations (15) and (16) , then return to the thirteenth step to judge the next boundary element of the viscous dynamic computational domain; after traversing all the boundary elements of the viscous dynamic computational domain, execute the fourteenth step. The viscous term of the flow governing equation is elliptical. Therefore, all the adjacent elements of the viscous perturbed element can be added to the viscous dynamic computational domain without enlarging the convective dynamic computational domain.

第十四步:更新粗网格动态计算域。Step 14: Update the coarse grid dynamic computational domain.

根据细网格的对流动态计算域更新所有层粗网格的对流动态计算域,根据细网格的粘性动态计算域更新所有层粗网格的粘性动态计算域。对于结构网格,粗网格与细网格间动态计算域的对应关系如图2所示,若某一细网格单元包含于细网格的动态计算域中,则对应于该单元的粗网格单元也会被加入粗网格的动态计算域中。According to the convective dynamic computational domain of the fine grid, the convective dynamic computational domain of the coarse grid of all layers is updated, and the viscous dynamic computational domain of the coarse grid of all layers is updated according to the viscous dynamic computational domain of the fine grid. For the structured grid, the corresponding relationship between the coarse grid and the fine grid of the dynamic computational domain is shown in Figure 2. If a fine grid cell is included in the dynamic computational domain of the fine grid, it corresponds to the coarse grid of the cell. Grid cells are also added to the dynamic computational domain of the coarse grid.

本发明在第十一步中,针对细网格,根据扰动传播特性仅将受扰动单元加入对流动态计算域中;在第十三步中,通过衡量粘性效应的影响,仅将受粘性效应主导的网格单元加入细网格的粘性动态计算域中;在第十二步中,利用求解收敛特征,在不影响收敛速率的前提下,及时将已收敛单元从对流、粘性两类动态计算域中移除;在第十四步中,根据细网格的动态计算域,更新所有粗网格的对流、粘性两类动态计算域。因此,本发明可始终保持在所有网格层仅对求解未收敛的受扰单元求解,仅在粘性效应主导单元中考虑粘性,从而可有效避免现有多重网格方法求解中的无效计算。In the eleventh step of the present invention, for fine grids, only the disturbed units are added to the convective dynamic calculation domain according to the disturbance propagation characteristics; in the thirteenth step, only the viscous effect dominates The grid elements are added to the viscous dynamic computational domain of the fine grid; in the twelfth step, using the solution convergence feature, without affecting the convergence rate, the converged elements are transferred from the convection and viscous dynamic computational domains in time. In the fourteenth step, according to the dynamic computational domain of the fine grid, the two types of dynamic computational domains of convection and viscosity of all coarse grids are updated. Therefore, the present invention can always keep only the unconverged disturbed element to be solved in all grid layers, and only consider the viscosity in the viscous effect dominant element, thereby effectively avoiding the invalid calculation in the existing multi-grid method.

第十五步:重新分配第二类数据的存储空间。Step 15: Reallocate the storage space of the second type of data.

根据对流动态计算域,重新分配所有粗网格和细网格中第二类数据的存储空间。According to the convective dynamic computational domain, reallocate the storage space for the second type of data in all coarse and fine grids.

第十六步:将守恒量插值至更粗网格,返回第六步,进入更粗网格的计算。Step 16: Interpolate the conserved quantity to a coarser grid, return to step 6, and enter the calculation of a coarser grid.

守恒量从细网格至粗网格的插值方式与第三步中根据给定流场初始化的插值方式相同。The interpolation of the conserved quantities from the fine grid to the coarse grid is the same as the interpolation in the third step for the given flow field initialization.

第十七步:判断是否需要计算更粗网格;若需要计算更粗网格,则返回第十六步;若不需要计算更粗网格,则执行第十八步;Step 17: Determine whether a coarser grid needs to be calculated; if a coarser grid needs to be calculated, go back to step 16; if there is no need to calculate a coarser grid, execute step 18;

多重网格常用的循环方式如图3所示,左侧子图代表“V”型循环,右侧子图代表“W”型循环,图3中h代表细网格,2h代表第1层粗网格,4h代表第2层粗网格。以“V”型循环为例,若当前计算层为细网格,则第十七步便选择需要计算更粗网格;若当前计算层为最粗网格,则第十七步便选择无需计算更粗网格。The commonly used circulation mode of multi-grid is shown in Figure 3. The left subgraph represents a "V" type cycle, and the right subgraph represents a "W" type cycle. In Figure 3, h represents a fine grid, and 2h represents the first layer of coarse mesh. grid, 4h represents layer 2 coarse grid. Taking the "V"-shaped loop as an example, if the current calculation layer is a fine grid, the seventeenth step is to select a coarser grid to be calculated; if the current calculation layer is the coarsest grid, then the seventeenth step is to select no need. Calculates a coarser grid.

第十八步:将守恒量修正量插值回更细网格,并判断更细网格是否为细网格;若更细网格是细网格,则返回第六步,进入下一个迭代步的计算;若更细网格不是细网格,则返回第六步,进入更细网格的计算。Step 18: Interpolate the conservation amount correction back to a finer grid, and determine whether the finer grid is a fine grid; if the finer grid is a fine grid, return to step 6 and enter the next iteration step If the finer mesh is not a fine mesh, go back to step 6 and enter into the calculation of the finer mesh.

第十九步:输出结果。Step 19: Output the results.

下面分别对马赫数0.755的NACA0012翼型无粘绕流问题和马赫数0.729的跨声速RAE2822翼型湍流绕流问题进行模拟。The inviscid flow around the NACA0012 airfoil with Mach number 0.755 and the turbulent flow around the transonic RAE2822 airfoil with Mach number 0.729 are simulated as follows.

(1)对马赫数0.755的NACA0012翼型无粘绕流问题进行模拟。(1) Simulate the inviscid flow around the NACA0012 airfoil with a Mach number of 0.755.

图4示意了本发明提供的上述多重网格扰动域更新加速方法的对流动态计算域演化过程。图4中,ηc表示对流动态计算域与预设计算域的网格量之比,Nmax代表达到收敛态所需的总迭代步数。流场以来流条件初始化,扰动从壁面边界产生。在求解初期,对流动态计算域从壁面逐渐扩展至整个预设计算域。在求解中后期,对流动态计算域开始随着求解收敛逐渐从远场向壁面、从上游向下游收缩,最后仅剩近壁区及其下游区域。从图5a所示的多重网格扰动域更新加速方法的对流动态计算域网格量变化曲线可知,本发明提供的上述多重网格扰动域更新加速方法适用于不同网格层数的多重网格,不同网格层数情况下的对流动态计算域演化过程也基本相同。FIG. 4 illustrates the evolution process of the convective dynamic computational domain of the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention. In Fig. 4, η c represents the ratio of the grid size of the convective dynamic computational domain to the preset computational domain, and N max represents the total number of iteration steps required to reach a convergent state. The flow field is initialized with flow conditions and disturbances are generated from the wall boundary. At the beginning of the solution, the convective dynamic computational domain gradually expands from the wall to the entire preset computational domain. In the middle and late stages of the solution, the convective dynamic computational domain begins to gradually shrink from the far field to the wall and from the upstream to the downstream as the solution converges, and finally only the near-wall region and its downstream region remain. It can be seen from the change curve of the grid quantity in the convective dynamic calculation domain of the multi-grid disturbance domain update acceleration method shown in FIG. 5a that the multi-grid disturbance domain update acceleration method provided by the present invention is suitable for multi-grids with different grid layers. , the evolution process of the convective dynamic computational domain is basically the same under the condition of different grid layers.

图5b对比了本发明提供的上述多重网格扰动域更新加速方法与现有全局更新的多重网格方法的收敛曲线。对比结果表明,对于不同层数的多重网格方法,相比于现有方法,本发明提供的上述多重网格扰动域更新加速方法通过减少单个迭代步的计算量和总迭代步数,均可显著提升数值模拟效率。其中,当本发明采用1层粗网格时,相比于现有多重网格方法可节省52.1%的计算时间;当本发明采用2层粗网格时,相比于现有多重网格方法可节省45.6%的计算时间。Fig. 5b compares the convergence curve of the multi-grid disturbance domain update acceleration method provided by the present invention and the existing global update multi-grid method. The comparison results show that, for the multi-grid methods with different layers, compared with the existing methods, the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention can reduce the calculation amount of a single iteration step and the total number of iteration steps. Significantly improve the efficiency of numerical simulation. Among them, when the present invention adopts 1 layer of coarse grid, compared with the existing multi-grid method, 52.1% of the calculation time can be saved; when the present invention adopts 2 layers of coarse grid, compared with the existing multi-grid method 45.6% of computing time can be saved.

(2)对马赫数0.729的跨声速RAE2822翼型湍流绕流问题进行模拟。(2) Simulate the turbulent flow around the transonic RAE2822 airfoil with Mach number 0.729.

图6给出了本发明提供的上述多重网格扰动域更新加速方法的动态计算域演化过程。流场以来流条件初始化,因此,对流动态计算域和粘性动态计算域仍依据壁面建立。对流动态计算域的增大速率远大于粘性动态计算域的增大速率。由图6中1%Nmax的情况可知,当对流动态计算域已覆盖整个预设计算域时,粘性动态计算域仍保持在壁面附近,仅包含近壁区及其下游区域。由图7a给出的多重网格扰动域更新加速方法的对流、粘性两类动态计算域网格量变化曲线可知,本发明提供的上述多重网格扰动域更新加速方法适用于粘性流动,对流动态计算域和粘性动态计算域均遵循先增大后收缩的演化规律。FIG. 6 shows the dynamic calculation domain evolution process of the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention. The flow field is initialized from the flow conditions, therefore, the convective dynamic computational domain and the viscous dynamic computational domain are still established based on the wall surface. The growth rate of the convective dynamic computational domain is much greater than that of the viscous dynamic computational domain. It can be seen from the situation of 1% N max in Figure 6 that when the convective dynamic computational domain has covered the entire preset computational domain, the viscous dynamic computational domain still remains near the wall, including only the near-wall area and its downstream area. 7a, it can be seen from the grid quantity change curves of the two types of dynamic calculation domains, convective and viscous, that the multi-grid disturbance domain update acceleration method provided by the present invention is suitable for viscous flow, convective dynamic Both the computational domain and the viscous dynamic computational domain follow the evolution law of first increasing and then shrinking.

图7b对比了本发明提供的上述多重网格扰动域更新加速方法与现有全局更新的多重网格方法的收敛曲线。对比结果表明,本发明提供的上述多重网格扰动域更新加速方法与现有的多重网格方法一样,均可获得收敛解。得益于对流、粘性两类动态计算域的缩小和仅在局部考虑粘性,本发明提供的上述多重网格扰动域更新加速方法获得了显著的计算效率提升,可节省40.7%的计算时间。Fig. 7b compares the convergence curves of the multigrid disturbance domain update acceleration method provided by the present invention and the existing global update multigrid method. The comparison results show that the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention is the same as the existing multi-grid method, and a convergent solution can be obtained. Benefiting from the reduction of two types of dynamic computational domains of convection and viscosity and considering only local viscosity, the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention obtains a significant increase in computational efficiency and can save 40.7% of computational time.

本发明提供的上述飞行器绕流数值模拟的多重网格扰动域更新加速方法,残差估计仅在各层网格的对流动态计算域内计算残差的无粘项,仅在各层网格的粘性动态计算域内计算残差的粘性项,时间积分也仅在各层网格的对流动态计算域内执行,因此,本发明通过在各层网格采用对流、粘性两类动态计算域,可同时降低多重网格方法模拟飞行器绕流问题时各层网格的计算量,显著提升多重网格方法针对飞行器绕流数值模拟的计算效率,从而有效解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类;其中,对于第一类数据,以静态数据结构存储预设计算域中所有单元的信息;对于第二类数据,以动态数据结构仅存储各网格层对流动态计算域中网格单元的信息,并在每个迭代步及时根据各网格层的对流动态计算域调整第二类数据的存储空间,因此,相比于现有多重网格方法,本发明可有效降低飞行器绕流数值模拟的存储空间需求。针对飞行器流场的细网格,根据扰动传播特性仅将受扰动单元加入对流动态计算域中;通过衡量粘性效应的影响,仅将受粘性效应主导的网格单元加入细网格的粘性动态计算域中;利用求解收敛特征,在不影响收敛速率的前提下,将已收敛单元从对流、粘性两类动态计算域中移除;根据细网格的动态计算域,更新飞行器流场的所有粗网格的对流、粘性两类动态计算域;因此,本发明可始终保持在飞行器流场的所有网格层仅对求解未收敛的受扰单元求解,仅在粘性效应主导单元中考虑粘性,从而可有效避免现有多重网格方法在飞行器绕流数值模拟中的无效计算,解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。In the multi-grid disturbance domain update acceleration method for the numerical simulation of flow around an aircraft provided by the present invention, the residual estimation only calculates the inviscid term of the residual in the In the dynamic calculation domain, the viscous term of the residual is calculated, and the time integration is only performed in the convective dynamic calculation domain of each layer grid. Therefore, the present invention adopts two types of dynamic calculation domains, convective and viscous, in each layer grid. When the grid method simulates the flow around the aircraft, the calculation amount of each layer of grids significantly improves the calculation efficiency of the multi-grid method for the numerical simulation of the flow around the aircraft, thus effectively solving the problem that the calculation efficiency of numerical simulation is difficult to meet the increasing practical application of aircraft aerodynamic design. question of needs. The data to be stored is divided into two categories: the inherent information of all units and the information related to solving and updating; wherein, for the first type of data, the information of all units in the preset computational domain is stored in a static data structure; for the second type of data data, only the information of grid cells in the convective dynamic computing domain of each grid layer is stored in a dynamic data structure, and the storage space of the second type of data is adjusted according to the convective dynamic computing domain of each grid layer in time at each iteration step. Therefore, , compared with the existing multi-grid method, the invention can effectively reduce the storage space requirement of the numerical simulation of the flow around the aircraft. For the fine grid of the aircraft flow field, only the disturbed elements are added to the convective dynamic calculation domain according to the disturbance propagation characteristics; by weighing the influence of the viscous effect, only the grid elements dominated by the viscous effect are added to the viscous dynamic calculation of the fine grid In the domain; using the solution convergence feature, without affecting the convergence rate, the converged elements are removed from the two types of dynamic calculation domains of convection and viscosity; Convective and viscous two types of dynamic computational domains of the grid; therefore, the present invention can always keep all grid layers in the flow field of the aircraft to solve only the unconverged disturbed element, and only consider the viscosity in the viscous effect dominant element, thus The ineffective calculation of the existing multi-grid method in the numerical simulation of the flow around the aircraft can be effectively avoided, and the problem that the numerical simulation calculation efficiency is difficult to meet the increasing practical application requirements of the aerodynamic design of the aircraft can be solved.

显然,本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样,倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内,则本发明也意图包含这些改动和变型在内。It will be apparent to those skilled in the art that various modifications and variations can be made in the present invention without departing from the spirit and scope of the invention. Thus, provided that these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include these modifications and variations.

Claims (7)

1. A multiple grid disturbance domain updating and accelerating method for aircraft bypass flow numerical simulation is characterized by comprising the following steps:
s1: reading data, including fine grids of an aircraft flow field, a preset calculation domain, boundary conditions and calculation settings;
s2: generating a coarse grid with a specified number of layers by thinning the fine grid;
s3: initializing flow fields in preset calculation domains on the coarse grid and the fine grid according to an incoming flow condition or a given flow field;
s4: establishing a convection dynamic calculation domain and a viscosity dynamic calculation domain according to a flow field initialization mode;
s5: dividing data to be stored for solving all layers of grids into two types; the first kind of data is inherent information of a unit, and a static data structure is adopted to store information of all grid units in the preset calculation domain, wherein the information comprises grid coordinates and flow field variables; the second kind of data is information related to solving and updating, and only the information of grid cells in the convection dynamic calculation domain is stored by adopting a dynamic data structure, wherein the information comprises a conservative updating quantity and a local iteration step length;
S6: assigning a value to the conservation quantity of the boundary virtual grid of the current layer grid according to the type of the boundary condition;
s7: dividing residual terms of a flow control equation into two types of non-sticky terms and sticky terms; wherein the inviscid terms comprise forcing functions of the convection flux and the coarse grid, the inviscid terms being computed in the convection dynamics computation domain; the viscosity term comprises a viscosity flux and turbulence model equation source term, the viscosity term being calculated in a viscosity dynamics computation domain;
s8: in the convection dynamic calculation domain, solving a conservative updating quantity and updating a flow field variable;
s9: judging whether the current grid layer is a fine grid or not; if yes, go to step S10; if not, go to step S21;
s10: judging whether the fine grid reaches a convergence condition; if not, go to step S11; if yes, go to step S23;
s11: judging whether all boundary units of the convection dynamic calculation domain are subjected to inviscid disturbance one by one; if yes, go to step S12; if not, returning to the step S11, and judging the next boundary unit of the convection dynamic calculation domain; after traversing all the boundary cells of the convection dynamic calculation domain, executing step S13;
s12: measuring the propagation direction of the inviscid disturbance, bringing the cells positioned in the propagation direction in the adjacent cells of the boundary cell into the convection dynamic calculation domain, returning to the step S11, and judging the next boundary cell in the convection dynamic calculation domain;
S13: judging whether all boundary units of the convection dynamic calculation domain meet four conditions of converged, in compressible flow, at the most upstream and no longer influenced by other units in the convection dynamic calculation domain; if yes, go to step S14; if not, returning to the step S13, and judging the next boundary unit of the convection dynamic calculation domain; after traversing all the boundary cells of the convection dynamic calculation domain, executing step S16;
s14: removing the boundary unit from the convection dynamic calculation domain, and judging whether the boundary unit exists in the viscosity dynamic calculation domain; if yes, go to step S15; if not, returning to the step S13, and judging the next boundary unit of the convection dynamic calculation domain;
s15: removing the boundary cell from the viscous dynamic calculation domain, returning to step S13, and determining the next boundary cell in the convective dynamic calculation domain;
s16: judging whether all boundary units of the viscosity dynamic calculation domain are dominated by viscosity effect one by one; if yes, go to step S17; if not, returning to the step S16, and judging the next boundary unit of the viscosity dynamic calculation domain; after traversing all boundary cells of the sticky dynamic computation domain, executing step S18;
S17: bringing the cell located in the convection dynamic calculation domain in the immediate cell of the boundary cell into the viscous dynamic calculation domain, returning to step S16, and determining the next boundary cell of the viscous dynamic calculation domain;
s18: updating the convection dynamic calculation domain and the viscosity dynamic calculation domain of all layers of the coarse grids according to the convection dynamic calculation domain and the viscosity dynamic calculation domain of the fine grids;
s19: reallocating storage spaces of second-class data in all the coarse grids and the fine grids according to the convection dynamic calculation domain;
s20: interpolating the conservative quantity to a coarser grid, returning to the step S6, and entering the calculation of the coarser grid;
s21: judging whether a coarser grid needs to be calculated; if yes, return to step S20; if not, go to step S22;
s22: interpolating the conservative quantity correction back to the finer grid, and judging whether the finer grid is a fine grid or not; if yes, returning to the step S6, and entering the calculation of the next iteration step; if not, returning to the step S6 to enter the calculation of the finer grid;
s23: and outputting the result.
2. The method for accelerating updating of multiple grid disturbance domains according to claim 1, wherein step S3, initializing flow fields in preset computation domains on the coarse grid and the fine grid according to an incoming flow condition or according to a given flow field, specifically comprises:
Initializing according to incoming flow conditions, and assigning the conservative quantities of all grid units in a preset calculation domain on the coarse grid and the fine grid as incoming flow values;
and initializing according to a given flow field, assigning the conservation quantities of all grid units in a preset calculation domain on the fine grid to be given flow field values, and obtaining the conservation quantities of all grid units in the preset calculation domain on the coarse grid according to the conservation quantity interpolation of all grid units on the fine grid.
3. The method for accelerating updating of multiple grid disturbance domains according to claim 1, wherein step S4 is to establish a convection dynamic calculation domain and a viscosity dynamic calculation domain according to a flow field initialization method, and specifically includes:
initializing according to incoming flow conditions, regarding a convection dynamic calculation domain, taking 1 layer of cells adjacent to a wall surface as initial cells by the coarsest layer of grids, taking L layers of adjacent cells of the wall surface as initial cells by the fine grids, wherein L satisfies the following conditions:
L=2l(1)
wherein l represents the number of coarse grid layers, is a positive integer, and is given by the calculation setting read in step S1; generating initial units of the coarse layer grids according to the initial units of the fine grids; for the viscous dynamic calculation domain, all grid layers take 1 layer of units adjacent to the wall surface as initial units;
According to initialization of a given flow field, for a convection dynamic calculation domain, an initial unit of a fine grid is an interfered unit with flow characteristics in the given flow field inconsistent with incoming flow conditions, and the following conditions are met:
||W-W||/||ΔW(1)||maxa,c(2)
wherein W represents a conservative amount; wA conservation quantity representing an incoming flow condition; | Δ W(1)||maxRepresenting the maximum value of the conservative updating quantity of all units in the convection dynamic calculation domain in the 1 st iteration step;a,cthe newly added threshold value of convection is taken as 10-4a,c≤10-6(ii) a For the viscous dynamic calculation domain, the initial unit of the fine grid is a viscous effect dominant unit in the fine grid convection dynamic calculation domain, and the following conditions are met:
Figure FDA0002534371270000041
Figure FDA0002534371270000042
wherein Ψ represents a viscous effect measurement parameter as a ratio of viscous disturbance to non-viscous disturbance mass flow;
Figure FDA0002534371270000043
a viscosity effect measurement parameter representing the 1 st iteration; i, J and K respectively represent grid directions;
Figure FDA0002534371270000044
represents the spectral radius of the convection flux Jacobian matrix along the i direction,
Figure FDA0002534371270000045
representing the spectral radius of the viscous flux Jacobian matrix along the i direction;a,vindicating the newly added threshold value of viscosity, take 10-3a,v≤10-2(ii) a And the convection dynamic calculation domains of all the coarse grid layers are generated according to the convection dynamic calculation domain of the fine grid, and the viscosity dynamic calculation domains of all the coarse grid layers are generated according to the viscosity dynamic calculation domain of the fine grid.
4. The multiple mesh perturbation domain updating acceleration method according to claim 3 characterized in that step S7, the residual terms of the flow control equation are divided into two types, non-sticky terms and sticky terms; wherein the inviscid terms comprise forcing functions of the convection flux and the coarse grid, the inviscid terms being computed in the convection dynamics computation domain; the viscosity term comprises a viscosity flux and a turbulence model equation source term, and is calculated in a viscosity dynamic calculation domain, and specifically comprises the following steps:
The flow control equation is expressed as:
Figure FDA0002534371270000046
wherein t represents time, FcDenotes the convective flux, FvDenotes viscous flux, QFForcing function, Q, representing a coarse gridTRepresenting the source term of the turbulence model equation, | Ω | representing the volume of the grid cell, NfThe number of the unit surfaces is shown, and Delta S represents the area of the unit surfaces; the right-hand term of equation (5) is collectively referred to as the residual term of the flow control equation,
Figure FDA0002534371270000051
is a non-sticky item, and is provided with a sticky part,
Figure FDA0002534371270000052
is a viscosity term; solving the non-sticky terms of the residual terms in the flow dynamics computation domain and solving the sticky terms of the residual terms in the sticky dynamics computation domain.
5. The method as claimed in claim 4, wherein the step S11 of determining whether all boundary cells of the convective dynamic computation domain have been disturbed without sticking one by one includes:
judging whether all boundary units of the convection dynamic calculation domain meet the following requirements one by one:
||ΔW||/||ΔW(1)||maxa,c(6)
wherein, | Δ W | | | represents a modulus of a constant update quantity of a boundary unit of the convection dynamic calculation domain in the current iteration step;
step S12, measuring the propagation direction without viscous disturbance, and incorporating the cell located in the propagation direction in the cells immediately adjacent to the boundary cell into the convection dynamic calculation domain, which specifically includes:
The inviscid disturbance is transmitted at the sound velocity relative to the flow, and the transmission speed of the inviscid disturbance is positive along the transmission direction of the inviscid disturbance; let q denote the unit direction vector, then the propagation of the inviscid perturbation along the q direction is expressed as:
u·q+a>0 (7)
wherein u represents a flow velocity vector, and a represents a sound velocity; and q is a unit vector of a grid point of the boundary unit of the convection dynamic calculation domain pointing to a certain grid point of the boundary unit, and if q satisfies the formula (7), an adjacent unit sharing the grid point with the boundary unit of the convection dynamic calculation domain is brought into the convection dynamic calculation domain.
6. The multiple-grid-disturbance-domain updating acceleration method according to claim 5, wherein step S13 is implemented to determine whether four conditions of converged, in compressible flow, in most upstream, and no longer influenced by other units in the convective dynamic computation domain are simultaneously satisfied for all boundary cells of the convective dynamic computation domain one by one; if yes, go to step S14; if not, returning to the step S13, and judging the next boundary unit of the convection dynamic calculation domain; the method specifically comprises the following steps:
s131: judging whether all boundary units of the convection dynamic calculation domain meet the following requirements one by one:
||ΔW||/||ΔW(1)||maxd(8)
wherein,dindicating a deletion threshold, take d≤10-7(ii) a If yes, go to step S132; if not, returning to the step S131, and judging a next boundary unit of the convection dynamic calculation domain;
s132: judging whether the Mach number of the boundary unit is larger than 0.3; if yes, go to step S133; if not, returning to the step S131, and judging a next boundary unit of the convection dynamic calculation domain;
s133: judging whether the boundary unit satisfies the following conditions:
Figure FDA0002534371270000061
wherein, q is a unit vector of grid center of the convection dynamic calculation domain boundary unit pointing to the adjacent cell center of the boundary unit, thetadRepresenting the tolerance angle of the upstream unit, and taking theta more than or equal to 5 degreesdLess than or equal to 10 degrees; if yes, go to step S134; if not, returning to the step S131, and judging a next boundary unit of the convection dynamic calculation domain;
s134: for supersonic non-viscous flow, step S14 is directly executed without further judgment;
for subsonic and viscous flows, it is determined whether the boundary cell satisfies:
||ΔW||+||Δ(ΔW)||<d(10)
Figure FDA0002534371270000062
wherein, | Δ (Δ W) | | represents a module value of the boundary unit conservative update quantity affected by the adjacent unit in the convection dynamic calculation domain of the boundary unit; Δ t denotes the iteration step, CCFLNumber of CFLs, Δ R, representing a time advance formatiRepresenting the influence of the residual error items of the adjacent units on the boundary units along the i direction in the flow dynamic calculation domain;
For subsonic inviscid cells, the effect Δ R in the i direction of the residual terms of the adjacent cells to the boundary cell in the convection dynamics computation domainiExpressed as:
Figure FDA0002534371270000063
wherein, Δ FcRepresenting the amount of change in convective flux; the subscript i +1 represents the immediate vicinity of the boundary cell of the convective dynamics computation domain along the positive i direction, the subscript i-1 represents the immediate vicinity of the boundary cell of the convective dynamics computation domain along the negative i direction, Δ Si+1/2The area of a boundary unit of the convection dynamic calculation domain along a positive i direction unit surface is represented; delta Si-1/2The area of a boundary unit of the convection dynamic calculation domain along a negative i direction unit surface is represented;
for sticky cells, the effect Δ R in the i direction of the residual term of the immediate cell on the boundary cell in the flow dynamics computation domainiExpressed as:
Figure FDA0002534371270000071
if yes, go to step S14; if not, the process returns to step S131 to determine the next boundary unit of the convection dynamic calculation domain.
7. The method for accelerating updating of multiple grid disturbance domains according to claim 6, wherein step S16, determining whether all boundary cells of the sticky dynamic computation domain are dominated by a sticky effect one by one, specifically comprises:
judging whether all boundary units of the viscosity dynamic calculation domain meet the following conditions one by one:
Figure FDA0002534371270000072
Figure FDA0002534371270000073
Wherein phi represents the scaling factor of the newly added viscosity threshold, | | Δ W | | magnetismmaxRepresenting the maximum value of the conservative update quantity of all units in the convection dynamic calculation domain in the current iteration step (| | Δ W | | countmax)minRepresents the minimum value of the maximum value of the conservation constant updating quantity of the single step from the 1 st iteration step to the current iteration step.
CN202010528125.7A 2020-06-11 2020-06-11 A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft Active CN111859529B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010528125.7A CN111859529B (en) 2020-06-11 2020-06-11 A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010528125.7A CN111859529B (en) 2020-06-11 2020-06-11 A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft

Publications (2)

Publication Number Publication Date
CN111859529A true CN111859529A (en) 2020-10-30
CN111859529B CN111859529B (en) 2022-04-08

Family

ID=72987229

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010528125.7A Active CN111859529B (en) 2020-06-11 2020-06-11 A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft

Country Status (1)

Country Link
CN (1) CN111859529B (en)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113239462A (en) * 2021-05-25 2021-08-10 江苏普旭科技股份有限公司 Simulation method for aircraft turbulent environment simulation
CN113962030A (en) * 2021-12-20 2022-01-21 北京航空航天大学 An Overlapping Grid Disturbance Domain Updating Method for Aircraft Multibody Separation Simulation
CN114218878A (en) * 2022-02-17 2022-03-22 北京航空航天大学 Dynamic grid disturbance domain updating method for aircraft maneuvering process simulation
CN116127877A (en) * 2023-04-04 2023-05-16 中国空气动力研究与发展中心计算空气动力研究所 Acceleration method, device, terminal equipment and storage medium for multiple grids
CN116757110A (en) * 2023-06-05 2023-09-15 浙江大学 Full-speed domain numerical simulation nonlinear numerical iteration solving method for complex aircraft
CN117669429A (en) * 2024-01-31 2024-03-08 陕西空天信息技术有限公司 Fluid simulation method and device for rotary machine and computer storage medium

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120245903A1 (en) * 2011-03-23 2012-09-27 Desktop Aeronautics, Inc. Generating inviscid and viscous fluid flow simulations over a surface using a quasi-simultaneous technique
CN106354903A (en) * 2016-08-18 2017-01-25 中国人民解放军国防科学技术大学 Determination method of the computational domain outer boundary for solving the steady circumfluence value of flying objects
CN108563843A (en) * 2018-03-26 2018-09-21 北京航空航天大学 The disturbance region update method of steady compressible flowing
CN108595788A (en) * 2018-04-05 2018-09-28 西北工业大学 A kind of flow field Accelerated Convergence Method based on mode multi grid
CN110852005A (en) * 2019-10-21 2020-02-28 北京理工大学 A Numerical Simulation Method for Adaptively Expanding Computational Domain for Massively Parallel Computing

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120245903A1 (en) * 2011-03-23 2012-09-27 Desktop Aeronautics, Inc. Generating inviscid and viscous fluid flow simulations over a surface using a quasi-simultaneous technique
CN106354903A (en) * 2016-08-18 2017-01-25 中国人民解放军国防科学技术大学 Determination method of the computational domain outer boundary for solving the steady circumfluence value of flying objects
CN108563843A (en) * 2018-03-26 2018-09-21 北京航空航天大学 The disturbance region update method of steady compressible flowing
CN108595788A (en) * 2018-04-05 2018-09-28 西北工业大学 A kind of flow field Accelerated Convergence Method based on mode multi grid
CN110852005A (en) * 2019-10-21 2020-02-28 北京理工大学 A Numerical Simulation Method for Adaptively Expanding Computational Domain for Massively Parallel Computing

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SHUYAO HU 等: "Disturbance region update method for steady compressible flows", 《COMPUTER PHYSICS COMMUNICATIONS》 *
SHUYAO HU 等: "Zonal disturbance region update method for steady compressible viscous flows", 《COMPUTER PHYSICS COMMUNICATIONS》 *
李杰 等: "结构化多重网格粘性流场数值模拟", 《计算机仿真》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113239462A (en) * 2021-05-25 2021-08-10 江苏普旭科技股份有限公司 Simulation method for aircraft turbulent environment simulation
CN113239462B (en) * 2021-05-25 2025-06-17 江苏普旭科技股份有限公司 A simulation method for simulating aircraft turbulent environment
CN113962030A (en) * 2021-12-20 2022-01-21 北京航空航天大学 An Overlapping Grid Disturbance Domain Updating Method for Aircraft Multibody Separation Simulation
CN114218878A (en) * 2022-02-17 2022-03-22 北京航空航天大学 Dynamic grid disturbance domain updating method for aircraft maneuvering process simulation
CN116127877A (en) * 2023-04-04 2023-05-16 中国空气动力研究与发展中心计算空气动力研究所 Acceleration method, device, terminal equipment and storage medium for multiple grids
CN116127877B (en) * 2023-04-04 2023-09-22 中国空气动力研究与发展中心计算空气动力研究所 Acceleration method, device, terminal equipment and storage medium for multiple grids
CN116757110A (en) * 2023-06-05 2023-09-15 浙江大学 Full-speed domain numerical simulation nonlinear numerical iteration solving method for complex aircraft
CN116757110B (en) * 2023-06-05 2025-04-25 浙江大学 Full-speed domain numerical simulation nonlinear numerical iteration solving method for complex aircraft
CN117669429A (en) * 2024-01-31 2024-03-08 陕西空天信息技术有限公司 Fluid simulation method and device for rotary machine and computer storage medium
CN117669429B (en) * 2024-01-31 2024-05-10 陕西空天信息技术有限公司 Fluid simulation method and device for rotary machine and computer storage medium

Also Published As

Publication number Publication date
CN111859529B (en) 2022-04-08

Similar Documents

Publication Publication Date Title
CN111859529B (en) A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft
CN111859530B (en) Iterative propulsion disturbance domain updating method for aircraft dynamic aerodynamic characteristic simulation
CN113850008B (en) Adaptive grid disturbance domain update acceleration method for aerodynamic prediction of aircraft
CN111651831B (en) Update calculation method of partitioned disturbance domain for steady viscous compressible flow around an aircraft
Biancolini et al. Static aeroelastic analysis of an aircraft wind-tunnel model by means of modal RBF mesh updating
CN113609598B (en) RANS/LES disturbance domain update method for aerodynamic simulation of aircraft
CN111079228B (en) An aerodynamic shape optimization method based on flow field prediction
CN109190233B (en) A structure topology optimization method
CN113392472B (en) OpenMP parallel disturbance domain updating method for aircraft aerodynamic characteristic simulation
Pan et al. Computation of incompressible flows with immersed bodies by a simple ghost cell method
CN108563843A (en) The disturbance region update method of steady compressible flowing
Shershnev et al. HyCFS, a high-resolution shock capturing code for numerical simulation on hybrid computational clusters
CN110032756A (en) Flow boundary layer numerical analysis method based on the transformation of stream function fractional coordinates system
CN115828782A (en) Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm
Lu et al. Flow simulation system based on high order space-time extension of flux reconstruction method
JP3587827B2 (en) Airfoil performance estimation method and airfoil performance estimation program
CN115017604B (en) An efficient disturbance domain propulsion method for compressible/incompressible mixed flows
CN114218878B (en) A dynamic grid disturbance domain update method for aircraft maneuvering process simulation
CN112182771B (en) Data processing method based on numerical simulation, storage medium and electronic device
Remigi et al. Evaluation of tandem transonic compressor performances with mixing plane and metric-based mesh adaptation
CN114757115A (en) Method, device, equipment and storage medium for numerical determination of flow field of compressor
Tota et al. Meshfree Euler solver using local radial basis functions for inviscid compressible flows
Kolesnik et al. Testing of various schemes with quasi-one-dimensional reconstruction of gasdynamic variables in the case of unstructured-grid calculations
Paraschivoiu et al. Multi-domain multi-model formulation for compressible flows-Conservative interface coupling and parallel implicit solvers for 3D unstructured meshes
Perić Finite-Volume Methods for Navier-Stokes Equations

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