CN115017604A - Efficient disturbance domain propulsion method for compressible/incompressible mixed flow - Google Patents
Efficient disturbance domain propulsion method for compressible/incompressible mixed flow Download PDFInfo
- Publication number
- CN115017604A CN115017604A CN202210469835.6A CN202210469835A CN115017604A CN 115017604 A CN115017604 A CN 115017604A CN 202210469835 A CN202210469835 A CN 202210469835A CN 115017604 A CN115017604 A CN 115017604A
- Authority
- CN
- China
- Prior art keywords
- domain
- grid
- dynamic
- parameters
- compressible
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 238000004364 calculation method Methods 0.000 claims abstract description 89
- 238000007781 pre-processing Methods 0.000 claims abstract description 20
- 230000008859 change Effects 0.000 claims abstract description 10
- 230000004907 flux Effects 0.000 claims description 15
- 230000003247 decreasing effect Effects 0.000 claims description 9
- 230000000694 effects Effects 0.000 claims description 7
- 230000010354 integration Effects 0.000 claims description 5
- 238000011144 upstream manufacturing Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 abstract description 4
- 239000012530 fluid Substances 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 9
- 238000005192 partition Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 5
- 230000009467 reduction Effects 0.000 description 4
- 238000004088 simulation Methods 0.000 description 3
- 238000005293 physical law Methods 0.000 description 2
- 230000035939 shock Effects 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Analysis (AREA)
- Evolutionary Computation (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Aviation & Aerospace Engineering (AREA)
- Computational Mathematics (AREA)
- Automation & Control Theory (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Physics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The disclosure relates to the technical field of computational fluid mechanics, and provides a high-efficiency disturbance domain propulsion method for compressible/non-compressible mixed flow. The efficient disturbance domain updating method comprises the following steps: data reading and flow field initialization; establishing an initial dynamic calculation domain according to a flow field initialization mode, and performing residual estimation and time integral iterative solution based on the dynamic calculation domain; updating the dynamic calculation domain; based on the relation between the global residual and a given convergence residual threshold, if yes, judging whether the calculation is finished, if yes, outputting, and if not, skipping to solving the residual and time integral iteration; and outputting a calculation result. The invention solves the problems of change of the characteristic value of the Jacobian matrix and inapplicability of the criterion of the low Mach number area, provides a viscous leading judgment basis suitable for solving the compressible/non-compressible flow field by preprocessing, avoids a large amount of invalid calculations and greatly improves the calculation efficiency.
Description
Technical Field
The disclosure relates to the technical field of computational fluid mechanics, in particular to a high-efficiency disturbance domain propulsion method for compressible/incompressible mixed flow.
Background
In the aerospace field, many important problems are compressible/incompressible mixed flow, and when a compressible calculation method is directly adopted to solve the compressible/incompressible mixed flow, the difference of characteristic values causes serious pathological condition numbers, and the efficiency is reduced and the calculation is difficult. The traditional low-speed flow is usually expanded based on a semi-implicit method of a pressure coupling equation set, the efficiency of the expansion is higher when the low-speed problem is calculated, but the pressure and speed conservation of the separated solving is poor compared with the density base, and the traditional shock wave capturing format is difficult to directly popularize. While compressible flow solving methods characterized by a time propulsion method and a shock wave capture format have been developed, when low-speed flow is directly solved, the system rigidity is too high and the convergence speed is too low due to too large difference between the sound velocity and the flow velocity. The advent of preprocessing methods has made the use of compressible equations to solve for compressible/non-compressible mixed flow fields an optimal approach.
However, although the preprocessing method improves the convergence rate of the compressible equation in solving the compressible/non-compressible flow field, in practical engineering application, higher requirements are put forward on the calculation efficiency of the compressible equation. The disturbance domain updating method is a novel numerical simulation acceleration method provided for compressible flow solving, and effectively improves the calculation efficiency of numerical simulation. However, the traditional disturbance domain updating method cannot be applied to the solution of the preprocessing method to the compressible/non-compressible flow, and the defects are as follows:
firstly, the preprocessing method changes the eigenvalue of the convective Jacobian matrix, the eigenvalue is the local speed and the local sound velocity which are used as disturbance propagation criteria in the traditional disturbance domain updating method, and the change influences a plurality of criteria in the disturbance domain updating method; secondly, the compressible/incompressible mixed flow is usually low in Mach number in most areas, the density and temperature change is small, and the judgment failure can be caused by taking the total temperature or the total enthalpy as a viscosity partition criterion.
Therefore, a method for effectively improving the calculation efficiency of compressible/non-compressible flow and serving numerical simulation calculation in the related field is lacking at present.
Disclosure of Invention
In view of this, the present disclosure provides a high-efficiency disturbance domain propulsion method for compressible/incompressible mixed flow, so as to solve the technical problem that the calculation efficiency is greatly improved while the error caused by the finite of the calculation domain is not affected in the prior art.
The present disclosure provides a method for propelling a high-efficiency disturbance domain of compressible/non-compressible mixed flow, comprising:
s1, reading data and initializing a flow field;
s2, establishing an initial dynamic calculation domain according to the flow field initialization mode, wherein the dynamic calculation domain comprises a convection dynamic domain and a viscosity dynamic domain;
s3, in the dynamic calculation domain, solving the residual error of the compressible Navier-Stokes equation at the current time step, and carrying out time integration iteration to solve the grid parameters;
s4, updating the dynamic calculation domain based on the residual error and the grid parameters obtained in S3;
s5, judging whether the global residual is smaller than a given convergence residual threshold, if so, finishing the calculation and executing S6, and if not, returning to S3;
s6 outputs the calculated values of the mesh parameters.
Further, the S1 includes:
reading in computational grids, boundary conditions, grid parameters and computational settings;
based on the grid parameters, initializing the flow field according to the incoming flow conditions or the given flow field, and determining initial values of physical quantities of which boundaries and domains are not defined before calculation and formula parameters in the selected physical model, wherein the formula parameters in the physical model comprise a gas constant, a viscosity coefficient and a turbulence model coefficient.
Further, the S2 includes:
when the dynamic viscous dynamic domain is initialized according to the incoming flow conditions, all grid parameters are set as incoming flow parameters, the convective dynamic domain is established according to the wall surface conditions, and only 1 layer of grid units adjacent to the wall surface is taken as an initial viscous dynamic domain;
when the method is initialized according to a given flow field, all grid parameters are set according to the given flow field, the viscous dynamic domain is established according to the difference between the given flow field and an incoming flow, grid units which are not consistent with the flowing characteristics of the incoming flow are used as an initial convection dynamic domain, and a viscous leading grid unit in the convection dynamic domain is selected as a viscous dynamic domain.
Further, the S3 includes:
initializing all given grid parameters based on a flow field, solving a compressible Navier-Stokes equation, selecting a preprocessed Roe format, calculating convection flux, and obtaining a residual error of a current time step;
and based on the residual error result, carrying out time integral iterative solution, and updating to obtain the grid parameters of the current time step.
Further, the S4 includes:
based on the grid parameters of the current time step, judging whether the convection dynamic domain needs to be increased or not according to the speed and the sound velocity after the preprocessing, and then judging whether the convection dynamic domain needs to be decreased or not by combining the convergence condition of the residual error;
and judging whether the current grid unit is dominated by a viscosity effect or not, if so, adding all the adjacent grid units of the current grid unit to a viscosity dynamic domain, and if not, deleting the current grid unit from the viscosity dynamic domain.
Further, the increasing and decreasing of the convection dynamic domain in S4 specifically includes:
judging whether the current grid unit is disturbed without adhesion by measuring the grid parameter updating amount calculated in S3, if so, determining the adjacent grid unit which can be influenced, and adding all the adjacent grid units into the convection dynamic domain, wherein for the unit in subsonic flow, all the adjacent grid units can be influenced by the disturbance; for the unit in supersonic flow, determining the propagation direction of disturbance according to the characteristic value of convection flux, and adding the grid unit at the downstream of the propagation direction into a convection dynamic domain;
judging whether the current grid unit is converged, if so, judging whether the current grid unit is positioned at the most upstream of the convection dynamic domain, if so, judging whether other grid units in the convection dynamic domain no longer influence the current grid unit, and if not, deleting the current grid unit from the convection dynamic domain.
Further, the preprocessed speed and sound velocity are obtained based on the local speed and sound velocity, and specifically include:
obtaining a preprocessed speed and a preprocessed sound speed based on the local speed u and the sound speed c;
obtaining a preprocessed characteristic value based on the preprocessed speed and the sound velocity,
the calculation formulas of the speed and the sound velocity after the pretreatment are respectively as follows:
λ＝[u,u,u,u′+c′,u′-c′]
U r 2 ＝min(c 2 ,max(|u| 2 ,K|u| 2 ))
wherein, λ is a characteristic value used in calculating the convection flux, U is a local speed, U 'is a preprocessed speed, c is a sound velocity, c' is a preprocessed sound velocity, α and β are parameters used for solving the preprocessed characteristic value, and U is a local speed r For the pre-processing of the parameters, given by the pre-processing method, ρ p ,ρ T Respectively density, density derivative to pressure, density derivative to temperature, calculated from flow field parameters, C p K is a small value (10) for the specific heat constant at constant pressure -5 )。
Further, the increasing and decreasing of the sticky dynamic domain in S4 specifically includes:
dividing the threshold of the total enthalpy change amount by the Mach number square as a viscosity judgment basis, and increasing and reducing the viscosity dynamic domain, specifically as follows:
if |1-h 0,∞ /h 0 |/Ma 2 ≥ε d,v Then add all immediately adjacent grid cells of the current grid cell to the sticky dynamic field;
if |1-h 0,∞ /h 0 |/Ma 2 <ε d,v If yes, deleting the current grid unit from the sticky dynamic domain;
wherein h is 0,∞ Total enthalpy of incoming flow, h 0 The total enthalpy of the current grid unit, Ma is the Mach number of the current grid unit, epsilon d,v A threshold is determined for a given viscosity domain.
Further, the grid parameters include:
density, velocity, pressure and temperature.
Compared with the prior art, the beneficial effects of this disclosure are:
1. the invention establishes a disturbance propagation calculation formula of a compressible/incompressible disturbance domain updating method based on the preprocessed local sound velocity and the preprocessed local speed, thereby establishing criteria of the compressible/incompressible disturbance domain propagation and influence, solving the problems of change of the characteristic value of the Jacobian matrix and inapplicability of the criterion of a low Mach number area, and contracting a calculation domain;
2. the viscous leading judgment basis suitable for solving the compressible/non-compressible flow field through pretreatment is provided, the threshold of the total enthalpy change is divided by the square of the Mach number to serve as the judgment threshold of the viscous partition, and the judgment threshold of the viscous partition serves as the judgment threshold of the viscous partition, so that the dependency of the viscous partition on the Mach number when the total enthalpy change is used as a criterion is avoided;
3. a large amount of invalid calculations are avoided, and the calculation efficiency is greatly improved.
Drawings
In order to more clearly illustrate the technical solutions of the present disclosure, the drawings needed for the description of the embodiments or the prior art will be briefly described below, and it is apparent that the drawings in the following description are only some embodiments of the present invention, and that other drawings can be obtained by those skilled in the art without inventive efforts.
FIG. 1 is a schematic diagram of a high-efficiency turbulent field propulsion method for mixed flow under pressure/non-pressure provided by the present invention;
FIG. 2 is a schematic illustration of the ratio of convective dynamic field to viscous-dominated region provided by the present invention;
FIG. 3a is a schematic diagram illustrating a comparison of calculation results of a perturbation domain updating method without using compressible/non-compressible flow according to an embodiment of the present invention; FIG. 3b is a schematic diagram illustrating a comparison of calculation results of a perturbation domain updating method using compressible/non-compressible flow according to an embodiment of the present invention;
fig. 4 is a schematic diagram comparing convergence situations of a perturbation domain updating method for determining whether to use compressible/non-compressible flow according to an embodiment of the present invention.
Detailed Description
In the following description, for purposes of explanation and not limitation, specific details are set forth such as particular system structures, techniques, etc. in order to provide a thorough understanding of the disclosed embodiments. However, it will be apparent to one skilled in the art that the present disclosure may be practiced in other embodiments that depart from these specific details. In other instances, detailed descriptions of well-known systems, devices, circuits, and methods are omitted so as not to obscure the description of the present disclosure with unnecessary detail.
An efficient turbulent field propulsion method of a compressible/incompressible mixing flow according to the present disclosure will be described in detail with reference to the accompanying drawings.
Fig. 1 is a schematic diagram of an efficient disturbance domain propulsion method for a compressible/incompressible mixed flow provided by the present invention.
As shown in fig. 1, the efficient perturbation domain updating method includes:
and S1, reading data and initializing a flow field.
S1 includes:
s11, reading in the computational grid, boundary conditions, and computational settings.
And S12, based on the grid parameters, carrying out flow field initialization according to the incoming flow conditions or the given flow field, and determining initial values of physical quantities of which boundaries and domains are not defined before calculation and formula parameters in the selected physical model.
And initializing a flow field based on the grid parameters. The flow field initialization can determine initial values of physical parameters of which boundaries and domains are not defined before calculation and formula parameters in a selected physical model, wherein the formula parameters in the physical model comprise gas constants, viscosity coefficients and turbulence model coefficients.
And S2, establishing an initial dynamic calculation domain according to the flow field initialization mode, wherein the dynamic calculation domain comprises a convection dynamic domain and a viscosity dynamic domain.
And S21, when initializing according to the inflow condition, setting all grid parameters as the inflow parameters, establishing a convection dynamic domain according to the wall condition, and simultaneously taking only 1 layer of grid cells adjacent to the wall as an initial viscous dynamic domain.
When initializing according to the inflow condition, setting all grid parameters as the inflow parameters, establishing a convection dynamic domain according to the wall condition, generally taking 10 layers of grid cells adjacent to the wall boundary as initial grid cells of the convection dynamic domain, and taking the 1 st layer of grid cells adjacent to the wall boundary as initial viscous dynamic domain in the viscous dynamic domain.
S22, when initializing according to the given flow field, all grid parameters are set according to the given flow field, a viscous dynamic domain is established according to the difference between the given flow field and the incoming flow, grid cells which are not in accordance with the flowing characteristics of the incoming flow are used as an initial convection dynamic domain, and a viscous leading grid cell in the convection dynamic domain is selected as a viscous dynamic domain.
And S3, in the dynamic calculation domain, solving the residual error of the compressible Navier-Stokes equation at the current time step, and carrying out time integration iteration to solve the grid parameters.
And solving the residual error of the compressible Navier-Stokes equation at the current time step, performing time advancing, and updating the grid parameters of the dynamic calculation domain.
S3 includes:
and S31, initializing all given grid parameters based on the flow field, solving a compressible Navier-Stokes equation, selecting a preprocessed Roe format, calculating the convection flux, and obtaining the residual error of the current time step.
And initializing parameters of all given grid units based on the flow field, and performing residual estimation to obtain updated grid parameters. Firstly, an initial value of the grid parameter is given, and the updating amount of the grid parameter is obtained after residual estimation. The residual is obtained by the following expression:
residual error is initial value-next value
And S32, based on the residual error result, carrying out time integration iterative solution, and updating to obtain the grid parameters of the current time step.
And S4, updating the dynamic calculation domain based on the residual error and the grid parameters obtained in the S3.
Based on the grid parameters of the current time step, judging whether the convection dynamic domain needs to be increased or not according to the speed and the sound velocity after the preprocessing, and then judging whether the convection dynamic domain needs to be decreased or not by combining the convergence condition of the residual error;
and judging whether the current grid unit is dominated by a viscosity effect or not, if so, adding all the adjacent grid units of the current grid unit to a viscosity dynamic domain, and if not, deleting the current grid unit from the viscosity dynamic domain.
The increasing and decreasing of the convection dynamic domain in S4 specifically includes:
judging whether the current grid cell is disturbed without adhesion by measuring the grid parameter updating amount obtained by calculation in S3, if so, determining the adjacent grid cells which can be influenced, and adding all the adjacent grid cells into a convection dynamic domain, wherein for the cells in subsonic flow, all the adjacent grid cells can be influenced by the disturbance; for the unit in supersonic flow, determining the propagation direction of disturbance according to the characteristic value of convection flux, and adding the grid unit at the downstream of the propagation direction into a convection dynamic domain;
and judging whether the current grid unit is converged, if so, judging whether the current grid unit is positioned at the most upstream of the convection dynamic domain, if so, judging whether other grid units in the convection dynamic domain have no influence on the current grid unit, and if not, deleting the current grid unit from the convection dynamic domain.
The preprocessed speed and sound velocity are obtained based on the local speed and sound velocity, and specifically include:
obtaining a preprocessed speed and a preprocessed sound velocity based on the local speed u and the sound velocity c;
obtaining a preprocessed characteristic value based on the preprocessed speed and the sound velocity,
the calculation formulas of the speed and the sound velocity after the pretreatment are respectively as follows:
λ＝[u,u,u,u′+c′,u′-c′]
U r 2 ＝min(c 2 ,max(|u| 2 ,K|u| 2 ))
wherein, λ is a characteristic value used in calculating the convection flux, U is a local speed, U 'is a preprocessed speed, c is a sound velocity, c' is a preprocessed sound velocity, α and β are parameters used for solving the preprocessed characteristic value, and U is a local speed r For the pre-processing of the parameters, given by the pre-processing method, ρ p ,ρ T Respectively density, density derivative to pressure, density derivative to temperature, calculated from flow field parameters, C p K is a small value (10) for the specific heat constant at constant pressure -5 )。
The increasing and decreasing of the viscous dynamic domain in S4 specifically includes:
dividing the threshold of the total enthalpy change amount by the Mach number square as a viscosity judgment basis, and increasing and reducing a viscosity dynamic domain, wherein the specific method is as follows:
if |1-h 0,∞ /h 0 |/Ma 2 ≥ε d,v Then all immediately adjacent grid cells of the current grid cell are added toA viscous dynamic domain;
if |1-h 0,∞ /h 0 |/Ma 2 <ε d,v If yes, deleting the current grid cell from the sticky dynamic domain;
wherein h is 0,∞ For the total enthalpy of the incoming flow, h 0 The total enthalpy of the current grid unit, Ma is the Mach number of the current grid unit, epsilon d,v A threshold is determined for a given viscosity domain.
The grid parameters include:
density, velocity, pressure and temperature.
Illustratively, in the present embodiment, the sticky domain determination threshold is 0.001, which is suitable for the requirement of the perturbation domain update method for compressible flow calculation. For this embodiment, the calculation result of determining the sticky partition by dividing the mach number by the square is compared. Taking a grid cell close to a wall as an example, the grid must be in the viscosity-dominant region. However, the total enthalpy of the incoming flow obtained by non-dimensionalization calculation is 25000.5, the total enthalpy of the current grid unit is 25002.0064556779, and the original judgment is based on |1-h 0,∞ /h 0 The result of the | calculation was 6.025 × 10 -5 Much less than a given viscosity decision threshold. After dividing by the Mach number square, the grid cell is determined according to |1-h 0,∞ /h 0 |/Ma 2 0.6025, the grid unit can be judged to be in the viscosity dominant region, which proves that the technical means can effectively judge the viscosity subarea of the non-compressible flow, and in the compressible flow with higher speed, the square of the Mach number is close to 1, and the original judgment basis is not changed for the applicability of the compressible flow, namely the technical means adopted by the invention can be applied to the calculation of the compressible/non-compressible mixed flow.
In general, the viscosity-dominated region of this embodiment remains at 1.4% (i.e., the initial viscosity domain) until the division by the Mach number square is employed, losing the correctness of the calculation. After the Mach number square is divided by the viscosity-dominated region, the occupation ratio of the viscosity-dominated region is quite obvious as shown in FIG. 2, the viscosity-dominated region is more in line with physical laws, and accurate calculation results can be obtained.
The reduction of the sticky dynamic field is determined and performed. The reduction of the sticky dynamic field corresponds to the expansion of the sticky dynamic field in the previous step, and when the sticky effect does not occupy the dominant factor, the grid unit can be deleted from the sticky dynamic field, namely the range of the sticky dynamic field is reduced.
S5, judging whether the global residual is less than the given convergence residual threshold, if yes, finishing the calculation and executing S6, and if not, returning to S3.
And judging whether the calculation is finished or not according to the size relation between the global residual and the given convergence residual threshold, namely whether the global residual is smaller than the given convergence residual threshold or not, if so, performing S6, and otherwise, skipping to S3.
(1) If the global residual is smaller than the given threshold, the calculation is completed, and the process proceeds to S6.
R<ε lim
Wherein R is the global residual, ε lim Given a convergence residual threshold.
(2) When the global residual is equal to or greater than the given convergence residual threshold, the process proceeds to S3.
R≥ε lim
And S6, outputting the calculated network parameter value.
After the calculation is completed, the numerical values of the physical quantities obtained by the calculation are output.
Example 1
Fig. 1 is a schematic diagram of an efficient disturbance domain propulsion method for a compressible/non-compressible mixed flow according to the present invention, which is described in detail below. The invention provides a disturbance domain updating method suitable for compressible/non-compressible flow, which takes the solution of a two-dimensional airfoil NACA0012 flow field as an example, and the implementation mode of applying a preprocessing method to the solution of a compressible Navier-Stokes equation comprises the following steps:
the first step is as follows: data reading and flow field initialization.
(1) Reading in specific calculation parameters of the NACA0012 airfoil profile, wherein the specific parameters of the embodiment are as follows: mach number 0.01, Reynolds number 10000, laminar viscous flow;
(2) reading the calculation grid of the NACA0012 airfoil profile and setting of boundary conditions;
(3) in the embodiment, a method of initializing according to the incoming flow conditions is adopted, and the flow field parameters on all the computational grids are initialized by adopting the incoming flow boundary conditions.
The second step is that: an initial dynamic computing domain is established.
(1) Because a method of initializing according to the incoming flow condition is adopted, the convection and viscous dynamic domains are established only according to the wall boundary;
the third step: and iteratively solving the residual estimation and the time integration.
(1) And according to different calculation conditions, selecting a proper reconstruction format and a proper convection format, and solving the residual error of the compressible Navier-Stokes equation at the current time step. In the embodiment, a preprocessed Roe format is selected, the convection flux is calculated, and the residual error of the current time step is obtained; characteristic values are used when the convection flux is calculated, and the characteristic values comprise local speed and sound velocity.
(2) And after obtaining the residual, carrying out time advance by adopting an explicit or implicit time advance solving method, and updating to obtain the physical quantity of the current time step calculation domain. In this embodiment, an implicit LU-SGS format is selected to perform time-marching solution, and the physical quantity update quantity of the current time step is obtained.
The fourth step: dynamic domain update
(1) And judging and executing the increase of the drainage basin. Firstly, whether the current grid unit is disturbed without adhesion is judged by measuring the updated physical quantity calculated in the third step. If inviscid, the immediate vicinity of the cell that will be affected needs to be determined. If the unit is in subsonic flow, all the adjacent units are affected by disturbance, if the unit is in supersonic flow, the propagation direction of the disturbance is determined according to the characteristic value of convection flux, so that the disturbed unit is determined, and the grid units at the downstream of the propagation direction are added into a convection dynamic domain.
The characteristic value used in the calculation of the convection flux needs to adopt a preprocessed characteristic value, and the preprocessed speed and the speed of sound adopted in the calculation of the characteristic value are obtained on the basis of the local speed and the speed of sound.
In this step, the eigenvalue used in calculating the convection flux needs to adopt the preprocessed eigenvalue, and the solving method is as follows:
λ＝[u,u,u,u′+c′,u′-c′]
the calculation of the preprocessed speed and sound velocity is as follows:
wherein, λ is a characteristic value used in calculating the convection flux, U is a local velocity, U 'is a velocity after preprocessing, c is a sound velocity, c' is a preprocessing sound velocity, U is a local velocity, and r for the pre-processing of the parameters, given by the pre-processing method, ρ p ,ρ T Respectively density, density derivative to pressure, density derivative to temperature, calculated from flow field parameters, C p Is a constant specific heat capacity constant.
(2) And judging and executing the reduction of the basin. Firstly, judging whether the current unit is converged, if so, judging whether the current grid unit is positioned at the most upstream of the convection dynamic calculation domain, if so, judging whether other units in the convection dynamic domain have no influence on the unit, and deleting the current grid unit from the convection dynamic domain.
If other units no longer affect the unit, the watershed can be reduced.
(3) The increase of the sticky domain is judged and performed. And judging whether the unit is dominated by the viscous effect, and if the unit is dominated by the viscous unit, adding all the adjacent units to the viscous domain. The invention provides a viscosity leading judgment basis suitable for solving a compressible/non-compressible flow field by preprocessing:
|1-h 0,∞ /h 0 |×Ma 2 ≥ε d,v
wherein h is 0,∞ For the total enthalpy of the incoming flow, h 0 Is the total enthalpy of the current cell, Ma is the Mach number of the current cell, epsilon d,v For a given viscosity domain decision threshold, 0.001 is taken in this example, which is a value that is applicable to the requirements of the perturbation domain update method for compressible flow calculations.
For this embodiment, the sticky partition determination is made by comparing whether to divide by Mach number squaredAnd calculating a result of the break. Taking a grid cell close to a wall as an example, the grid must be in the viscosity-dominant region. However, the total enthalpy of the incoming flow obtained by non-dimensionalization calculation is 25000.5, the total enthalpy of the current unit is 25002.0064556779, and the original judgment is based on |1-h 0,∞ /h 0 The result of the | calculation was 6.025 × 10 -5 Much less than a given viscosity decision threshold. After dividing by the Mach number square, the unit determines the basis |1-h 0,∞ /h 0 |×Ma 2 0.6025, the cell is in the viscosity dominant region, which proves that the technical means can effectively judge the viscosity subarea of the non-compressible flow, while in the compressible flow with higher speed, the square of the Mach number is close to 1, and the original judgment basis is not changed, namely the technical means can be applied to the calculation of the compressible flow.
FIG. 2 is a schematic illustration of the ratio of convective dynamic field to viscous-dominated region provided by the present invention.
In general, if the Mach number square is not divided, the viscosity-dominated region of the present embodiment is always kept at 1.4% (i.e., the initial viscosity domain), and the calculation is lost. After the Mach number square is divided, the occupation ratio of the viscous dominant region is as shown in FIG. 2, which obviously conforms to the physical law better and can obtain a correct calculation result.
As shown in fig. 2, size _ a and size _ v refer to the proportion of the convection dynamic domain and the viscous dynamic domain, respectively, in all the calculation domains.
(4) The reduction of the sticky domain is determined and performed. The narrowing of the sticky field corresponds to the widening of the sticky field in the previous step, and when the sticky effect does not dominate the effect, the unit can be deleted from the convection field.
The fifth step: and judging whether the calculation is finished.
(1) And judging the size relation between the global residual and a given convergence residual threshold, if the global residual is smaller than the given threshold, finishing the calculation, and entering the sixth step.
R<ε lim
Wherein R is the global residual, ε lim For a given convergence residual threshold, i.e. if the global residual is less thanAnd (5) determining a threshold value, finishing the calculation and entering the sixth step.
(2) Otherwise, if the global residual is still larger than the given threshold, the process loops to the third step.
And a sixth step: and outputting a calculation result.
After the calculation is completed, the numerical values of the physical quantities obtained by the calculation are output.
FIG. 3a is a schematic diagram illustrating a comparison of calculation results of a disturbance domain updating method without using a compressible/non-compressible flow according to an embodiment of the present invention; fig. 3b is a schematic diagram comparing calculation results of a perturbation domain updating method using compressible/non-compressible flow according to an embodiment of the present invention.
Fig. 4 is a schematic diagram comparing convergence situations of a perturbation domain updating method for determining whether to use compressible/non-compressible flow according to an embodiment of the present invention.
As shown in fig. 4, DRUM refers to the method of the present invention, GUM refers to the case where the method of the present invention is not employed, the abscissa is the number of iteration steps, and the ordinate is the maximum residual convergence.
After the calculation is completed in this embodiment, for example, as shown in fig. 3a, fig. 3b and fig. 4, it can be seen whether the calculation result and the calculation convergence condition of the disturbance domain updating method of compressible/non-compressible flow provided by the present invention are adopted, that is, the pressure distribution on the airfoil surface is completely consistent, and the total residual convergence condition is basically consistent, that is, the method provided by the present invention is adopted, so that the calculation accuracy is ensured. After the calculation method provided by the invention is applied, the total calculation time of the embodiment is reduced from 40183.3s to 26197.6s, and the time is saved by 35%, so that the calculation efficiency of the compressible/non-compressible mixed flow can be greatly improved.
According to the technical scheme provided by the embodiment of the disclosure, a disturbance propagation calculation formula of the compressible/incompressible disturbance domain updating method is established based on the preprocessed local sound velocity and the preprocessed local speed, so that a criterion of the compressible/incompressible disturbance domain propagation and influence is established, the problems of change of the characteristic value of the Jacobian matrix and inapplicability of the criterion of the low-Mach number area are solved, and the calculation domain is contracted; the viscous leading judgment basis suitable for solving the compressible/non-compressible flow field through pretreatment is provided, the judgment threshold value of the viscous subarea is obtained by dividing the threshold value of the total enthalpy variable quantity by the Mach number square, the dependency of the viscous subarea on the Mach number when the total enthalpy variable quantity is used as the criterion is avoided, and the calculation efficiency is greatly improved.
All the above optional technical solutions may be combined arbitrarily to form optional embodiments of the present application, and are not described herein again.
It should be understood that, the sequence numbers of the steps in the foregoing embodiments do not imply an execution sequence, and the execution sequence of each process should be determined by its function and inherent logic, and should not constitute any limitation on the implementation process of the embodiments of the present disclosure.
The above examples are only intended to illustrate the technical solutions of the present disclosure, not to limit them; although the present disclosure has been described in detail with reference to the foregoing embodiments, it should be understood by those of ordinary skill in the art that: the technical solutions described in the foregoing embodiments may still be modified, or some technical features may be equivalently replaced; such modifications and substitutions do not substantially depart from the spirit and scope of the embodiments of the present disclosure, and are intended to be included within the scope of the present disclosure.
Claims (9)
1. A high-efficiency disturbance domain propulsion method for compressible/non-compressible mixed flow is characterized by comprising the following steps:
s1, reading data and initializing a flow field;
s2, establishing an initial dynamic calculation domain according to the flow field initialization mode, wherein the dynamic calculation domain comprises a convection dynamic domain and a viscosity dynamic domain;
s3, in the dynamic calculation domain, solving the residual error of the compressible Navier-Stokes equation at the current time step, and carrying out time integration iteration to solve the grid parameters;
s4, updating the dynamic calculation domain based on the residual error and the grid parameters obtained in S3;
s5, judging whether the global residual is smaller than a given convergence residual threshold, if so, finishing the calculation and executing S6, and if not, returning to S3;
s6 outputs the calculated values of the mesh parameters.
2. The efficient perturbation domain updating method according to claim 1, wherein the step S1 comprises:
reading in computational grids, boundary conditions, grid parameters and computational settings;
based on the grid parameters, initializing the flow field according to the incoming flow conditions or the given flow field, and determining initial values of physical quantities of which boundaries and domains are not defined before calculation and formula parameters in the selected physical model, wherein the formula parameters in the physical model comprise a gas constant, a viscosity coefficient and a turbulence model coefficient.
3. The efficient perturbation domain updating method according to claim 2, wherein the step S2 comprises:
when the dynamic viscous dynamic domain is initialized according to the incoming flow conditions, all grid parameters are set as incoming flow parameters, the convective dynamic domain is established according to the wall surface conditions, and only 1 layer of grid units adjacent to the wall surface is taken as an initial viscous dynamic domain;
when the method is initialized according to a given flow field, all grid parameters are set according to the given flow field, the viscous dynamic domain is established according to the difference between the given flow field and an incoming flow, grid units which are not consistent with the flowing characteristics of the incoming flow are used as an initial convection dynamic domain, and a viscous leading grid unit in the convection dynamic domain is selected as a viscous dynamic domain.
4. The efficient perturbation domain updating method according to claim 1, wherein the step S3 comprises:
initializing all given grid parameters based on a flow field, solving a compressible Navier-Stokes equation, selecting a preprocessed Roe format, calculating convection flux, and obtaining a residual error of the current time step;
and based on the residual error result, carrying out time integral iterative solution, and updating to obtain the grid parameters of the current time step.
5. The efficient perturbation domain updating method according to claim 4, wherein the S4 comprises:
based on the grid parameters of the current time step, judging whether the convection dynamic domain needs to be increased or not according to the speed and the sound velocity after the preprocessing, and then judging whether the convection dynamic domain needs to be decreased or not by combining the convergence condition of the residual error;
and judging whether the current grid unit is dominated by a viscosity effect or not, if so, adding all the adjacent grid units of the current grid unit to a viscosity dynamic domain, and if not, deleting the current grid unit from the viscosity dynamic domain.
6. The efficient disturbance domain updating method according to claim 5, wherein the increasing and decreasing of the convection dynamic domain in S4 specifically includes:
judging whether the current grid unit is disturbed without adhesion by measuring the grid parameter updating amount calculated in S3, if so, determining the adjacent grid unit which can be influenced, and adding all the adjacent grid units into the convection dynamic domain, wherein for the unit in subsonic flow, all the adjacent grid units can be influenced by the disturbance; for the unit in supersonic flow, determining the propagation direction of disturbance according to the characteristic value of convection flux, and adding the grid unit at the downstream of the propagation direction into a convection dynamic domain;
and judging whether the current grid unit is converged, if so, judging whether the current grid unit is positioned at the most upstream of the convection dynamic domain, if not, judging whether other grid units in the convection dynamic domain no longer influence the current grid unit, and if not, deleting the current grid unit from the convection dynamic domain.
7. The efficient disturbance domain updating method according to claim 6, wherein the preprocessed speed and sound speed are obtained based on local speed and sound speed, and specifically comprises:
obtaining a preprocessed speed and a preprocessed sound speed based on the local speed u and the sound speed c;
obtaining a preprocessed characteristic value based on the preprocessed speed and the sound velocity,
the calculation formulas of the speed and the sound velocity after the pretreatment are respectively as follows:
λ＝[u,u,u,u′+c′,u′-c′]
u′＝u(1-α),
α＝(1-βU r 2 )/2,
U r 2 ＝min(c 2 ,max(|u| 2 ,K|u| 2 ))
wherein, λ is a characteristic value used in calculating the convection flux, U is a local speed, U 'is a preprocessed speed, c is a sound velocity, c' is a preprocessed sound velocity, α and β are parameters used for solving the preprocessed characteristic value, and U is a local speed r For the pre-processing of the parameters, given by the pre-processing method, ρ p ,ρ T Respectively density, density derivative to pressure, density derivative to temperature, calculated from flow field parameters, C p K is a small value (10) for the specific heat constant at constant pressure -5 )。
8. The efficient perturbation domain updating method according to claim 5, wherein the increasing and decreasing of the sticky dynamic domain in S4 specifically comprises:
dividing the threshold of the total enthalpy change amount by the Mach number square as a viscosity judgment basis, and increasing and reducing the viscosity dynamic domain, specifically as follows:
if |1-h 0,∞ /h 0 |/Ma 2 ≥ε d,v Adding all the immediate neighbor grid cells of the current grid cell to the sticky dynamic field;
if |1-h 0,∞ /h 0 |/Ma 2 <ε d,v If yes, deleting the current grid cell from the sticky dynamic domain;
wherein h is 0,∞ For the total enthalpy of the incoming flow, h 0 The total enthalpy of the current grid unit, Ma is the Mach number of the current grid unit, epsilon d,v A threshold is determined for a given viscosity domain.
9. The efficient perturbation domain updating method according to claim 6, wherein the grid parameters comprise: density, velocity, pressure and temperature.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210469835.6A CN115017604B (en) | 2022-04-28 | 2022-04-28 | Efficient disturbance domain propulsion method for compressible/non-compressible mixed flow |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210469835.6A CN115017604B (en) | 2022-04-28 | 2022-04-28 | Efficient disturbance domain propulsion method for compressible/non-compressible mixed flow |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115017604A true CN115017604A (en) | 2022-09-06 |
| CN115017604B CN115017604B (en) | 2024-06-07 |
Family
ID=83067905
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210469835.6A Active CN115017604B (en) | 2022-04-28 | 2022-04-28 | Efficient disturbance domain propulsion method for compressible/non-compressible mixed flow |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115017604B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117787147A (en) * | 2024-02-28 | 2024-03-29 | 中国空气动力研究与发展中心计算空气动力研究所 | Model building method for solving full-speed domain flow problem and related products |
Citations (4)
| 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 |
| CN111651831A (en) * | 2020-04-13 | 2020-09-11 | 北京航空航天大学 | Partition disturbance domain updating calculation method for constant-viscosity compressible streaming of aircraft |
| CN113392472A (en) * | 2021-08-17 | 2021-09-14 | 北京航空航天大学 | OpenMP parallel disturbance domain updating method for aircraft aerodynamic characteristic simulation |
| CN113609598A (en) * | 2021-10-09 | 2021-11-05 | 北京航空航天大学 | RANS/LES disturbance domain updating method for aircraft aerodynamic characteristic simulation |
-
2022
- 2022-04-28 CN CN202210469835.6A patent/CN115017604B/en active Active
Patent Citations (4)
| 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 |
| CN111651831A (en) * | 2020-04-13 | 2020-09-11 | 北京航空航天大学 | Partition disturbance domain updating calculation method for constant-viscosity compressible streaming of aircraft |
| CN113392472A (en) * | 2021-08-17 | 2021-09-14 | 北京航空航天大学 | OpenMP parallel disturbance domain updating method for aircraft aerodynamic characteristic simulation |
| CN113609598A (en) * | 2021-10-09 | 2021-11-05 | 北京航空航天大学 | RANS/LES disturbance domain updating method for aircraft aerodynamic characteristic simulation |
Non-Patent Citations (2)
| Title |
|---|
| M. KRAPOSHIN 等: "Numerical algorithm based on regularized equations for incompressible flow modeling and its implementation in OpenFOAM", 《COMPUTER SCIENCE》, 1 October 2021 (2021-10-01) * |
| 刘中玉;张明锋;郑冠男;杨国伟;: "基于预处理HLLEW格式的全速域数值算法", 计算物理, no. 03, 25 May 2016 (2016-05-25), pages 25 - 34 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117787147A (en) * | 2024-02-28 | 2024-03-29 | 中国空气动力研究与发展中心计算空气动力研究所 | Model building method for solving full-speed domain flow problem and related products |
| CN117787147B (en) * | 2024-02-28 | 2024-05-03 | 中国空气动力研究与发展中心计算空气动力研究所 | Model building method for solving full-speed domain flow problem and related products |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115017604B (en) | 2024-06-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bijl et al. | Implicit time integration schemes for the unsteady compressible Navier–Stokes equations: laminar flow | |
| CN105718634B (en) | A kind of aerofoil profile Robust Optimal Design based on non-probability interval analysis model | |
| CN111859530B (en) | Iterative propulsion disturbance domain updating method for aircraft dynamic aerodynamic characteristic simulation | |
| Biancolini et al. | Static aeroelastic analysis of an aircraft wind-tunnel model by means of modal RBF mesh updating | |
| Sørensen et al. | Computation of wind turbine wakes using combined Navier-Stokes actuator-line Methodology | |
| Verstraete et al. | Stability analysis of partitioned methods for predicting conjugate heat transfer | |
| Imamura et al. | Acceleration of steady-state lattice Boltzmann simulations on non-uniform mesh using local time step method | |
| CN111079228B (en) | Pneumatic shape optimization method based on flow field prediction | |
| CN112016167A (en) | Aircraft aerodynamic shape design method and system based on simulation and optimization coupling | |
| CN111651831B (en) | Partition disturbance domain updating calculation method for constant-viscosity compressible streaming of aircraft | |
| CN111859529B (en) | Multi-grid disturbance domain updating acceleration method for aircraft streaming numerical simulation | |
| CN108563843A (en) | The disturbance region update method of steady compressible flowing | |
| CN115329689A (en) | High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion | |
| CN115017604A (en) | Efficient disturbance domain propulsion method for compressible/incompressible mixed flow | |
| Peng et al. | Nested Cartesian grid method in incompressible viscous fluid flow | |
| Zhang et al. | Research on data assimilation strategy of turbulent separated flow over airfoil | |
| Meng et al. | Fast flow prediction of airfoil dynamic stall based on Fourier neural operator | |
| CN116628854A (en) | Wing section aerodynamic characteristic prediction method, system, electronic equipment and storage medium | |
| Kerestes et al. | LES Modeling of High-Lift High-Work LPT Blades: Part II—Validation and Application | |
| CN115470653A (en) | Rigid chemical reaction flow semi-hidden semi-obvious space-time multiple time step propulsion method | |
| Buccafusca et al. | Maximizing power in wind turbine arrays with variable wind dynamics | |
| De Boer et al. | Radial basis functions for interface interpolation and mesh deformation | |
| CN112632825B (en) | Electrostatic field smooth finite element numerical algorithm based on finite element super-convergence | |
| CN117077298B (en) | Aircraft robust optimization design method based on gradient enhancement random Co-Kriging model | |
| CN114117950B (en) | Flutter judgment method for shuttle aircraft based on acting principle |
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 |