CN114091376B - High-precision reconstruction correction shock wave capturing method based on subunit weighting format - Google Patents

High-precision reconstruction correction shock wave capturing method based on subunit weighting format Download PDF

Info

Publication number
CN114091376B
CN114091376B CN202210069705.3A CN202210069705A CN114091376B CN 114091376 B CN114091376 B CN 114091376B CN 202210069705 A CN202210069705 A CN 202210069705A CN 114091376 B CN114091376 B CN 114091376B
Authority
CN
China
Prior art keywords
flux
shock wave
format
points
point
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.)
Active
Application number
CN202210069705.3A
Other languages
Chinese (zh)
Other versions
CN114091376A (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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202210069705.3A priority Critical patent/CN114091376B/en
Publication of CN114091376A publication Critical patent/CN114091376A/en
Application granted granted Critical
Publication of CN114091376B publication Critical patent/CN114091376B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force 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)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a high-precision reconstruction correction shock wave capturing method based on a subunit weighting format, which belongs to the technical field of shock wave capturing and comprises the following steps: s1, calculating the distribution of solution points and flux points of the high-order CPR and the flux derivative; s2, detecting the units with shock wave by shock wave detection method; s3, dividing the unit into subunits based on the flux points, and constructing a format with shock wave capturing capability on the subunits; s4, determining CPR, HNNW and calculation units of second-order format based on the value of the shock wave detection factor in the shock wave detection method; s5, the Riemann flux at the cell interface is calculated by interpolation polynomials corresponding to the calculation formats of the adjacent cells. The invention not only saves the calculation amount, but also can remarkably improve the capturing capability of the format on the shock wave, and can ensure that the final algorithm still has higher resolution near the shock wave while reducing the requirement of the algorithm on the detection accuracy.

Description

High-precision reconstruction correction shock wave capturing method based on subunit weighting format
Technical Field
The invention relates to the technical field of shock wave capturing, in particular to a high-precision reconstruction correction shock wave capturing method based on a subunit weighting format.
Background
Computational Fluid Dynamics (CFD) is one of the important means for developing the research of Fluid mechanics mechanism, and plays an increasingly important role in the design and performance evaluation of aerospace aircrafts. With continuous refinement of design in engineering application, people have higher and higher precision requirements on CFD calculation results, and particularly when the precision flow field structure is simulated, a second-order precision algorithm is difficult to meet the requirements. Higher order precision algorithms require less computation than second order precision algorithms to achieve the same error level. The high-precision algorithm has the advantages of simulating a fine flow field structure, such as the advantages in various aspects of turbulence large vortex simulation, separation vortex simulation, direct simulation, pneumatic acoustic calculation and the like. Because of the method of Correction process via Reconstruction (CPR), when selecting a special Correction function, the method can be equivalent to a Discrete Galileo (DG) method, a Spectral Difference (SD) method or a Spectral Volume (SV) method, and the calculation amount is relatively small, which has received great attention in recent years.
In recent years, CPR methods have been developed in terms of free flow conservation, nonlinear stability, large vortex simulation, shock wave capture, and the like. However, the CPR method still has problems in shock wave capture, and further development is required. One method that has gained much attention is to mix the higher order format with the robust weighted class/lower order shock wave format, which can take full advantage of each format.
The difficulty of this shock wave capture method of CPR method is represented by: solution points (coordinate points stored by physical quantity when an equation is discrete) of the CPR method are distributed in a non-equidistant mode (distances between adjacent solution points at different positions are unequal), when the solution points are mixed with a traditional high-order finite difference format with better shock wave capturing capability and developed based on the equidistant solution points, the solution points of the two types of formats are not matched, and data conversion is required to be carried out repeatedly during format switching; secondly, when the CPR method is mixed with a second-order shock wave capturing format, the resolution ratio is greatly reduced due to the fact that the detection is not accurate enough near the shock wave, and how to guarantee the shock wave capturing and the high-resolution characteristic is difficult.
Disclosure of Invention
The invention aims to overcome the defects of the prior art, provides a high-precision reconstruction correction shock wave capturing method based on a subunit weighting format, solves the problems in the background art, realizes the numerical simulation of the high-order CPR method on the flow of shock waves in computational fluid mechanics, not only saves the calculated amount, but also can obviously improve the capturing capability of the format on the shock waves, and can ensure that the final algorithm still has higher resolution near the shock waves while reducing the requirement of the algorithm on the detection accuracy.
The purpose of the invention is realized by the following scheme:
a high-precision reconstruction correction shock wave capturing method based on a subunit weighting format comprises the following steps:
s1, distributing solution points on grid cells as Gauss-Legendre integral points, wherein flux points are one more than solution points according to the number of flux points in each dimension, the distance between adjacent flux points is Gauss integral weight, and the flux points at two ends are distributed on a cell interface; the flux derivative of the high-order CPR method is obtained by derivation based on Lagrange interpolation polynomial of the flux point, and the flux derivative of the solution point is corrected through Radau polynomial by the difference between the Riemann flux and the unit discontinuous flux;
s2, detecting the units with shock wave by shock wave detection method;
s3, dividing the unit into sub-units based on the flux points, and constructing a format with shock wave capturing capacity on the sub-units: the first format is a high-order hybrid non-equidistant nonlinear weighting format HNNW, and flux derivatives are calculated based on a differential operator of a solution point flux point hybrid; the flux at the flux point is Riemann flux, and the left and right values of the physical quantity at the flux point are obtained through non-equidistant nonlinear weighted interpolation based on the solution point; the second format is a second order format, calculating flux derivatives based on two flux points adjacent to the solution point; the left and right values of the physical quantity at the flux point are obtained through restricted first-order polynomial distribution;
s4, determining CPR, HNNW and calculation units of second-order format based on the value of the shock wave detection factor in the shock wave detection method;
s5, the Riemann flux at the cell interface is calculated by interpolation polynomials corresponding to the calculation formats of the adjacent cells.
Further, in step S1, the mesh cells include quadrilateral mesh cells of a two-dimensional structural mesh.
Further, in step S2, the shock wave detection method is specifically a shock wave detection method based on the highest modal attenuation.
Further, in step S3, the calculating the flux derivative based on the difference operator of the solution point flux point hybrid type includes the sub-steps of: two flux points adjacent to the solution point and the solution point outside the flux point are used to form the template point for calculating the flux derivative.
Further, in step S4, the calculation unit includes a problem unit, and the problem unit adopts a high-order non-equidistant non-linear weighting format HNNW and/or a second-order format calculation.
Further, comprising the sub-steps of: and adopting a second-order format for the problem units in the small area containing the shock waves, and adopting a high-order HNNW format for the rest problem units.
The beneficial effects of the invention include:
because the solution point of the HNNW and the second-order format is coincident with the solution point of the CPR, the shock wave capturing format of the CPR method based on the sub-unit HNNW and the second-order format limit does not need data conversion, thereby saving the calculation amount.
The HNNW format with strong shock wave capturing capability and the second-order format with strong shock wave capturing capability are adopted, so that the shock wave capturing capability of the formats can be obviously improved, and meanwhile, the HNNW still has high-order precision in a smooth area, so that the requirement of the algorithm on the detection accuracy is greatly reduced, and finally, the algorithm still has high resolution near the shock wave.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to these drawings without creative efforts.
FIG. 1 is a diagram of the distribution of solution points and flux points on a quadrilateral grid cell in an embodiment of the present invention;
fig. 2 is a template point for the fifth-order HNNW5 to calculate the flux derivative in an embodiment of the present invention;
FIG. 3 is a block diagram of the template points for second order format calculation of flux derivatives in an embodiment of the invention.
Detailed Description
All features disclosed in all embodiments in this specification, or all methods or process steps implicitly disclosed, may be combined and/or expanded, or substituted, in any way, except for mutually exclusive features and/or steps.
The technical problems, technical concepts, working principles and working processes to be solved by the present invention are further described in detail and fully with reference to the accompanying drawings 1 to 3.
The design concept of the invention is as follows: the scheme provided by the invention aims to provide a shock wave capturing technical scheme based on a subunit non-equidistant weighting format of a high-order CPR method, and can realize numerical simulation of the high-order CPR method on flow including shock waves in computational fluid mechanics. In order to achieve the purpose, the scheme provided by the invention firstly constructs a Hybrid non-equidistant non-linear Weighted (HNNW) format and a second-order format of high-order solution point flux points based on solution points which are distributed in a non-equidistant way in a calculation space, and then divides a calculation unit into a CPR unit, an HNNW unit and a second-order format unit by a shock detection method, so that a problem unit in a small region containing shock waves is calculated by adopting the second-order format, and the rest problem units are calculated by adopting the HNNW format with high resolution characteristic and shock wave capturing capability, thereby realizing effective shock wave capturing and high resolution.
In specific application, the technical scheme of the invention is explained in detail as follows:
(1) on a quadrilateral grid unit of a two-dimensional structure grid, the distribution solution points are Gauss-Legendre integral points, flux points are one more than solution points according to the number of the flux points in each dimension, the distance between adjacent flux points is Gauss integral weight, the flux points at two ends are distributed on a unit interface, the condition that the number of solution points in each direction is 5 is given in figure 1, the solution points are circular points in the figure, and the flux points are square points in the figure. The flux derivative of the high-order CPR method is obtained by derivation based on Lagrange interpolation polynomial of the flux point, and the flux derivative of the solution point is corrected through Radau polynomial according to the difference between the Riemann flux and the unit discontinuous flux.
(2) The shock detection technique selects the detection method based on the highest modal attenuation (MD). Units in which shock waves are likely to exist are detected by a shock wave detection technology, and are called problem units.
(3) The cell is split into sub-cells based on the flux points, as is the grey background cell in fig. 1. A format with shock wave capture capability is constructed on the subunit. The first format is a high-order non-equidistant non-linear weighting format (HNNW), which calculates the flux derivative based on a difference operator of a solution point flux point hybrid, that is, two flux points adjacent to the solution point and a solution point outside the flux point are used to form a template point for calculating the flux derivative. Using fifth order CPR as an example, fig. 2 shows the template points of the calculated flux derivative for the second solution point of a cell. The flux at the flux point is Riemann flux, and the left and right values of the physical quantity at the flux point are obtained through non-equidistant nonlinear weighted interpolation based on the solution point. The second format is a second order format, where the flux derivative is calculated based on two flux points adjacent to the solution point, as shown in fig. 3. The left and right values of the physical quantity at the flux point are obtained by a restricted first order polynomial distribution.
(4) The problem unit is computed using a high-order HNNW and/or a second-order format. In a specific embodiment, the calculation can be divided into three regions from small to large based on the value of the shock detection factor, and the calculation is performed by using CPR, HNNW and second-order formats respectively. The principle is that a second-order format is adopted for problem units in a small area containing shock waves, and a high-order HNNW format is adopted for other problem units.
(5) The Riemann flux at the cell interface is calculated by interpolation polynomials corresponding to the calculation formats of adjacent cells.
In the embodiment of the invention, the provided brand-new mixed high-order non-equidistant nonlinear weighting format of the subunits improves and guarantees the high resolution of the final format, and data conversion is not needed between the main unit and the subunits. And the problem unit is divided into two formats for calculation based on the shock wave detection factor, a second-order format is adopted near a shock wave region to ensure the robustness of the final format capture, and a high-resolution characteristic of the final format is ensured by the problem unit outside the shock wave by adopting a high-order HNNW.
Example 1: a high-precision reconstruction correction shock wave capturing method based on a subunit weighting format comprises the following steps:
s1, distributing solution points on grid cells as Gauss-Legendre integral points, wherein flux points are one more than solution points according to the number of flux points in each dimension, the distance between adjacent flux points is Gauss integral weight, and the flux points at two ends are distributed on a cell interface; the flux derivative of the high-order CPR method is obtained by derivation based on Lagrange interpolation polynomial of the flux point, and the flux derivative of the solution point is corrected through Radau polynomial by the difference between the Riemann flux and the unit discontinuous flux;
s2, detecting the units with shock wave by shock wave detection method;
s3, dividing the unit into sub-units based on the flux points, and constructing a format with shock wave capturing capacity on the sub-units: the first format is a high-order non-equidistant non-linear weighting format HNNW, and flux derivatives are calculated based on a differential operator of a solution point flux point mixed type; the flux at the flux point is Riemann flux, and the left and right values of the physical quantity at the flux point are obtained through non-equidistant nonlinear weighted interpolation based on the solution point; the second format is a second order format, calculating flux derivatives based on two flux points adjacent to the solution point; the left and right values of the physical quantity at the flux point are obtained through restricted first-order polynomial distribution;
s4, determining CPR, HNNW and calculation units of second-order format based on the value of the shock wave detection factor in the shock wave detection method;
s5, the Riemann flux at the cell interface is calculated by interpolation polynomials corresponding to the calculation formats of the adjacent cells.
Example 2: on the basis of embodiment 1, in step S1, the mesh cells include quadrangular mesh cells of a two-dimensional structural mesh.
Example 3: in step S2, based on embodiment 1, the shock wave detection method is specifically a shock wave detection method based on the highest modal attenuation.
Example 4: on the basis of embodiment 1, in step S3, the calculating the flux derivative based on the difference operator of the solution point flux point mixture type includes the sub-steps of: two flux points adjacent to the solution point and the solution point outside the flux point are used to form the template point for calculating the flux derivative.
Example 5: on the basis of embodiment 1, in step S4, the calculation unit includes a problem unit, and the problem unit adopts a high-order non-equidistant non-linear weighting format HNNW and/or a second-order format calculation.
Example 6: on the basis of the embodiment 5, the method comprises the following substeps: and adopting a second-order format for the problem units in the small area containing the shock waves, and adopting a high-order HNNW format for the rest problem units.
The parts not involved in the present invention are the same as or can be implemented using the prior art.
The above-described embodiment is only one embodiment of the present invention, and it will be apparent to those skilled in the art that various modifications and variations can be easily made based on the application and principle of the present invention disclosed in the present application, and the present invention is not limited to the method described in the above-described embodiment of the present invention, so that the above-described embodiment is only preferred, and not restrictive.
Other embodiments than the above examples may be devised by those skilled in the art based on the foregoing disclosure, or by adapting and using knowledge or techniques of the relevant art, and features of various embodiments may be interchanged or substituted and such modifications and variations that may be made by those skilled in the art without departing from the spirit and scope of the present invention are intended to be within the scope of the following claims.

Claims (4)

1. A high-precision reconstruction correction shock wave capturing method based on a subunit weighting format is characterized by comprising the following steps:
s1, distributing solution points on grid cells as Gauss-Legendre integral points, wherein flux points are one more than solution points according to the number of flux points in each dimension, the distance between adjacent flux points is Gauss integral weight, and the flux points at two ends are distributed on a cell interface; the flux derivative of the high-order CPR method is obtained by derivation based on Lagrange interpolation polynomial of the flux point, and the flux derivative of the solution point is corrected through Radau polynomial by the difference between the Riemann flux and the unit discontinuous flux;
s2, detecting the units with shock wave by shock wave detection method;
s3, dividing the unit into sub-units based on the flux points, and constructing a format with shock wave capturing capacity on the sub-units: the first format is a high-order hybrid non-equidistant nonlinear weighting format HNNW, and flux derivatives are calculated based on a differential operator of a solution point flux point hybrid; the flux at the flux point is Riemann flux, and the left and right values of the physical quantity at the flux point are obtained through non-equidistant nonlinear weighted interpolation based on the solution point; the second format is a second order format, calculating flux derivatives based on two flux points adjacent to the solution point; the left and right values of the physical quantity at the flux point are obtained through restricted first-order polynomial distribution;
s4, determining CPR, HNNW and calculation units of second-order format based on the value of the shock wave detection factor in the shock wave detection method; in step S4, the calculating unit includes a problem unit, and the problem unit adopts a high-order non-equidistant non-linear weighting format HNNW and/or a second-order format calculation; and comprises the substeps of: adopting a second-order format for the problem units in the region with small shock waves, and adopting a high-order HNNW format for the rest problem units;
s5, the Riemann flux at the cell interface is calculated by interpolation polynomials corresponding to the calculation formats of the adjacent cells.
2. The method for high-precision reconstruction-modified shock wave capturing based on sub-unit weighting format according to claim 1, wherein in step S1, the grid cells comprise quadrilateral grid cells of a two-dimensional structural grid.
3. The method for high-precision reconstruction-modified shock wave capturing based on subunit weighting format according to claim 1, wherein in step S2, the shock wave detection method is a shock wave detection method based on highest modal attenuation.
4. The method for high-precision reconstruction-modified shock wave capturing based on subunit weighting format according to claim 1, wherein in step S3, the calculating the flux derivative based on the difference operator of solution point flux point hybrid type includes the sub-steps of: two flux points adjacent to the solution point and the solution point outside the flux point are used to form the template point for calculating the flux derivative.
CN202210069705.3A 2022-01-21 2022-01-21 High-precision reconstruction correction shock wave capturing method based on subunit weighting format Active CN114091376B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210069705.3A CN114091376B (en) 2022-01-21 2022-01-21 High-precision reconstruction correction shock wave capturing method based on subunit weighting format

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210069705.3A CN114091376B (en) 2022-01-21 2022-01-21 High-precision reconstruction correction shock wave capturing method based on subunit weighting format

Publications (2)

Publication Number Publication Date
CN114091376A CN114091376A (en) 2022-02-25
CN114091376B true CN114091376B (en) 2022-04-12

Family

ID=80309034

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210069705.3A Active CN114091376B (en) 2022-01-21 2022-01-21 High-precision reconstruction correction shock wave capturing method based on subunit weighting format

Country Status (1)

Country Link
CN (1) CN114091376B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114282462B (en) * 2022-03-03 2022-05-17 中国空气动力研究与发展中心计算空气动力研究所 Method for capturing shock wave by coupling flow and grid distribution reconstruction limiter
CN116542184B (en) * 2023-07-05 2023-09-19 中国空气动力研究与发展中心计算空气动力研究所 Method and device for calculating viscosity flux, terminal equipment and storage medium
CN116882322B (en) * 2023-09-06 2024-02-13 中国空气动力研究与发展中心计算空气动力研究所 Calculation method and device for non-sticky flux, terminal equipment and storage medium

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108197072B (en) * 2017-12-27 2021-02-02 中国空气动力研究与发展中心计算空气动力研究所 High-precision intermittent Galerkin artificial viscous shock wave capturing method based on weighted conservative variable step
CN108153984B (en) * 2017-12-27 2021-04-13 中国空气动力研究与发展中心计算空气动力研究所 High-precision Galois gold artificial viscous shock wave capturing method based on flow field density step
CN108197367B (en) * 2017-12-27 2021-07-27 中国空气动力研究与发展中心计算空气动力研究所 High-precision intermittent Galerkin artificial viscous shock wave capturing method
CN111159956B (en) * 2019-12-10 2021-10-26 北京航空航天大学 Feature-based flow field discontinuity capturing method

Also Published As

Publication number Publication date
CN114091376A (en) 2022-02-25

Similar Documents

Publication Publication Date Title
CN114091376B (en) High-precision reconstruction correction shock wave capturing method based on subunit weighting format
CN110069800B (en) Three-dimensional structure topology optimization design method and equipment with smooth boundary expression
CN112100835B (en) High-efficiency high-precision airfoil-shaped flow numerical simulation method suitable for complex flow
CN116151084B (en) Simulation method and device based on structural grid, terminal equipment and storage medium
CN106682262B (en) Numerical simulation method for obtaining aircraft flow field
CN110929443B (en) Two-dimensional flood simulation method based on high-precision terrain generalization
CN109726465B (en) Three-dimensional non-adhesive low-speed streaming numerical simulation method based on non-structural curved edge grid
CN107038308B (en) A kind of regular grid terrain modeling method based on linear interpolation
CN114970865B (en) Quantum circuit processing method and device on quantum chip and electronic equipment
CN109738972B (en) Air pollutant forecasting method and device and electronic equipment
CN111507026A (en) Dual-grid multi-scale finite element method for simulating node Darcy permeation flow rate
CN115688607A (en) Band-shaped cross-sea quasi-geoid refinement method based on multi-channel spectrum combination
CN110147646B (en) Over-current processing method for linear water retaining structure under numerical simulation framework
CN113378440A (en) Fluid-solid coupling numerical simulation calculation method, device and equipment
CN116542184B (en) Method and device for calculating viscosity flux, terminal equipment and storage medium
CN111159956B (en) Feature-based flow field discontinuity capturing method
CN107563080A (en) Two-phase medium stochastic model parallel generation method, electronic equipment based on GPU
CN107797132A (en) A kind of inversion method of three dimensional radiation field dosage
CN116451493A (en) Fluid simulation method for complex structure
CN111898819B (en) Space grid dividing method and device
CN108197368A (en) It is a kind of to be suitable for the geometrical constraint of aircraft complexity aerodynamic configuration and weight function Two Simple Methods
CN114611423A (en) Three-dimensional multiphase compressible fluid-solid coupling rapid calculation method
CN110674581B (en) Method and system for accurately judging consistency of digital twin model
CN106569027A (en) Accurate algorithm for correcting electric quantity and power factor under three-phase three-wire electric-energy-meter false wiring
CN112632825B (en) Electrostatic field smooth finite element numerical algorithm based on finite element super-convergence

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