Ship stack mold resistance CFD calculation method suitable for different grid masses

Abstract

The invention discloses a ship stack mold resistance CFD calculation method adapting to different grid masses, and relates to the technical field of ship CFD. According to the method, a non-orthogonal correction super-relaxation method is adopted for a diffusion term in a momentum discrete equation, the influence of poor grid orthogonality on a calculation result is considered, a convection term adopts a second-order windward format containing a gradient limiter, and a non-orthogonal correction super-relaxation method is adopted for a non-orthogonal correction super-relaxation method. The problem of extreme value of interface speed reconstruction caused by uneven grid density and large slope rate can be effectively avoided, in addition, during solving based on the SIMPLE algorithm, a dual sub-relaxation method combining global implicit sub-relaxation with local explicit sub-relaxation is adopted, the problem of divergence of an iteration process under the condition of low-quality grids can be effectively avoided, and the method is suitable for solving the interface speed. The method can effectively avoid the influence of low-mass grids such as violent scale change, large slope rate and poor orthogonality, ensures the calculation stability under the condition of the low-mass grids, and is beneficial to improving the stability and calculation precision of ship stack mold resistance CFD calculation.

Description

Technical Field

[0001] This application relates to the technical field of ship CFD, and in particular to a ship superposition mold resistance CFD calculation method adapted to different mesh qualities. Background Art

[0002] CFD (Computational Fluid Dynamics) technology can obtain the discrete distribution of the flow field of fluid flow on a continuous region by numerically solving the differential equations that control fluid flow, thereby approximately simulating the fluid flow situation. The ship superposition mold resistance CFD calculation can use CFD technology to achieve resistance calculation, which is an important technology in the field of ship research.

[0003] The ship superposition mold resistance CFD calculation depends on mesh generation. Mesh quality is the basis of the ship superposition mold resistance CFD calculation, which will have a direct and important impact on the calculation accuracy and efficiency, and even become the key to the success or failure of the ship superposition mold resistance CFD calculation. However, due to the very complex geometric shape of the object processed in the ship superposition mold resistance CFD calculation, there will inevitably be low-quality meshes during mesh generation:

[0004] When using traditional structured meshes and unstructured meshes, taking surface ships as an example, due to the complex three-dimensional space curve of the hull geometry, there will inevitably be low-quality meshes with high distortion rates and low orthogonality during the mesh generation process. These meshes usually appear in the areas near the bow and stern of the ship, where the curvature changes greatly, the spatial flow is complex, the physical quantities such as velocity and pressure change violently in space and time, and have a huge or even decisive impact on the calculation results. The hybrid mesh generated by the Cartesian mesh through the cutting element method combined with the body-fitted prism layer mesh to process the physical surface boundary layer is also becoming more and more popular, representing the current trend of mesh technology. However, at the junction of the Cartesian mesh and the prism layer mesh, there are often problems such as a sharp change in the scale of a small number of mesh elements and a large mesh inclination rate.

[0005] Therefore, in the ship superposition mold resistance CFD calculation scenario, due to the difficulty in ensuring the quality of the divided mesh elements and the inevitable existence of low-quality mesh elements, this not only affects the numerical calculation accuracy and stability, but also causes calculation divergence, affecting the ship superposition mold resistance CFD calculation effect. There are some practices that will appropriately reduce the mesh size and increase the mesh density in places with more complex geometries. This method solves the above problems by improving the mesh quality. This method will greatly increase the workload of mesh generation, with a large implementation difficulty and an unsatisfactory actual improvement effect. Summary of the Invention

[0006] In view of the above problems and technical requirements, this application proposes a CFD calculation method for ship stack resistance that is suitable for different grid qualities. The technical solution of this application is as follows:

[0007] A CFD calculation method for ship stacking resistance adapted to different grid qualities, the ship stacking resistance CFD calculation method comprising:

[0008] The incompressible RANS equations including the continuity equation and the momentum equation are used as the control equations. The control volume unit is taken as the grid unit obtained by dividing the computational domain. The control equations are discretized using the finite volume method to obtain the discrete control equations of each control volume unit.

[0009] The second-order upwind scheme is used for the convection term in the discrete control equation of the control volume unit, and the super-relaxation non-orthogonal correction method is used for the diffusion term, and the discrete control algebraic equation group of the control volume unit is obtained;

[0010] The turbulence model is used to close the discrete control algebraic equations, and based on the SIMPLE algorithm, the global implicit sub-relaxation method combined with the local explicit sub-relaxation method is used to perform numerical calculations on the discrete control algebraic equations of any control volume unit in the computational domain to obtain the velocity field, pressure field and flow rate on the surface of the control volume unit.

[0011] The beneficial technical effects of this application are:

[0012] The present application discloses a CFD calculation method for ship stack resistance that is suitable for different mesh qualities. The method adopts an incompressible RANS equation group as a control equation and uses a finite volume method to discretize the calculation domain. A non-orthogonal modified super-relaxation method is used for the diffusion term in the momentum discrete equation, and the influence of poor mesh orthogonality on the calculation result is considered. The convection term adopts a second-order upwind scheme with a gradient limiter, which can effectively avoid the extreme value problem of interface velocity reconstruction caused by uneven mesh density and large inclination. In addition, when solving based on the SIMPLE algorithm, a double sub-relaxation method of global implicit sub-relaxation combined with local explicit sub-relaxation is adopted, which can effectively avoid the divergence problem of the iterative process in the case of low-quality meshes, improve the calculation stability, and calculate the prediction speed, which will not affect the final calculation accuracy. The method can effectively avoid the influence of low-quality meshes such as drastic scale changes, large inclination and poor orthogonality, ensure the stability of the calculation in the case of low-quality meshes, and is conducive to improving the stability and calculation accuracy of the CFD calculation of ship stack resistance, and can effectively avoid the divergence of the iterative process and accelerate convergence, so that the accuracy and efficiency of numerical simulation can be optimized without additional mesh encryption processing on the basis of existing mesh division. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1It is a schematic flow chart of the CFD calculation method for the resistance of stacked ship molds in an embodiment of the present application. Detailed implementation manners

[0014] The following further describes the detailed implementation manners of the present application with reference to the accompanying drawings.

[0015] The present application discloses a CFD calculation method for the resistance of stacked ship molds adaptable to different grid qualities. This method uses the incompressible RANS equations as the governing equations. After discretizing the computational domain based on the finite volume method, a turbulence model is used to close the RANS equations and the SIMPLE algorithm is used to implement the numerical calculation of the flow field to obtain the velocity field, pressure field, and flow rate at the surface. On the framework of this method, the present application mainly optimizes the governing equations and the iterative calculation process based on the SIMPLE algorithm, so as to be adaptable to different grid qualities and have good numerical calculation effects in low-quality grid scenarios. The introduction is as follows:

[0016] I. Discretize the governing equations

[0017] The present application uses the incompressible RANS equations as the governing equations as follows:

[0018]

[0019] Among them, is called the continuity equation, is called the momentum equation.

[0020] In the above formula, v is the velocity vector of the flow field, p is the pressure, ρ is the fluid density, μ is the fluid viscosity, and μ t is the turbulent viscosity. is the gradient operator.

[0021] Take the control volume element as the grid element itself obtained by dividing the computational domain, and combine the collocated grid technique, that is, store both the scalar field and the vector field in the flow field at the center position of the grid element. Use the finite volume method to discretize the governing equation (1) to obtain the discrete governing equations of each control volume element.

[0022] The discretization of this step includes discretizing the momentum equation and the continuity equation respectively:

[0023] (a) Use the finite volume method to discretize the momentum equation to obtain the momentum discrete governing equation of the control volume element P, written as:

[0024]

[0025] Among them, is the unsteady term, is the convection term, is the diffusion term, p f Sf is the source term.

[0026] In the above equation, nb(P) is the set composed of the surfaces of control volume cell P, and f represents any surface of control volume cell P. (ρvv) f represents the value of ρvv at the surface f of control volume cell P, S f represents the area of the surface f of control volume cell P, p f represents the pressure at the surface f of control volume cell P, represents at control volume cell P value, V P represents the volume of control volume cell P.

[0027] (b) Discretize the continuity equation using the finite volume method to obtain the discrete control equation of continuity for control volume cell P, written as:

[0028]

[0029] II. Rearrange the discrete control equation of the control volume cell to obtain the discrete control algebraic equation set of the control volume cell. This step introduces (2) and (3) in the discrete control equation as follows:

[0030] (a) Rearrange the momentum discrete control equation shown in the above equation (2) to obtain the velocity algebraic equation set:

[0031] In order to make this method better adapt to grid cells of various qualities and eliminate the influence of low-quality grids as much as possible, a key point in this step when obtaining the discrete control algebraic equation set is to use the second-order upwind scheme for the convection term in the momentum discrete control equation and the over-relaxation non-orthogonal correction method for the diffusion term. Specifically:

[0032] ① Represent the convection term in the momentum discrete control equation of control volume cell P using the second-order upwind scheme as:

[0033]

[0034] Among them, the velocity field at the surface f of control volume cell P at the (n + 1)-th iteration step represents the velocity field of the upstream cell UP at the (n + 1)-th iteration step, represents the velocity field of the upstream cell UP of control volume cell P at the n-th iteration step, represents the connection vector between the surface center of the surface f of control volume cell P and the center of the upstream cell UP; represents the interface flux of the surface f of control volume cell P at the n-th iteration step. Among them, the upstream cell UP is determined according to the upwind scheme according to the current general method, which will not be elaborated here.

[0035] ② The diffusion term in the momentum discrete control equation for the control volume cell P is expressed by the over-relaxation non-orthogonal correction method as:

[0036]

[0037] Among them, is the orthogonal linearization term, is the non-orthogonal non-linearization term; using the over-relaxation non-orthogonal correction method, Ε f can be expressed as T f = S f - E f , where e is the unit vector connecting the control volume cell P and its adjacent control volume cell N.

[0038] In addition, corresponding treatments are also performed on the unsteady term and the source term, including:

[0039] ③ The unsteady term in the momentum discrete control equation for the control volume cell P is expressed by the first-order implicit format as:

[0040]

[0041] Among them, ρ P is the fluid density at the control volume cell P, v n+1 is the velocity of the flow field at the (n + 1)-th iteration step, v n is the velocity of the flow field at the n-th iteration step, and Δt is the time interval between two adjacent iteration steps.

[0042] ④ The pressure p f in the source term of the momentum discrete control equation for the control volume cell P is expressed by the second-order central difference format as:

[0043]

[0044] Among them, p P represents the pressure of the control volume cell P, p N represents the pressure of the control volume cell N, d Pf represents the connection vector between the surface center of the surface f of the control volume cell P and the center of the control volume cell P, and d Nf represents the connection vector between the surface center of the surface f of the control volume cell P and the center of the control volume cell N.

[0045] Finally, by integrating the above four representation methods, the velocity algebraic equations of the control volume cell P obtained by organizing the momentum discrete control equation (2) of the control volume cell P are as follows:

[0046]

[0047] Among them, represents the main diagonal coefficient of the coefficient matrix of the velocity algebraic equations, represents the off-diagonal coefficient of the coefficient matrix of the velocity algebraic equations. NB(P) is the set composed of the control volume cells adjacent to the control volume cell P, and N represents any control volume cell adjacent to the control volume cell P. v P represents the velocity field of the control volume cell P, v N represents the velocity field of the control volume cell N, represents the source term vector of the velocity algebraic equations.

[0048] (b) Rearranging the continuity discrete control equation shown in the above formula (3) to obtain the pressure correction algebraic equations, including adopting the over-relaxation non-orthogonal correction method for the diffusion discrete terms, and the pressure correction algebraic equation for the control volume cell P is sorted out as:

[0049]

[0050] Among them, represents the main diagonal coefficient of the coefficient matrix of the pressure correction algebraic equations, represents the off-diagonal coefficient of the coefficient matrix of the pressure correction algebraic equations. p′ P represents the pressure correction term of the control volume cell P, p′ N represents the pressure correction term of the control volume cell N, represents the source term of the pressure correction algebraic equations

[0051] III. Closing the discrete control algebraic equations by using a turbulence model

[0052] In this application, the RNG k-ε model is selected, and its turbulent viscosity μ t is defined as:

[0053]

[0054] Among them, C μ is a constant, generally taken as 0.09. k represents the turbulent kinetic energy, ε represents the dissipation rate, and the equations of the turbulent kinetic energy k and the dissipation rate ε are written as:

[0055]

[0056] In the above formula, G k = μ t S, η = Sk / ε, S is the modulus of the mean strain rate tensor S, σ k represents the turbulent Prandtl number of the turbulent kinetic energy k and σ k = 0.7194, σ ε represents the turbulent Prandtl number of the dissipation rate ε and σε = 0.7194. The constant C ε1 = 1.42, C ε2 = 1.68, η 0 = 4.38, β = 0.012.

[0057] IV. Based on the SIMPLE algorithm, the global implicit sub-relaxation method is combined with the local explicit sub-relaxation method to numerically calculate the discrete control algebraic equations of any control volume element in the computational domain, and the velocity field, pressure field, and flow rate at the surface of the control volume element are obtained. Please combine with Figure 1 the process schematic diagram, including the following steps:

[0058] Step 1: Initialize the flow field, set the initial velocity field, initial pressure field, and turbulent characteristic quantities. Among them, the advection characteristic quantities include turbulent kinetic energy k and dissipation rate ε.

[0059] Step 2: Determine the velocity guess value and pressure guess value

[0060] of the control volume element P at the nth iteration step of any tth time step. Among them, when t = n = 1, the initial velocity field set during initialization is used as the velocity guess value

[0061] and the initial pressure field is used as the pressure guess value When n ≥ 2, the calculation result of the velocity field at the (n - 1)th iteration step of the tth time step is used as the velocity guess value

[0062] and the calculation result of the pressure field at the (n - 1)th iteration step of the tth time step is used as the pressure guess value When n = 1 and t ≥ 2, the calculation result of the velocity field at the last iteration step of the (t - 1)th time step is used as the velocity guess value

[0063] Step 3: According to the velocity guess value and pressure guess value at the nth iteration step of the tth time step, calculate the coefficient matrix and source term vector of the velocity algebraic equations by combining the global implicit sub-relaxation method with the velocity algebraic equations (8) of the control volume element P, and then iteratively solve the velocity algebraic equations to obtain a new velocity guess value

[0064] As described above, when obtaining the velocity algebraic equation set (8) from the momentum discrete control equation (2) in the present application, the second-order upwind scheme with a gradient limiter is adopted for the convection term, which can effectively avoid the extreme value problem in the interface velocity reconstruction caused by uneven grid density and large inclination rate, improve the calculation stability, and enhance the calculation accuracy. For the diffusion term, a super-relaxation non-orthogonal correction method is adopted, which takes into account the influence of poor grid orthogonality on the calculation results and is more in line with the physical reality. Therefore, the results calculated according to the velocity algebraic equation set (8) are more accurate.

[0065] Furthermore, the global implicit sub-relaxation method is adopted to ensure that the coefficient matrix of the velocity algebraic equation set (8) is strictly diagonally dominant, control the variation degree of the flow field velocity between each iteration step, and avoid the divergence of the iteration process, thereby further improving the calculation stability. For any (n + 1)-th iteration step, for the velocity v calculated by using the velocity algebraic equation set (8) in this iteration step n+1 The formula for implementing the global implicit sub-relaxation method can be expressed as:

[0066] v = v n + λ v (v n+1 - v n ) (12)

[0067] where λ v represents the sub-relaxation factor of the velocity and 0 < λ v < 1, and its typical value is 0.7. v n is the velocity at the n-th iteration step, and v is the predicted velocity after relaxation of v n+1 .

[0068] Then, combining the global implicit sub-relaxation method, the velocity algebraic equation set (8) of the control volume cell P can be modified into the following form:

[0069]

[0070] According to the velocity guess value calculate the coefficient matrix and the source term vector in the modified velocity algebraic equation set (13), and then iteratively solve the velocity algebraic equation set (13) to obtain a new velocity guess value

[0071] Step 4, adopt the local explicit sub-relaxation method to update the new velocity guess value , including: for the velocity guess value in the region where it is greater than ξ times the oncoming flow velocity in the new velocity guess value use the explicit sub-relaxation technique to update it to where denotes the velocity at the previous iteration step, where λ is the relaxation coefficient and its typical value is 0.15. The parameter ξ > 1 and is generally taken as 2.5. The local explicit sub-relaxation method can limit the drastic change of velocity and further improve the computational stability.

[0072] Step 5: Use the updated velocity guess Combined with the pressure guess Update to obtain the flow rate at surface f of control volume cell P Including updating to obtain the flow rate at surface f of control volume cell P according to the Rhie-Chow momentum interpolation method It is:

[0073]

[0074] Step 6: According to the flow rate Calculate the coefficient matrix and source term vector in the pressure correction algebraic equation (9), and iteratively solve the pressure correction algebraic equation to obtain the pressure correction term p′ P .

[0075] Step 7: According to the pressure correction term p′ P Correct to obtain the velocity field, pressure field and flow rate at the surface of control volume cell P at the nth iteration step of the tth time step, including obtaining the calculation result v of the velocity field of control volume cell P at the nth iteration step of the tth time step P , the calculation result p of the pressure field P And the flow rate m at surface f f , and the calculation formula is:

[0076]

[0077] Where, is the predicted pressure value. p′ f Denotes the pressure correction term at surface f, and is obtained by interpolating the pressure correction terms of control volume cell P and control volume cell N on both sides of surface f.

[0078] Step 8: Judge whether the iteration process of the current time step converges according to the calculated residual and convergence criterion. When it is determined that the convergence state is not reached, use the obtained velocity field, pressure field and flow rate at the surface as the guess for the next iteration step and enter the calculation of the (n + 1)th iteration step of the tth time step, and repeat steps 2 to 8. When it is determined that the convergence state is reached, use the obtained velocity field, pressure field and flow rate at the surface as the guess for the first iteration step of the next time step, update the flow field flow time at the same time, and enter the calculation of the first iteration step of the (t + 1)th time step.

[0079] The above is only a preferred embodiment of the present application, and the present application is not limited to the above embodiments. It is understood that other improvements and changes directly derived or associated by those skilled in the art without departing from the spirit and concept of the present application should be considered to be included in the protection scope of the present application.

Claims

1. A CFD calculation method for ship stack resistance that is suitable for different grid qualities, characterized in that: The ship stacking resistance CFD calculation method includes: An incompressible RANS equation group including the continuity equation and the momentum equation is used as the control equation, the control volume unit is taken as the grid unit itself obtained by dividing the computational domain, and the control equation is discretized using the finite volume method to obtain the discrete control equation of each control volume unit; The second-order upwind scheme is used for the convection term in the discrete control equation of the control volume unit, and the super-relaxation non-orthogonal correction method is used for the diffusion term, and the discrete control algebraic equation group of the control volume unit is obtained; The turbulence model is used to close the discrete control algebraic equations, and based on the SIMPLE algorithm, the global implicit sub-relaxation method combined with the local explicit sub-relaxation method is used to perform numerical calculations on the discrete control algebraic equations of any control volume unit in the computational domain to obtain the velocity field, pressure field and flow rate at the surface of the control volume unit.

2. The CFD calculation method for ship stacking resistance according to claim 1 is characterized in that: The discrete control equations of the control volume unit P obtained by discretizing the control equation using the finite volume method include: For the momentum equation Discretization is performed to obtain the momentum discrete control equation of the control volume unit P: in, is an unsteady term, is the convection term, is the diffusion term, p f S f is the source term; v is the velocity vector of the flow field, p is the pressure, ρ is the fluid density, μ is the fluid viscosity, and μ t is the turbulent viscosity, nb(P) is the set of surfaces of the control volume unit P, f represents any surface of the control volume unit P, (ρvv) f represents the value of ρvv at the surface f of the control volume unit P, S f represents the area of ​​the surface f of the control volume unit P, p f represents the pressure at the surface f of the control volume element P, represents the control volume element P Value, V P represents the volume of the control unit P; And the continuity equation The continuous discrete control equation of the control volume unit P obtained by discretization is:

3. The CFD calculation method for ship stacking resistance according to claim 2 is characterized in that: The convection term in the momentum discrete control equation of the control volume unit P is expressed in the second-order upwind format as follows: Among them, the surface f of the control volume unit P is in the velocity field of the n+1th iteration step represents the velocity field of the upstream unit UP at the n+1th iteration step, It represents the velocity field of the upstream unit UP of the control volume unit P in the nth iteration step, represents the connection vector between the surface center of the surface f of the control volume element P and the center of the upstream element UP; represents the interface flux of the surface f of the control volume element P at the nth iteration step; The diffusion term in the momentum discrete control equation of the control volume unit P is expressed by the super-relaxation non-orthogonal correction method as follows: in, is the orthogonal linearization term, is a non-orthogonal nonlinear term; d PN is the distance between the center of the control volume unit P and the center of its adjacent control volume unit N, v P is the velocity at the center of the control volume element P, v N The velocity at the center of the control unit N adjacent to the control unit P; T f =S f -E f , where e is the unit vector connecting the control volume unit P and its adjacent control volume unit N.

4. The CFD calculation method for ship stacking resistance according to claim 3 is characterized in that: The discrete control algebraic equations of the control volume unit P obtained by sorting include: The first-order implicit format is used for the unsteady terms in the momentum discrete control equation of the control volume element P, and the pressure p in the source term in the momentum discrete control equation of the control volume element P is f The velocity algebraic equations of the control volume unit P are obtained by rearranging the momentum discrete control equations of the control volume unit P using the second-order central difference format. And, the pressure correction algebraic equation group of the control volume unit P is obtained by using the super-relaxation non-orthogonal correction method to sort out the diffusion discrete term in the continuous discrete control equation of the control volume unit P: in, The main diagonal coefficients of the coefficient matrix representing the velocity algebraic equations, represents the non-diagonal coefficients of the coefficient matrix of the velocity algebraic equations; NB(P) is the set of control units adjacent to the control unit P, and N represents any control unit adjacent to the control unit P; v P represents the velocity field of the control volume unit P, v N represents the velocity field of the control volume unit N, The source vector representing the velocity algebraic equations; The main diagonal coefficients of the coefficient matrix representing the pressure correction algebraic equation system, The non-diagonal coefficients of the coefficient matrix representing the pressure correction algebraic equation system; p′ P represents the pressure correction term of the control volume element P, p′ N represents the pressure correction term for the control volume unit N, represents the source term vector representing the pressure-corrected algebraic equation system.

5. The CFD calculation method for ship stacking resistance according to claim 4 is characterized in that: The unsteady term in the momentum discrete control equation of the control volume unit P is expressed in the first-order implicit format as follows: Among them, ρ P is the fluid density at the control volume element P, v n+1 is the velocity of the flow field at the n+1th iteration step, v n is the velocity of the flow field at the nth iteration step, Δt is the time interval between two adjacent iteration steps, V P represents the volume of the control unit P; The pressure p in the source term of the momentum discretization control equation for the control volume element P is f The second-order central difference format is expressed as: Among them, p P represents the pressure of the control volume unit P, p N represents the pressure of the control unit N adjacent to the control unit P, d Pf represents the connection vector between the surface center of the surface f of the control volume unit P and the center of the control volume unit P, d Nf represents the connection vector between the surface center of the surface f of the control volume unit P and the center of the control volume unit N; The velocity algebraic equations of the control volume unit P are sorted out 6. The CFD calculation method for ship stacking resistance according to claim 4 is characterized in that: Based on the SIMPLE algorithm, the global implicit sub-relaxation method combined with the local explicit sub-relaxation method is used to numerically calculate the discrete control algebraic equations of the control volume unit P, including: Initialize the flow field, set the initial velocity field, initial pressure field and turbulence characteristic quantities; Determine the velocity guess value v of the control unit P at the nth iteration step at any tth time step * P And the estimated pressure Among them, when t = n = 1, the initial velocity field is used as the velocity guess value Take the initial pressure field as the pressure guess value When n≥2, the velocity field calculation result of the n-1th iteration step of the tth time step is the velocity guess value The pressure field calculation result of the n-1th iteration step at the tth time step is taken as the pressure guess value When n=1 and t≥2, the velocity field calculation result of the last iteration step of the t-1th time step is used as the velocity guess value. The pressure field calculation result of the last iteration step of the t-1th time step is taken as the pressure guess value Guess the value based on the speed of the nth iteration step at the tth time step And the estimated pressure According to the velocity algebraic equations of the control volume unit P, the coefficient matrix and source term vector of the velocity algebraic equations are calculated by combining the global implicit sub-relaxation method, and then the velocity algebraic equations are iteratively solved to obtain a new velocity guess value. The new velocity guess value is obtained by using the local explicit sub-relaxation method. Make updates; Use the updated speed guess value Combined pressure guess Update the flow at the surface f of the control volume unit P According to the flow Calculate the coefficient matrix and source term vector in the pressure correction algebraic equations, and iteratively solve the pressure correction algebraic equations to obtain the pressure correction term p′ P ; According to the pressure correction term p′ P The velocity field, pressure field and flow rate on the surface of the control volume unit P at the nth iteration step of the tth time step are corrected; When it is determined that the convergence state has not been reached, the calculation of the n+1th iteration step of the tth time step is entered; when it is determined that the convergence state is reached, the calculation of the first iteration step of the t+1th time step is entered.

7. The CFD calculation method for ship stacking resistance according to claim 6 is characterized in that: The calculation based on the velocity algebraic equations of the control volume element P combined with the global implicit sub-relaxation method includes: Combined with the global implicit sub-relaxation method, the velocity algebraic equations of the control volume unit P are modified to: Guess value based on speed Calculate the coefficient matrix and source term vector in the modified velocity algebraic equations, and iteratively solve the pressure correction algebraic equations to obtain the new velocity guess value 8. The CFD calculation method for ship stacking resistance according to claim 6 is characterized in that: The new velocity guess value is obtained by using the local explicit sub-relaxation method. Updates include: New speed guess The estimated value of the velocity in the area greater than ξ times the incoming flow velocity Using the explicit sub-relaxation technique, it is updated to in, represents the speed of the previous iteration step, λ is the relaxation coefficient, and the parameter ξ>

1.

9. The CFD calculation method for ship stacking resistance according to claim 6 is characterized in that: Update the flow at the surface f of the control volume unit P include: According to the Rhie-Chow momentum interpolation method, the flow at the surface f of the control volume unit P is updated for: Among them, V N Represents the volume of the control volume unit N adjacent to the control volume unit P.

10. The CFD calculation method for ship stacking resistance according to claim 6, characterized in that: According to the pressure correction term p′ P The velocity field, pressure field and flow rate on the surface of the modified control volume element P at the nth iteration step of the tth time step include: Get the velocity field calculation result of the control volume unit P at the nth iteration step in the tth time step And get the pressure field calculation result of the control volume unit P at the nth iteration step of the tth time step in, is the predicted pressure value; And get the flow of the control volume unit P surface f at the nth iteration step at the tth time step Among them, p′ f represents the pressure correction term at the surface f, and is obtained by interpolating the pressure correction terms of the control volume element P and the control volume element N on both sides of the surface f.

Images