CN116244853A - Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine - Google Patents

Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine Download PDF

Info

Publication number
CN116244853A
CN116244853A CN202310145756.4A CN202310145756A CN116244853A CN 116244853 A CN116244853 A CN 116244853A CN 202310145756 A CN202310145756 A CN 202310145756A CN 116244853 A CN116244853 A CN 116244853A
Authority
CN
China
Prior art keywords
time step
calculation
diagonal
dominant
convergence
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
CN202310145756.4A
Other languages
Chinese (zh)
Other versions
CN116244853B (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 CN202310145756.4A priority Critical patent/CN116244853B/en
Publication of CN116244853A publication Critical patent/CN116244853A/en
Application granted granted Critical
Publication of CN116244853B publication Critical patent/CN116244853B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/12—Geometric CAD characterised by design entry means specially adapted for CAD, e.g. graphical user interfaces [GUI] specially adapted for CAD
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Fluid Mechanics (AREA)
  • Architecture (AREA)
  • Human Computer Interaction (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a harmonic balance method accelerating convergence method for predicting unsteady flow of a turbine, which comprises the following main implementation processes: step one, performing space and time dispersion on a control equation; step two, calculating a local time step; step three, calculating parameters used by a time pushing format; step four, calculating a main diagonal dominant threshold d; step five, solving a new time step Dii according to the dominant threshold d of the main diagonal; step six, the new time step Dii is carried into an equation, and flow field calculation is carried out. According to the invention, the local time step is optimized based on the local flow field characteristics, so that the rapid convergence of the harmonic balancing method is realized while the calculation accuracy is ensured, the numerical stability is ensured, the convergence speed is also increased, the calculation consumption is effectively saved, and the problem of overlong calculation period in engineering application is solved.

Description

Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine
Technical Field
The invention relates to the field of engineering computational fluid dynamics, in particular to a harmonic balance method accelerating convergence method for predicting unsteady flow of a turbine.
Background
There are two main approaches to hydrodynamic studies, experimental and computational fluid dynamics (Computational Fluid Dynamics, CFD for short), respectively; the experimental measurement is affected by the arrangement of the measuring points, and global data of the flow field cannot be obtained; along with the rapid development of computer computing power in recent years, CFD technology gradually becomes an important technical means for theoretical research and engineering application in the fields of aerodynamics, hydrodynamics and the like.
In complex flows of the turbine stage, there is an inherent static-rotating interference unsteady flow, which has an important influence on turbine performance and needs to be accurately predicted in the design process. For a multistage turbine, the conventional steady mixing surface simulation method ignores the steady flow of static-rotating interference in the turbine, and the simulation deviation is increased along with the increase of the turbine stage number; the traditional unsteady simulation method is too high in calculation cost and difficult to apply to daily engineering design. The harmonic balancing method is a method for converting a balanced simulated turbine into static interference unsteady flow, the method utilizes a steady simulation method to calculate flow field data at a plurality of moments, then fits flow fields at different moments through Fourier series, and then reconstructs a complete unsteady flow field; compared with the traditional unsteady method, the harmonic balancing method can be better in balancing precision and calculated amount, so that the method has higher engineering application value.
Although the harmonic balance method saves computational resources compared with the unsteady simulation method, and the convergence speed is much faster than that of the unsteady simulation method; however, the harmonic balance method still has problems of poor convergence and slow convergence speed compared with the steady simulation method. In order to ensure convergence, the existing harmonic balancing method adopts a mode of reducing the number of CFLs (Courant-Friedrichs-Lewy numbers) and sacrificing the calculation speed to ensure the convergence. In order to further accelerate the convergence speed of the harmonic balancing method so as to save the computing resources, development of a method for accelerating the convergence of the harmonic balancing method is needed.
Disclosure of Invention
First, the technical problem to be solved
The invention aims to provide a harmonic balance method accelerating convergence method for predicting unsteady flow of a turbine, which aims to realize rapid convergence of the harmonic balance method by optimizing local time step based on local flow field characteristics, thereby saving calculation consumption and solving the problem of overlong calculation period in engineering application.
(II) technical scheme
In order to solve the technical problems, the invention provides a harmonic balance method accelerating convergence method for predicting unsteady flow of a turbine, which comprises the following steps:
step one, performing space and time dispersion on a control equation;
step two, calculating a local time step;
step three, calculating parameters used by a time pushing format;
step four, calculating a main diagonal dominant threshold d;
step five, solving a new time step D according to the dominant threshold D of the main diagonal ii ；
Step six, a new time step D ii Carrying out flow field calculation in the equation;
(1) the spatially and temporally discretizing the control equation includes:
the space flux is calculated by adopting a space flux calculation method of a traditional steady Reynolds average Navier-Stokes (RANS) method, the time pushing format adopts a BJ-SOR (Block Jacobi Successive Over Relaxation) method, and further, in order to ensure the convergence of calculation, the source item is subjected to implicit discrete, and the calculation method is given by the following formula:
Figure BDA0004089110530000011
AΔW＝-R+C(ΔW)
(L+D+U)ΔW＝-R+C(ΔW)
wherein V represents grid volume, deltat represents physical time step, and I is an identity matrix; r is a residual term, W is a conservation quantity, deltaW is a variation quantity of the conservation quantity, and C (DeltaW) is the dispersion of a harmonic source term; a is a jacobian matrix corresponding to the whole linear system, wherein the jacobian matrix A corresponding to the whole linear system is decomposed into a lower diagonal matrix L, a diagonal matrix D and an upper diagonal matrix U; the time discrete adopts a first-order Euler post-difference format;
the control equation adopted by the BJ-SOR method is given by:
(D+L)ΔW n+1 / 2 ＝-[R+C(ΔW)]-UΔW n
(D+U)ΔW n+1 ＝-[R+C(ΔW)]-LΔW n+1 / 2
solving a control equation adopted by the BJ-SOR method by adopting multiple cycles, and adopting 4-6 times of calculation in the calculation of the harmonic balancing method to achieve better convergence;
(2) the calculating the local time step includes:
calculating an original time step and a time step corrected by a harmonic balance method:
wherein the calculation of the original time step Δt is given by:
Figure BDA0004089110530000021
the time step delta t corrected by the harmonic balancing method HB The calculation of (2) is given by:
Figure BDA0004089110530000022
wherein lambda is the spectral radius of the jacobian matrix A corresponding to the whole linear system, omega is the angular frequency adopted in the harmonic balancing method, N is the order of the harmonic balancing method, the value of the CFL number is the same as that of the original algorithm, and omega is the grid volume;
(3) the parameters used in the calculation time advance format include:
calculating a diagonal matrix L, a diagonal matrix D and an upper diagonal matrix U which are generated after the BJ-SOR method adopted by the time pushing format is discretized; the diagonal matrix D is based on the original timeDiagonal matrix D (Δt) of inter-step Δt and time-step Δt corrected based on said harmonic balancing method HB Diagonal matrix D (Δt) HB ) The calculation method is given by the following formula:
Figure BDA0004089110530000023
where ΔS is the area of a single mesh plane, Λ c Spectral radius for non-sticking flux; the superscript I, J, K indicates different coordinate directions of the structural grid, the subscripts i, j, k indicate index values of grid coordinates, the subscript c indicates non-sticky flux, and the subscript v indicates sticky flux;
(4) the calculating the dominant diagonal threshold d includes:
the maximum time step when the main diagonal is dominant is obtained by selecting the main diagonal dominant threshold d of the jacobian matrix A corresponding to the whole linear system; the dominant diagonal threshold d is given by:
Figure BDA0004089110530000024
wherein C is d Taking a value between 0.68 and 0.8 as the dominant diagonal coefficient, and properly reducing the dominant diagonal coefficient C d The value of (2) can increase the convergence rate;
the maximum time step when the main diagonal is dominant is obtained and used for calculation, so that the problem that the convergence is poor due to the fact that the main diagonal dominant characteristic of the jacobian matrix A corresponding to the whole linear system cannot be ensured by the time step deltat adopted in the original harmonic balance method can be effectively solved, and the convergence is effectively accelerated;
(5) the new time step D is solved according to the dominant threshold D of the main diagonal ii Comprising the following steps:
solving for a new time step D using ii ；
Figure BDA0004089110530000025
(6) Said new time step D ii In the carry-over equation, performing flow field calculations includes:
the new time step D solved in (5) ii Carrying out flow field calculation in the calculation method in the step (1);
in particular, the new time step D ii According to different values of the main diagonal dominant threshold d, the time step delta t after the correction of the original time step delta t and the harmonic balance method is carried out HB The switching is carried out, so that most grids adopt the original time step delta t in the calculation process, and the small grids adopt the time step delta t corrected by the harmonic balancing method HB Thereby realizing the effect of accelerating convergence.
(III) beneficial effects
The harmonic balance method acceleration convergence method for predicting the unsteady flow of the turbine provided by the invention has the following beneficial effects: the maximum time step of the dominant diagonal of the jacobian matrix A corresponding to the whole linear system in the numerical calculation process is obtained by selecting the dominant diagonal threshold d, and the maximum time step is used for calculation, so that the problem that the dominant diagonal of the jacobian matrix A corresponding to the whole linear system cannot be guaranteed by the time step deltat adopted in the original harmonic balancing method, and therefore the convergence is poor is effectively accelerated;
the method does not involve parallel processing, but only involves local variables, so that the code is simpler and is easy to expand into the existing program;
the method has the advantages that the method shows better application effect in the test of a plurality of examples, the convergence speed is improved to a certain extent, and the effect is improved by more than 35% compared with the convergence speed of the obvious examples; the main diagonal dominant coefficient C related by the method of the invention d Can be adjusted according to specific examples in the use process, and properly reduce C while ensuring stability d The value of (2) can achieve the effect of accelerating the convergence speed;
the method can realize a larger gain effect under a smaller CFL number, thereby ensuring the numerical stability and accelerating the convergence rate.
Drawings
FIG. 1 is a flow chart of a harmonic balancing method accelerating convergence method for predicting unsteady flow of a turbine;
FIG. 2 is a grid schematic diagram of an embodiment Stage35 compressor of the present invention at 50% leaf height;
fig. 3 is a graph showing the comparison of the convergence speed of two-dimensional flow sheets at 50% of the height of the blade calculated by using different methods in the Stage35 compressor according to the embodiment of the invention.
Detailed Description
The following describes in further detail the embodiments of the present invention with reference to the drawings and examples. The following examples are only illustrative of the present invention and are not intended to limit the scope of the invention.
The invention provides a harmonic balance method accelerating convergence method for predicting unsteady flow of a turbine, which comprises the following steps:
step one, performing space and time dispersion on a control equation;
in this step, a space flux calculation method of a traditional steady Reynolds average Navier-Stokes (RANS) method is adopted to calculate the space flux, a BJ-SOR (Block Jacobi Successive Over Relaxation) method is adopted in a time pushing format, further, in order to ensure the convergence of calculation, the source item is subjected to implicit discrete, and the calculation method is given by the following formula:
Figure BDA0004089110530000031
AΔW＝-R+C(ΔW)
(L+D+U)ΔW＝-R+C(ΔW)
wherein V represents grid volume, deltat represents physical time step, and I is an identity matrix; r is a residual term, W is a conservation quantity, deltaW is a variation quantity of the conservation quantity, and C (DeltaW) is the dispersion of a harmonic source term; a is a jacobian matrix corresponding to the whole linear system. Decomposing a jacobian matrix A corresponding to the whole linear system into a lower diagonal matrix L, a diagonal matrix D and an upper diagonal matrix U; the time dispersion adopts a first-order Euler post-difference format.
The control equation for the BJ-SOR method employed by the time-marching format is given by:
(D+L)ΔW n+1 / 2 ＝-[R+C(ΔW)]
(D+U)ΔW n+1 ＝-[R+C(ΔW)]-LΔW n+1 / 2
the control equation adopted by the BJ-SOR method is solved by adopting multiple cycles, and 4-6 times of calculation are adopted in the calculation of the harmonic balancing method, so that better convergence can be achieved.
Step two, calculating a local time step;
in this step, the calculation of the original time step Δt is given by:
Figure BDA0004089110530000032
time step delta t corrected by harmonic balancing method HB The calculation of (2) is given by:
Figure BDA0004089110530000033
wherein lambda is the spectral radius of the jacobian matrix A corresponding to the whole linear system, omega is the angular frequency adopted in the harmonic balancing method, N is the order of the harmonic balancing method, the value of the CFL number is the same as that of the original algorithm, and omega is the grid volume;
by the expression, the original time step delta t and the time step delta t corrected by the harmonic balancing method are respectively calculated HB 。
Step three, calculating parameters used by a time pushing format;
in the step, calculating a diagonal matrix L, a diagonal matrix D and an upper diagonal matrix U which are generated after the BJ-SOR method adopted by the time pushing format is discretized; the diagonal matrix D comprises a diagonal matrix D (delta t) based on an original time step delta t and a time step delta t corrected based on a harmonic balance method HB Diagonal matrix D (Δt) HB ) The calculation method is given by the following formulaAnd (3) out:
Figure BDA0004089110530000041
where ΔS is the area of a single mesh plane, Λ c Spectral radius for non-sticking flux; the superscript I, J, K indicates different coordinate directions of the structural grid, the subscripts i, j, k indicate index values of grid coordinates, the subscript c indicates non-sticky flux, and the subscript v indicates sticky flux;
step four, calculating a main diagonal dominant threshold d;
in the step, the maximum time step when the main diagonal is dominant is obtained by selecting the main diagonal dominant threshold d of the jacobian matrix A corresponding to the whole linear system; wherein the dominant diagonal threshold d is given by:
Figure BDA0004089110530000042
C d taking a value between 0.68 and 0.8 as the dominant diagonal coefficient, and properly reducing the dominant diagonal coefficient C d The value of (2) can increase the convergence rate;
particularly, the maximum time step when the main diagonal is dominant is obtained and used for calculation, so that the problem that the convergence is poor due to the fact that the main diagonal dominant characteristic of the jacobian matrix A corresponding to the whole linear system cannot be ensured by the time step deltat adopted in the original harmonic balancing method can be effectively solved, and the convergence is effectively accelerated;
step five, solving a new time step D according to the dominant threshold D of the main diagonal ii ；
In this step, a new time step D is solved using the following formula ii ；
Figure BDA0004089110530000043
Step six, a new time step D ii Carrying out flow field calculation in the equation;
in this step, the new time step D solved in step five is applied ii Carrying out flow field calculation in the calculation method in the first step;
in particular, a new time step D ii According to different values of the main diagonal dominant threshold d, the time step delta t after the correction of the original time step delta t and the harmonic balance method HB The switching is carried out, so that most grids adopt the original time step delta t in the calculation process, and the small grids adopt the time step delta t corrected by the harmonic balance method HB Thereby realizing the effect of accelerating convergence.
As shown in FIG. 1, a flow chart of a method for accelerating convergence by a harmonic balance method for predicting unsteady flow of a turbine is implemented in codes, compared with the traditional method for determining a time step mode, the method adds a calculation and judgment process, and needs to calculate the main diagonal dominant threshold D of a jacobian matrix A corresponding to a whole linear system first, and then solve a new time step D according to the difference of the main diagonal dominant threshold D ii And carrying out flow field calculation by using the equation.
The invention relates to a method for accelerating convergence by a harmonic balance method for predicting unsteady flow of a turbine, which is added to a CFD solver of the subject group, wherein a two-dimensional flow sheet at 50% of the height of a blade of a Stage35 compressor is taken as an example, and verification calculation is carried out by adopting the solver added with the method for accelerating convergence by the harmonic balance method for predicting unsteady flow of the turbine. In order to better compare the advantages of the method, two methods are adopted for calculation, namely an original correction method, namely, fixed correction based on frequency and harmonic order is added only in a local time step, the method is called as an HB correction method for convenience of description, and the optimization method based on dominant diagonal is provided, and the Stage35 flow convergence speed calculated by different methods is compared and evaluated.
Specifically, a Stage35 air compressor is selected as a calculation example, a structured grid is adopted to carry out space dispersion on a calculation domain, as shown in fig. 2, the adopted grid is a quasi-three-dimensional grid, the grid topology adopts an O4H topological structure, 3 layers of grids are reserved in the radial direction, and the total number of the grids is 10 ten thousand. The harmonic balance method acceleration convergence method for predicting the unsteady flow of the turbine has no special requirement on a turbulence model, any turbulence model can be used, and the SA turbulence model and the CFL number are adopted in calculation of the embodiment.
In this example, the rotation speed of the rotor is 17189rpm, and the simulation is performed by using a 9-order harmonic balance method, i.e. the correction time step Deltat is calculated HB The value of omega is 1800rad/s, wherein the solver can further carry out dimensionless treatment, the value of N is 9, and the main diagonal is the dominant threshold C d 0.75. For comparison of acceleration effects at different CFL numbers, two CFL numbers of 6 and 12 were selected for calculation, respectively.
As shown in fig. 3, the numbers immediately following the "CFL" in the legend indicate the magnitude of the CFL number used in this example, the suffix indicates that a different correction method is used, where "-HB" indicates the calculation result using the original correction method, and "—new" indicates that the accelerated convergence method proposed by the present invention is added, and any suffix is not added, indicating that no correction method is added.
The calculation result is shown in fig. 3, and a two-dimensional flow sheet convergence speed comparison chart of 50% of the leaf height of the Stage35 compressor calculated by adopting different methods in the embodiment of the invention is provided. When the CFL number is 6, the original time step is directly adopted for calculation without considering any correction algorithm based on harmonic balance, so that calculation divergence is caused, and a convergence result cannot be obtained; when the HB correction method is adopted, and the CFL number is calculated as 6, the HB correction method adds a correction value with a fixed size in the local time step, so that the whole local time step becomes very small, and the convergence is slower; the harmonic balance method for predicting the unsteady flow of the turbine can effectively accelerate the convergence speed, and the convergence speed is increased by 40% compared with the HB correction method by taking the CFL number as 6. With the gradual increase of the CFL number, the convergence speed of the harmonic balance method accelerating convergence method for predicting the unsteady flow of the turbine does not change greatly with the increase of the CFL number. When the CFL number is 6, the highest convergence rate is reached; when the CFL number is increased to 12, the convergence speed of the HB correction method is the fastest, but under the condition that the CFL number is 6, the convergence speed of the HB correction method is still superior to that of the HB correction method, and the harmonic balance method accelerating convergence method for predicting the unsteady flow of the turbine has great advantages in calculating the harmonic balance method.
The foregoing description of the preferred embodiments of the present invention is not intended to be limiting, but rather, any modifications, equivalents, improvements, etc. that fall within the spirit and principles of the present invention are intended to be included within the scope of the present invention.
In summary, the present invention selects the dominant threshold C of the principal diagonal of the jacobian matrix A corresponding to the whole linear system in the numerical calculation process d The maximum time step when the main diagonal is dominant is obtained and used for calculation, so that the problem that the main diagonal dominant characteristic of the jacobian matrix A corresponding to the whole linear system cannot be guaranteed by the time step deltat adopted in the original harmonic balance method, and therefore the convergence is poor is effectively accelerated, a larger gain effect is achieved under a smaller CFL number, the numerical stability is guaranteed, the convergence speed is accelerated, and the accelerated convergence method is provided for turbine flow field prediction in the complex engineering field.

Claims (1)

1. A harmonic balance method accelerating convergence method for predicting unsteady flow of a turbine is characterized by comprising the following steps:
step one, performing space and time dispersion on a control equation;
step two, calculating a local time step;
step three, calculating parameters used by a time pushing format;
step four, calculating a main diagonal dominant threshold d;
step five, solving a new time step D according to the dominant threshold D of the main diagonal ii ；
Step six, the new timeStep length D ii Carrying out flow field calculation in the equation;
(1) the spatially and temporally discretizing the control equation includes:
the space flux is calculated by adopting a space flux calculation method of a traditional steady Reynolds average Navier-Stokes (RANS) method, the time pushing format adopts a BJ-SOR (Block Jacobi Successive Over Relaxation) method, and further, in order to ensure the convergence of calculation, the source item is subjected to implicit discrete, and the calculation method is given by the following formula:
Figure FDA0004089110520000011
AΔW＝-R+C(ΔW)
(L+D+U)ΔW＝-R+C(ΔW)
wherein V represents grid volume, deltat represents physical time step, and I is an identity matrix; r is a residual term, W is a conservation quantity, deltaW is a variation quantity of the conservation quantity, and C (DeltaW) is the dispersion of a harmonic source term; a is a jacobian matrix corresponding to the whole linear system, wherein the jacobian matrix A corresponding to the whole linear system is decomposed into a lower diagonal matrix L, a diagonal matrix D and an upper diagonal matrix U; the time discrete adopts a first-order Euler post-difference format;
the control equation adopted by the BJ-SOR method is given by:
(D+L)ΔWn+1/2＝-[R+C(ΔW)]
(D+U)ΔWn+1＝-[R+C(ΔW)]-LΔWn+1/2
solving a control equation adopted by the BJ-SOR method by adopting multiple cycles, and adopting 4-6 times of calculation in the calculation of the harmonic balancing method to achieve better convergence;
(2) the calculating the local time step includes:
calculating an original time step and a time step corrected by a harmonic balance method:
wherein the calculation of the original time step Δt is given by:
Figure FDA0004089110520000012
the time step delta t corrected by the harmonic balancing method HB The calculation of (2) is given by:
Figure FDA0004089110520000013
wherein lambda is the spectral radius of the jacobian matrix A corresponding to the whole linear system, omega is the angular frequency adopted in the harmonic balancing method, N is the order of the harmonic balancing method, the value of the CFL number is the same as that of the original algorithm, and omega is the grid volume;
(3) the parameters used in the calculation time advance format include:
calculating a diagonal matrix L, a diagonal matrix D and an upper diagonal matrix U which are generated after the BJ-SOR method adopted by the time pushing format is discretized; the diagonal matrix D comprises a diagonal matrix D (delta t) based on the original time step delta t and a time step delta t corrected based on the harmonic balance method HB Diagonal matrix D (Δt) HB ) The calculation method is given by the following formula:
Figure FDA0004089110520000014
where ΔS is the area of a single mesh plane, Λ c Spectral radius for non-sticking flux; the superscript I, J, K indicates different coordinate directions of the structural grid, the subscripts i, j, k indicate index values of grid coordinates, the subscript c indicates non-sticky flux, and the subscript v indicates sticky flux;
(4) the calculating the dominant diagonal threshold d includes:
the maximum time step when the main diagonal is dominant is obtained by selecting the main diagonal dominant threshold d of the jacobian matrix A corresponding to the whole linear system; the dominant diagonal threshold d is given by:
Figure FDA0004089110520000015
wherein C is d Taking a value between 0.68 and 0.8 as the dominant diagonal coefficient, and properly reducing the dominant diagonal coefficient C d The value of (2) can increase the convergence rate;
the maximum time step when the main diagonal is dominant is obtained and used for calculation, so that the problem that the convergence is poor due to the fact that the main diagonal dominant characteristic of the jacobian matrix A corresponding to the whole linear system cannot be ensured by the time step deltat adopted in the original harmonic balance method can be effectively solved, and the convergence is effectively accelerated;
(5) the new time step D is solved according to the dominant threshold D of the main diagonal ii Comprising the following steps:
solving for a new time step D using ii ；
Figure FDA0004089110520000021
(6) Said new time step D ii In the carry-over equation, performing flow field calculations includes:
the new time step D solved in (5) ii Carrying out flow field calculation in the calculation method in the step (1);
in particular, the new time step D ii According to different values of the main diagonal dominant threshold d, the time step delta t after the correction of the original time step delta t and the harmonic balance method is carried out HB The switching is carried out, so that most grids adopt the original time step delta t in the calculation process, and the small grids adopt the time step delta t corrected by the harmonic balancing method HB Thereby realizing the effect of accelerating convergence.
CN202310145756.4A 2023-02-22 2023-02-22 Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine Active CN116244853B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310145756.4A CN116244853B (en) 2023-02-22 2023-02-22 Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310145756.4A CN116244853B (en) 2023-02-22 2023-02-22 Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine

Publications (2)

Publication Number Publication Date
CN116244853A true CN116244853A (en) 2023-06-09
CN116244853B CN116244853B (en) 2023-10-20

Family

ID=86625663

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310145756.4A Active CN116244853B (en) 2023-02-22 2023-02-22 Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine

Country Status (1)

Country Link
CN (1) CN116244853B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8473533B1 (en) * 2010-06-17 2013-06-25 Berkeley Design Automation, Inc. Method and apparatus for harmonic balance using direct solution of HB jacobian
CN114021322A (en) * 2021-10-28 2022-02-08 华中科技大学 Modal mode shape analysis method and device of linear periodic time-varying system
CN115186608A (en) * 2022-07-12 2022-10-14 北京航空航天大学 Grid self-adaptive turbulence simulation method based on turbulence energy spectrum coupling RSM model
CN115329689A (en) * 2022-07-05 2022-11-11 北京航空航天大学 High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8473533B1 (en) * 2010-06-17 2013-06-25 Berkeley Design Automation, Inc. Method and apparatus for harmonic balance using direct solution of HB jacobian
CN114021322A (en) * 2021-10-28 2022-02-08 华中科技大学 Modal mode shape analysis method and device of linear periodic time-varying system
CN115329689A (en) * 2022-07-05 2022-11-11 北京航空航天大学 High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion
CN115186608A (en) * 2022-07-12 2022-10-14 北京航空航天大学 Grid self-adaptive turbulence simulation method based on turbulence energy spectrum coupling RSM model

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王子维;范召林;江雄;陈逖;: "用于多级压气机非定常流动的谐波平衡方法研究", 推进技术, no. 11 *

Also Published As

Publication number Publication date
CN116244853B (en) 2023-10-20

Similar Documents

Publication Publication Date Title
He et al. Robust aerodynamic shape optimization—from a circle to an airfoil
Morgenthal et al. An immersed interface method for the Vortex-In-Cell algorithm
Kersken et al. Time-linearized and time-accurate 3D RANS methods for aeroelastic analysis in turbomachinery
Bosnyakov et al. Computational tools for supporting the testing of civil aircraft configurations in wind tunnels
CN112231847B (en) Transition position determining method and device, electronic equipment and storage medium
Cheng et al. A surface parametric control and global optimization method for axial flow compressor blades
Li et al. Aeroelastic global structural optimization using an efficient CFD-based reduced order model
Li et al. Multipoint and multiobjective optimization of a centrifugal compressor impeller based on genetic algorithm
CN112417773B (en) Multidisciplinary optimization design method, device and equipment of multistage axial flow expander
Gang et al. Mesh deformation on 3D complex configurations using multistep radial basis functions interpolation
Ning MAP: A CFD package for turbomachinery flow simulation and aerodynamic design optimization
Van Zante et al. The influence of compressor blade row interaction modeling on performance estimates from time-accurate, multistage, navier–stokes simulations
Tang et al. Natural laminar flow shape optimization in transonic regime with competitive Nash game strategy
Schwerdt et al. A model reduction method for bladed disks with large geometric mistuning using a partially reduced intermediate system model
Pierret Multi-objective and multi-disciplinary optimization of three-dimensional turbomachinery blades
Hu et al. Adjoint-based an adaptive finite volume method for steady Euler equations with non-oscillatory k-exact reconstruction
Xu et al. Robust Newton–Krylov adjoint solver for the sensitivity analysis of turbomachinery aerodynamics
CN116244853B (en) Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine
Ellbrant et al. CFD optimization of a transonic compressor using multiobjective GA and metamodels
Liu et al. Redesign of axial fan using viscous inverse design method based on boundary vorticity flux diagnosis
Huang et al. Parallel preconditioned WENO scheme for three-dimensional flow simulation of NREL Phase VI Rotor
He et al. Parametric study of the aeroelastic response of mistuned bladed disks
CN114065423B (en) Method for rapidly evaluating flutter of fan blade of aircraft engine
Adibi et al. Predicting airfoil stalling dynamics using upwind numerical solutions to non-viscous equations
Fan et al. Power system equilibrium tracing and bifurcation detection based on the continuation of the recursive projection method

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