CN104392105A

CN104392105A - IB (Immersed Boundary) method

Info

Publication number: CN104392105A
Application number: CN201410602439.1A
Authority: CN
Inventors: 袁瑞峰; 钟诚文; 张贺
Original assignee: Northwestern Polytechnical University
Current assignee: Northwestern Polytechnical University
Priority date: 2014-10-27
Filing date: 2014-10-27
Publication date: 2015-03-04

Abstract

The invention discloses an IB (Immersed Boundary) method and belongs to a method for treating the boundary of a moving object with complex processing in computational fluid mechanics. According to the method, simple Cartesian grids are adopted; a series of Lagrangian points are used for expressing the boundary of the object; the influence of the boundary on fluid in a surrounding unit is expressed as the action of force; the force is solved by adopting an implicit iteration method, so that non-slip boundary conditions are guaranteed; a flow field is solved by using a BGK (Bhatnagar Gross Krook) format under the action of external force; a unit interface flux is considered under the influence of an immersed boundary force. According to the IB method, a good precision is obtained; non penetrability of the boundary can be better guaranteed; complexity and winding of the moving object in incompressible viscous flow can be conveniently simulated; any special processing on the grids is avoided.

Description

an IB approach

所属技术领域Technical field

本发明涉及一种计算流体力学中，在计算有粘不可压流动时，处理复杂、运动物体边界的方法。The invention relates to a method for dealing with complex and moving object boundaries when calculating viscous and incompressible flows in computational fluid dynamics.

背景技术Background technique

现代计算流体力学中，经常需要处理具有复杂外形、运动的物体的流动问题。对具有复杂曲线边界的物体，传统的处理方法是生成贴合物体边界的结构或者非结构网格，考虑到网格拓扑结构、网格质量等问题，要生成满意的网格需要较为复杂的算法，甚至需要有经验的人员反复调整，效率较低。近年来发展起来的基于笛卡尔网格的切割单元法，利用物体边界从背景网格切出计算网格，网格生成快捷，但是容易产生不规则、非凸、小体积网格单元等。对运动物体，传统的处理办法包括网格重长和网格运动。其中网格重长采用重新生成整个或者局部网格的方法来处理移动后的边界，其缺点是效率较低，且由于新生成网格与前一时刻网格不匹配还需要流场插值。网格运动包括动态嵌套网格和网格变形方法。对前者，在嵌套重叠区域插值是一个问题，而且在整个可能的运动区域通常需要较密网格导致网格数目较大。对网格变形，变形效率、变形后网格的质量和有效性较难控制。In modern computational fluid dynamics, it is often necessary to deal with the flow problems of objects with complex shapes and motions. For objects with complex curved boundaries, the traditional processing method is to generate a structured or unstructured grid that fits the boundaries of the object. Considering the grid topology and grid quality, a more complex algorithm is required to generate a satisfactory grid. , Even needs experienced personnel to adjust repeatedly, and the efficiency is low. The Cartesian grid-based cutting unit method developed in recent years uses the object boundary to cut out the calculation grid from the background grid. The grid generation is fast, but it is easy to generate irregular, non-convex, and small-volume grid units. For moving objects, the traditional processing methods include mesh re-length and mesh movement. Among them, the method of regenerating the whole or partial grid is used to deal with the moved boundary. The disadvantage is that the efficiency is low, and because the newly generated grid does not match the grid at the previous moment, flow field interpolation is required. Mesh motion includes dynamically nested meshes and mesh deformation methods. For the former, interpolation is a problem in nested overlapping regions, and generally requires a denser mesh resulting in a larger number of meshes throughout the region of possible motion. For mesh deformation, it is difficult to control the deformation efficiency, the quality and effectiveness of the deformed mesh.

近年来发展起来的IB方法(immersed boundary method，浸入边界法)，是处理复杂/运动边界的新兴热门方法。该方法在全流场使用笛卡尔网格，网格不需要贴合物体表面，网格生成极其简单。根据Mittal和Iaccarino于2005年在期刊《Annual Review of Fluid Mechanics》上发表的“Immersed boundary methods”一文，目前的IB方法大致可分为两类，一类是将边界对周围流体的力计算出来，均匀地分配到周围网格点，称为“连续力法”，该方法中力的计算不需要具体考察边界周围网格点的相对位置关系，更易应用于运动边界问题。另一类是直接插值边界周围网格点的物理量，称为“直接力法”，该方法需要较为精确地计算边界周围网格点的空间位置关系来进行插值，更适合于高雷诺数计算。目前的IB方法，总的说来较少应用于有限体积法，而且从没有一种IB方法计算了边界对周围有限体积单元的界面通量的影响，而这对保证边界的无穿透性是比较重要的。The IB method (immersed boundary method) developed in recent years is an emerging popular method for dealing with complex/moving boundaries. This method uses a Cartesian grid in the full flow field. The grid does not need to fit the surface of the object, and the grid generation is extremely simple. According to the article "Immersed boundary methods" published by Mittal and Iaccarino in the journal "Annual Review of Fluid Mechanics" in 2005, the current IB methods can be roughly divided into two categories, one is to calculate the force of the boundary on the surrounding fluid, Evenly distributed to the surrounding grid points, called the "continuous force method", the calculation of the force in this method does not need to specifically examine the relative position relationship of the grid points around the boundary, and it is easier to apply to the problem of moving boundaries. The other is to directly interpolate the physical quantities of the grid points around the boundary, which is called the "direct force method". This method needs to calculate the spatial position relationship of the grid points around the boundary more accurately for interpolation, and is more suitable for high Reynolds number calculations. The current IB methods are generally less applied to finite volume methods, and none of the IB methods has ever calculated the influence of the boundary on the interface flux of the surrounding finite volume elements, which is essential for ensuring the impenetrability of the boundary more important.

发明内容Contents of the invention

本发明的目的在于改进现有IB方法，将其应用于有限体积之BGK格式中，增强边界无穿透性，构造具有较高精度、能够方便处理有粘不可压流动中复杂边界/运动物体的方法。The purpose of the present invention is to improve the existing IB method, apply it to the BGK format of finite volume, enhance the non-penetration of the boundary, the structure has higher precision, and can conveniently handle complex boundaries/moving objects in viscous and incompressible flows method.

本发明用到的BGK格式，载于Tian等人于2007年在期刊《Journal of ComputationalPhysics》上发表的“A three-dimensional multidimensional gas-kinetic scheme for theNavier-Stokes equations under gravitational fields”一文，后文中提到的“Tian等人的BGK格式”均指这一格式。该格式为有限体积格式，同时考虑了外力的作用。其计算流程和传统有限体积格式一致，包括单元内物理量重构、单元界面上重构出的初始物理量随时间的演化以获得数值通量、通过界面通量计算单元内空间平均值的投影，一共3个步骤。其中单元内物理量重构采用简单的线性重构，而计算通量采用含有外力项的BGK通量求解器，能够考虑外力作用下界面通量的改变，而外力在整个格式的计算中直接体现为加速度的形式。具体的实施细节此处不再赘述。The BGK format used in the present invention is contained in the article "A three-dimensional multidimensional gas-kinetic scheme for the Navier-Stokes equations under gravitational fields" published by Tian et al. in the journal "Journal of ComputationalPhysics" in 2007. The "BGK format of Tian et al." mentioned here refers to this format. This scheme is a finite volume scheme, taking into account the effect of external forces. Its calculation process is consistent with the traditional finite volume format, including the reconstruction of the physical quantity in the unit, the evolution of the initial physical quantity reconstructed on the unit interface with time to obtain the numerical flux, and the projection of the spatial average value in the unit through the interface flux calculation. 3 steps. Among them, the physical quantity reconstruction in the unit adopts simple linear reconstruction, and the calculation flux adopts the BGK flux solver with external force term, which can consider the change of interface flux under the action of external force, and the external force is directly reflected in the calculation of the whole format as form of acceleration. The specific implementation details are not repeated here.

类似传统IB方法中的连续力法，在本发明中用一系列拉格朗日点表示物体边界，单元的中心则为欧拉点，边界对周围流体的影响通过其对周围流体的力来体现。方法的总体思路是在拉格朗日点上求得力，将其分布至周围欧拉点，利用Tian等人的BGK格式进行整个流场的计算更新，进而满足在拉格朗日点上的速度无滑移条件，也即在拉格朗日点上的速度需要等于物体表面该点实际运动速度。Similar to the continuous force method in the traditional IB method, in this invention, a series of Lagrangian points are used to represent the boundary of the object, and the center of the unit is the Euler point. The influence of the boundary on the surrounding fluid is reflected by its force on the surrounding fluid . The general idea of the method is to find the force at the Lagrangian point, distribute it to the surrounding Euler points, and use the BGK format of Tian et al. to update the calculation of the entire flow field, and then satisfy the velocity at the Lagrange point No slip condition, that is, the velocity at the Lagrangian point needs to be equal to the actual velocity of the point on the surface of the object.

本发明中，为了方便Tian等人的BGK格式的实施，力表现为加速度的形式。将拉格朗日点上的加速度分布至欧拉点的式子为In the present invention, in order to facilitate the implementation of the BGK scheme of Tian et al., the force is expressed in the form of acceleration. The formula for distributing the acceleration on the Lagrange point to the Euler point is

$\overset{&RightArrow; &Right Arrow;}{a a} (({\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j},, {t t}^{n no + + 11})) = = \underset{l l}{Σ Σ} {\overset{&RightArrow; &Right Arrow;}{a a}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})) {δ δ}_{h h} ((\frac{| | {\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j} - - {\overset{&RightArrow; &Right Arrow;}{x x}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})) | |}{h h})) {Δs Δs}_{l l},, - - - - - - ((11))$

这是个线性的关系式。式中，为拉格朗日点加速度、坐标，为欧拉点加速度、坐标，h为空间步长，tⁿ⁺¹为第n+1步时间，s_l为第l个拉格朗日点到起始点弧长，Δs_l为第l个拉格朗日点对应弧段长度。δ_h为一维离散Delta函数，可表示为，This is a linear relationship. In the formula, is the Lagrangian point acceleration and coordinates, is the acceleration and coordinates of the Euler point, h is the space step size, t ⁿ⁺¹ is the n+1th step time, s _l is the arc length from the lth Lagrangian point to the starting point, Δs _l is the lth Lagrangian point The Grangian points correspond to arc lengths. δ _h is a one-dimensional discrete Delta function, which can be expressed as,

${δ δ}_{h h} ((r r)) = = \{\begin{matrix} \frac{11}{44} ((11 + + cos cos ((\frac{π π | | r r | |}{22})))),, & | | r r | | \leq \leq 22 \\ 00,, & | | r r | | > > 22 \end{matrix} - - - - - - ((22))$

可见边界对周围单元的力的作用只限于2个单元长度的半径范围(也就是式(2)的非0范围)。若已知时间tⁿ⁺¹时所有欧拉点上的加速度(也就是所有有限体积单元中心的加速度)，则可将流场利用Tian等人的BGK格式从tⁿ更新至时间tⁿ⁺¹。而加速度需要达到的效果是使tⁿ⁺¹时速度无滑移边界条件得到满足It can be seen that the effect of the boundary on the force of the surrounding units is limited to the radius range of 2 unit lengths (that is, the non-zero range of formula (2)). If the acceleration at all Euler points at time t ⁿ⁺¹ (that is, the acceleration at the center of all finite volume elements) is known, the flow field can be updated from t ⁿ to time t ⁿ⁺¹ using the BGK format of Tian et al. . The effect that the acceleration needs to achieve is to satisfy the no-slip boundary condition of the velocity at t ⁿ⁺¹

${\overset{&RightArrow; &Right Arrow;}{u u}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})) = = {\overset{&RightArrow; &Right Arrow;}{U u}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})),, - - - - - - ((33))$

式中，为拉格朗日点速度，为实际的边界运动速度。拉格朗日点速度可以通过欧拉点(也就是单元中心)的速度向拉格朗日点插值得到，In the formula, is the Lagrangian point velocity, is the actual boundary motion speed. The Lagrange point velocity can be obtained by interpolating the velocity of the Euler point (that is, the cell center) to the Lagrange point,

${\overset{&RightArrow; &Right Arrow;}{u u}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})) = = \underset{i i,, j j}{Σ Σ} \overset{&RightArrow; &Right Arrow;}{u u} (({\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j},, {t t}^{n no + + 11})) {D D.}_{h h} (({\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j} - - {\overset{&RightArrow; &Right Arrow;}{x x}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})))) {h h}^{22},, - - - - - - ((44))$

这也是个线性的关系式。其中，D_h为离散二维Delta函数This is also a linear relationship. Among them, D _h is a discrete two-dimensional Delta function

${D D.}_{h h} (({\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j} - - {\overset{&RightArrow; &Right Arrow;}{x x}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})))) = = {δ δ}_{h h} ((\frac{{x x}_{i i,, j j} - - {x x}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11}))}{h h})) {δ δ}_{h h} ((\frac{{y the y}_{i i,, j j} - - {y the y}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11}))}{h h})),, - - - - - - ((55))$

在以上一系列计算中，拉格朗日点加速度分布到欧拉点、欧拉点速度插值到拉格朗日点都是线性关系，但是，由于本方法通量计算时使用的含外力项的BGK求解器的原因，欧拉点上加速度会影响其单元界面的通量，所以欧拉点加速度对欧拉点速度的影响极其复杂，导致拉格朗日点加速度对拉格朗日点速度的影响也极其复杂，很难写成一个解析表达式。所以，为了使无滑移边界条件式(3)得到满足，本方法在每一时间步先假定拉格朗日点加速度的初值之后再利用迭代法不断进行修正。从时间tⁿ到tⁿ⁺¹的整个时间步迭代步骤如下(下面的叙述中k代表第k步迭代)：In the above series of calculations, the Lagrangian point acceleration distribution to the Euler point, and the Euler point velocity interpolation to the Lagrange point are all linear relationships. The reason for the BGK solver is that the acceleration on the Euler point will affect the flux of its unit interface, so the influence of the Euler point acceleration on the Euler point velocity is extremely complicated, resulting in the Lagrangian point acceleration on the Lagrangian point velocity. The influence is also extremely complex, and it is difficult to write an analytical expression. Therefore, in order to satisfy the no-slip boundary condition (3), this method assumes the initial value of the Lagrangian point acceleration at each time step After that, iterative method is used to continuously make corrections. The entire time step iteration steps from time t ⁿ to t ⁿ⁺¹ are as follows (k in the following description represents the kth step iteration):

1.利用上一时间步的拉格朗日加速度最终迭代值作为本时间步迭代初值 ${\overset{&RightArrow;}{a}}_{L}^{(0)} (s_{l}, t^{n + 1}) .$ 1. Use the Lagrangian acceleration of the previous time step final iteration value As the initial value of this time step iteration ${\overset{&Right Arrow;}{a}}_{L}^{(0)} ({the s}_{l}, t^{no + 1}) .$

2.用第k步迭代的拉格朗日加速度值通过式(1)计算欧拉点加速度 ${\overset{&RightArrow;}{a}}^{(k)} ({\overset{&RightArrow;}{x}}_{i, j}, t^{n + 1}) .$ 2. Use the Lagrangian acceleration value of the k-th iteration Calculate the Euler point acceleration by formula (1) ${\overset{&Right Arrow;}{a}}^{(k)} ({\overset{&Right Arrow;}{x}}_{i, j}, t^{no + 1}) .$

3.用Tian等人的BGK格式将整个流场从tⁿ到更新tⁿ⁺¹(如果这不是第0步迭代，则只需更新边界周围受到加速度影响、不为0的单元，因为其余单元在第0步迭代中已经更新且边界于其上施加的加速度恒为0)。3. Use the BGK format of Tian et al. to update the entire flow field from t ⁿ to t ⁿ⁺¹ (if this is not the 0th iteration, you only need to update the boundary around the influence of acceleration, is not 0, because the rest of the cells have been updated in the 0th iteration and the acceleration applied by the boundary on them is always 0).

4.利用式(4)，从时间tⁿ⁺¹的欧拉点速度插值出拉格朗日点速度 4. Using formula (4), the Euler point velocity from time t ⁿ⁺¹ Interpolate Lagrangian point velocities

5.第4步计算的不一定满足无滑移边界条件式(3)，需要对拉格朗日加速度值进行修正。如前所述，由于单元上的加速度会影响其界面通量，所以拉格朗日点加速度对拉格朗日点速度的影响极其复杂。而如果忽略加速度对界面通量的影响，则加速度造成的单元的速度的改变量为加速度乘以时间步长度，而且加速度从拉格朗日点分布到欧拉点式(1)、速度从欧拉点插值到拉格朗日点式(4)均为线性关系。所以如果忽略加速度对界面通量的影响，结合式(1)与式(4)可以直接写出拉格朗日点上加速度与拉格朗日点速度改变量的关系式，再利用拉格朗日点速度和实际边界速度间差值，可以直接计算出拉格朗日点上加速度的修正值 $Δ {\overset{&RightArrow;}{a}}_{L}^{(k)} (s_{l}, t^{n + 1}) :$ 5. Calculated in step 4 does not necessarily satisfy the no-slip boundary condition (3), the Lagrangian acceleration value needs to be Make corrections. As mentioned earlier, the effect of Lagrangian point acceleration on Lagrangian point velocity is extremely complex because the acceleration on the element affects its interface flux. However, if the influence of acceleration on the interface flux is ignored, the velocity change of the unit caused by acceleration is the acceleration multiplied by the time step length, and the acceleration is distributed from the Lagrangian point to the Euler point (1), and the velocity is from the Euler point The interpolation from the pulling point to the Lagrangian point formula (4) is a linear relationship. Therefore, if the influence of acceleration on the interface flux is ignored, the relationship between the acceleration at the Lagrangian point and the velocity change at the Lagrangian point can be written directly by combining formula (1) and formula (4), and then using the Lagrangian Daily speed and the actual boundary speed The difference between them can directly calculate the acceleration on the Lagrangian point Correction value of $Δ {\overset{&Right Arrow;}{a}}_{L}^{(k)} ({the s}_{l}, t^{no + 1}) :$

$\begin{matrix} {\overset{&RightArrow; &Right Arrow;}{U u}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})) - - {\overset{&RightArrow; &Right Arrow;}{u u}}_{L L}^{((k k))} (({s the s}_{l l},, {t t}^{n no + + 11})) = = \\ \underset{i i,, j j}{Σ Σ} \underset{l l}{Σ Σ} Δ Δ {\overset{&RightArrow; &Right Arrow;}{a a}}_{L L}^{((k k))} (({s the s}_{l l},, {t t}^{n no + + 11})) Δt Δt {δ δ}_{h h} ((\frac{| | {\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j} - - {\overset{&RightArrow; &Right Arrow;}{x x}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})) | |}{h h})) {D D.}_{h h} (({\overset{&RightArrow; &Right Arrow;}{x x}}_{i i,, j j} - - {\overset{&RightArrow; &Right Arrow;}{x x}}_{L L} (({s the s}_{l l},, {t t}^{n no + + 11})))) Δ Δ {s the s}_{l l} {h h}^{22},, \end{matrix} - - - - - - ((66))$

注意由于这里使用了加速度不影响界面通量的假设，而该假设在实际利用Tian等人的BGK格式求解流场时是不适用的，经修正后的拉格朗日点加速度依然是不能保证完全满足无滑移边界条件式(3)的，所以需要进行这整个迭代程序。Note that the assumption that the acceleration does not affect the interface flux is used here, but this assumption is not applicable when the BGK scheme of Tian et al. is used to solve the flow field. Satisfy the no-slip boundary condition (3), so the whole iterative procedure needs to be carried out.

6.判定是否满足收敛判据6. Judgment Whether the convergence criterion is satisfied

$Δ Δ {\overset{&RightArrow; &Right Arrow;}{a a}}_{L L}^{((k k))} (({s the s}_{l l},, {t t}^{n no + + 11})) < < {u u}_{ref ref} \times \times c c,, - - - - - - ((77))$

式中u_ref为特定问题的参考速度，c为控制精度的参数，一般较保守地取10^-5即可。若式(7)得到满足，则迭代结束，为最终满足一定精度的拉格朗日加速度，得到的tⁿ⁺¹的流场为最终流场。反之，则计算出新的拉格朗日点加速度值：In the formula, u _ref is the reference speed of a specific problem, and c is a parameter of control accuracy, which is generally 10 ^-5 conservatively. If formula (7) is satisfied, the iteration ends, In order to finally meet the Lagrangian acceleration with a certain accuracy, the obtained flow field of t ⁿ⁺¹ is the final flow field. Otherwise, the new Lagrange point acceleration value is calculated:

${\overset{&RightArrow; &Right Arrow;}{a a}}_{L L}^{((k k + + 11))} (({s the s}_{l l},, {t t}^{n no + + 11})) = = {\overset{&RightArrow; &Right Arrow;}{a a}}_{L L}^{((k k))} (({s the s}_{l l},, {t t}^{n no + + 11})) + + ωΔ ωΔ {\overset{&RightArrow; &Right Arrow;}{a a}}_{L L}^{((k k))} (({s the s}_{l l},, {t t}^{n no + + 11})),, - - - - - - ((88))$

式中，ω为松弛系数一般取1.2。然后回到步骤2，进行第k+1步的迭代。In the formula, ω is the relaxation coefficient and generally takes 1.2. Then go back to step 2 and perform the iteration of step k+1.

在以上方法中，使用了迭代法隐式求解边界对周围单元的力，使无滑移边界条件得以很好地满足。In the above method, the iterative method is used to implicitly solve the force of the boundary on the surrounding elements, so that the no-slip boundary condition can be well satisfied.

总的说来，本发明的有益效果在于：Generally speaking, the beneficial effects of the present invention are:

相对于传统用于处理复杂边界问题的结构网格、非结构网格方法，本方法无需生成贴合物体边界的网格，只需用一系列点表示物体边界，而计算网格则是简单的笛卡尔网格，处理起来非常简便，处理方法统一。当处理运动边界时，与传统的网格重长与网格运动方法相比，本方法网格不需要作任何复杂的变动，只需表达物面边界的拉格朗日点移动即可，其余计算流程与计算静止物体时保持一致，非常简捷。Compared with the traditional structured grid and unstructured grid methods used to deal with complex boundary problems, this method does not need to generate a grid that fits the boundary of the object, but only needs to use a series of points to represent the boundary of the object, and the calculation grid is simple Cartesian grid is very easy to handle and the processing method is unified. When dealing with moving boundaries, compared with the traditional mesh length and mesh movement method, this method does not require any complicated changes to the mesh, only the Lagrangian point expressing the boundary of the object surface is moved, and the rest The calculation process is consistent with the calculation of stationary objects, which is very simple.

在具有以上优点的同时，计算精度方面，利用本发明计算了一系列算例并与其它大量实验、数值结果进行比对，证明本发明精度优良，丝毫不亚于、某些方面甚至超过别的边界处理方法。表1中展示了本发明计算的不同雷诺数下静止单圆柱绕流算例平均阻力系数与参考文献中实验、数值结果的对比，其中，表中引用的计算结果按从上到下顺序分别为：While having the above advantages, in terms of calculation accuracy, a series of calculation examples have been calculated by using the present invention and compared with a large number of other experiments and numerical results, which proves that the present invention has excellent accuracy, no less than, some aspects even surpass other Boundary handling methods. Table 1 shows the comparison of the average drag coefficient of the static single-cylinder flow calculation example calculated by the present invention with the experimental and numerical results in the references at different Reynolds numbers, wherein the calculation results quoted in the table are in order from top to bottom: :

1.Silva等人于2003年在期刊《Journal of Computational Physics》上发表的“Numericalsimulation of two-dimensional flows over a circular cylinder using the immersed boundarymethod”一文中的数值计算结果；1. The numerical calculation results in the article "Numerical simulation of two-dimensional flows over a circular cylinder using the immersed boundary method" published by Silva et al. in the journal "Journal of Computational Physics" in 2003;

2.Park等人于1998年在期刊《KSME International Journal》上发表的“Numerical solutionsof flow past a circular cylinder at Reynolds numbers up to 160”一文中的数值计算结果；2. The numerical calculation results in the article "Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160" published by Park et al. in the journal "KSME International Journal" in 1998;

3.Ye等人于1999年在期刊《Journal of Computational Physics》上发表的“An accurateCartesian grid method for viscous incompressible flows with complex immersedboundaries”一文中的数值计算结果；3. Numerical calculation results in the article "An accurate Cartesian grid method for viscous incompressible flows with complex immersive boundaries" published by Ye et al. in the journal "Journal of Computational Physics" in 1999;

4.Tritton于1959年在期刊《Journal of Fluid Mechanics》上发表的“Experiments on the flowpast a circular cylinder at low Reynolds numbers”一文中的实验结果。4. The experimental results in the article "Experiments on the flowpast a circular cylinder at low Reynolds numbers" published by Tritton in the journal "Journal of Fluid Mechanics" in 1959.

可见，本发明计算结果与参考文献吻合良好，所有数据均处于众多数值计算结果与实验结果之间(Tritton的结果为实验结果)。另外的一些关于静止、运动物体模拟的算例比对结果可见说明书附图部分。所有结果均显示本发明在模拟静止、运动物体绕流问题上具有良好的精度。It can be seen that the calculation results of the present invention are in good agreement with the references, and all data are between the numerical calculation results and the experimental results (Tritton's results are experimental results). Other comparison results of static and moving object simulation can be seen in the accompanying drawings of the specification. All the results show that the present invention has good precision in simulating the flow around static and moving objects.

雷诺数ReReynolds number Re 2020 4040 4646 4747 6060 8080 100100 本发明的计算结果Calculation result of the present invention 2.0712.071 1.5481.548 1.4651.465 1.4581.458 1.4191.419 1.3761.376 1.3521.352 Silva等人的结果The results of Silva et al. 2.042.04 1.541.54 -- 1.461.46 -- 1.401.40 1.391.39 Park等人的结果The results of Park et al. 2.012.01 1.511.51 -- -- 1.391.39 1.351.35 1.331.33 Ye等人的结果The results of Ye et al. 2.032.03 1.521.52 -- -- -- 1.371.37 -- Tritton的结果Tritton's results 2.1032.103 1.6051.605 1.5271.527 1.5161.516 1.3981.398 1.3161.316 1.2711.271

表1 不同雷诺数下静止单圆柱绕流算例平均阻力系数Table 1 The average drag coefficient of flow around a stationary single cylinder at different Reynolds numbers

另外，与已有的IB方法相比，本方法应用在有限体积之BGK格式中，计算通量时采用含有力项的BGK通量求解器，能够将边界的力的作用考虑在通量计算中，更好保证了边界的无穿透性，这在以前的IB方法中是从未尝试过的。关于计算通量时考虑边界的力的作用与不考虑力的作用对计算结果造成的影响，可见附图部分的计算结果对比。对比显示，若计算通量考虑力的作用，边界封闭性更好，流线更加贴合边界。In addition, compared with the existing IB method, this method is applied in the BGK format of finite volume, and the BGK flux solver with force term is used to calculate the flux, which can take the force of the boundary into consideration in the flux calculation , which better guarantees the impenetrability of the boundary, which has never been attempted in previous IB methods. Regarding the influence on the calculation results of considering the force of the boundary and not considering the force when calculating the flux, see the comparison of the calculation results in the attached figure. The comparison shows that if the force is considered in the calculation flux, the boundary closure is better and the streamline fits the boundary better.

附图说明Description of drawings

图1为在计算式(1)时，将拉格朗日点上加速度分布至欧拉点示意图。如图所示，物体边界由一系列拉格朗日点表示，欧拉点为单元中心，虚线圆为一维Delta函数非0影响范围，画斜线的单元即为受边界力影响的单元，加速度由拉格朗日点通过式(1)插值到画斜线的单元的中心。Fig. 1 is a schematic diagram of the distribution of the acceleration on the Lagrangian point to the Euler point when the formula (1) is calculated. As shown in the figure, the boundary of the object is represented by a series of Lagrange points, the Euler point is the center of the unit, the dotted circle is the non-zero influence range of the one-dimensional Delta function, and the unit with the oblique line is the unit affected by the boundary force. The acceleration is interpolated from the Lagrangian point to the center of the hatched cell by Equation (1).

图2为在计算式(4)时，将欧拉点速度往拉格朗日点上插值的示意图，虚框为二维离散Delta函数的非0影响范围，画斜线的单元为受其影响的单元，速度由画斜线的单元的中心通过式(4)插值到拉格朗日点上。Figure 2 is a schematic diagram of interpolating the velocity of the Euler point to the Lagrangian point when calculating formula (4). The dotted frame is the non-zero influence range of the two-dimensional discrete Delta function, and the units with oblique lines are affected by it , the velocity is interpolated to the Lagrangian point from the center of the obliquely drawn unit through formula (4).

图3中，(a)为雷诺数Re＝20，40时，本发明与参考文献计算的单圆柱绕流表面压强分布结果对比图，图中箭头代表Re增大方向，其中的参考文献为Fornberg于1980年发表在期刊《Journal of Fluid Mechanics》上的论文“A numerical study of steady viscous flow past a circularcylinder”。(b)为Re＝60，80，100时，本发明与参考文献计算的单圆柱绕流表面压强分布结果对比图，图中箭头代表Re增大方向，这里的参考文献为Park于1998年在期刊《KSME InternationalJournal》上发表的论文“Numerical solutions of flow past a circular cylinder at Reynolds numbersup to 160”。图中显示，各组结果吻合良好。Among Fig. 3, (a) is Reynolds number Re=20, when 40, the present invention and the single cylinder flow surface pressure distribution result comparison figure of reference calculation, arrow represents Re increasing direction among the figure, and reference document wherein is Fornberg The paper "A numerical study of steady viscous flow past a circular cylinder" published in the journal "Journal of Fluid Mechanics" in 1980. (b) is when Re=60,80,100, the present invention and the single cylinder surface pressure distribution result comparison figure of reference calculation, the arrow in the figure represents Re increasing direction, and the reference document here is Park in 1998 Paper "Numerical solutions of flow past a circular cylinder at Reynolds numbersup to 160" published in the journal "KSME International Journal". The figure shows that the results of each group are in good agreement.

图4为Re＝100，KC＝5条件下振动圆柱算例中，180°相位时，在圆柱振动平衡位置附近无量纲x坐标四个位置处剖面的无量纲速度沿纵轴y方向分布图。其中(a)为x方向无量纲速度分量分布，(b)为y方向无量纲速度分量分布。图中的曲线/空心符号/实心符号，分别对应本发明数值计算结果/参考文献的数值计算结果/参考文献实验结果。这里，参考文献为Dütsch等人于1998年在期刊《Journal of Fluid Mechanics》上发表的论文“Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenternumbers”。图中显示，各组结果吻合良好，部分数据(如处分布)本发明的计算结果甚至更接近实验结果。Figure 4 shows the dimensionless x-coordinate near the vibration equilibrium position of the cylinder in the case of a vibrating cylinder under the condition of Re=100 and KC=5, when the phase is 180° The dimensionless velocity profile along the vertical axis y-direction at the four positions. Where (a) is the dimensionless velocity component in the x direction distribution, (b) is the dimensionless velocity component in the y direction distributed. The curves/hollow symbols/solid symbols in the figure respectively correspond to the numerical calculation results of the present invention/the numerical calculation results of the references/the experimental results of the references. Here, the reference is the paper "Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenternumbers" published by Dütsch et al. in the journal "Journal of Fluid Mechanics" in 1998. The figure shows that the results of each group are in good agreement, and some data (such as place distribution) The calculated results of the present invention are even closer to the experimental results.

图5为Re＝100，KC＝5条件下振动圆柱算例中，在一个圆柱振动周期内，圆柱所受轴向力随时间变化图。图中，纵轴为无量纲轴向力横轴为无量纲时间t/T，其中T为圆柱振动周期。图中的三组结果分别为本发明计算结果、参考文献数值计算结果、Morison公式拟合结果。这里，参考文献为Dütsch等人于1998年在期刊《Journal of Fluid Mechanics》上发表的论文“Low-Reynolds-number flow around an oscillating circular cylinder at lowKeulegan-Carpenter numbers”。Morison公式是广泛应用于工程应用中，估计物体在振荡流体中所受轴向力的半经验公式。由数值计算结果拟合出Morison公式中两个待定参数后，轴向力随时间变化可计算出来。由图5可见，本发明计算结果与参考文数值计算的结果几乎重合，与Morison公式的拟合结果也吻合良好。Fig. 5 is a diagram of the variation of the axial force on the cylinder with time in a vibration cycle of the cylinder in the case of a vibrating cylinder under the conditions of Re = 100 and KC = 5. In the figure, the vertical axis is the dimensionless axial force The horizontal axis is the dimensionless time t/T, where T is the vibration period of the cylinder. The three groups of results in the figure are respectively the calculation results of the present invention, the numerical calculation results of reference documents, and the fitting results of Morison formula. Here, the reference is the paper "Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers" published by Dütsch et al. in the journal "Journal of Fluid Mechanics" in 1998. The Morison formula is a semi-empirical formula widely used in engineering applications to estimate the axial force on an object in an oscillating fluid. After fitting the two undetermined parameters in the Morison formula from the numerical calculation results, the change of the axial force with time can be calculated. It can be seen from Fig. 5 that the calculation result of the present invention almost coincides with the numerical calculation result of the reference text, and also coincides well with the fitting result of the Morison formula.

图6为计算通量时考虑边界的力的影响与忽略力的影响的计算结果对比，算例为Re＝40之圆柱绕流，(a)为通量计算考虑力的圆柱后缘边界附近流线，可见圆柱内部后缘形成涡，且边界附近流线完美贴合边界。(b)为通量计算忽略力的流线，虽然没有出现明显的流线穿透，圆柱内部后缘涡消散，边界附近流线贴合程度明显下降。这说明通量计算考虑力能够更好保证边界封闭性。Figure 6 is a comparison of the calculation results of considering the influence of boundary force and ignoring the influence of force when calculating the flux. The calculation example is the flow around a cylinder with Re=40. (a) is the flow near the boundary of the rear edge of the cylinder considering the force in the flux calculation. It can be seen that a vortex is formed at the rear edge of the cylinder, and the streamline near the boundary perfectly fits the boundary. (b) The streamline of the force is ignored for the flux calculation. Although there is no obvious streamline penetration, the trailing edge vortex inside the cylinder dissipates, and the degree of fit of the streamline near the boundary decreases significantly. This shows that the consideration of force in flux calculation can better guarantee the boundary closure.

具体实施方式Detailed ways

若给定某一特定静止/运动边界问题，本发明一种较优、较有代表性的具体实施方式为：If given a specific stationary/moving boundary problem, a better and more representative embodiment of the present invention is as follows:

1.生成整个计算域的笛卡尔网格，给定每个单元的初始物理量；1. Generate a Cartesian grid of the entire computational domain, given the initial physical quantity of each unit;

2.将物体边界表达为一系列已知坐标的离散点，点间距为1.4～2个网格长度；2. Express the object boundary as a series of discrete points with known coordinates, and the point spacing is 1.4 to 2 grid lengths;

3.设定本时间步tⁿ⁺¹-tⁿ的拉格朗日加速度初值，若为第一个时间步t¹-t⁰，可取任意值，若不为第一时间步，则取上一时间步拉格朗日加速度最终迭代值作为本时间步迭代初值 3. Set the Lagrangian acceleration of this time step t ⁿ⁺¹ -t ⁿ Initial value, if it is the first time step t ¹ -t ⁰ , it can take any value, if it is not the first time step, then take the Lagrangian acceleration of the previous time step final iteration value As the initial value of this time step iteration

4.用本时间步第k步迭代的拉格朗日加速度值通过式(1)计算欧拉点加速度 ${\overset{&RightArrow;}{a}}^{(k)} ({\overset{&RightArrow;}{x}}_{i, j}, t^{n + 1});$ 4. Use the Lagrangian acceleration value of the kth iteration of this time step Calculate the Euler point acceleration by formula (1) ${\overset{&Right Arrow;}{a}}^{(k)} ({\overset{&Right Arrow;}{x}}_{i, j}, t^{no + 1});$

5.用Tian等人的BGK格式将整个流场从tⁿ到更新tⁿ⁺¹(如果这不是第0步迭代，则只需更新边界周围受到加速度影响、不为0的单元，即图1所示的画斜线单元)；5. Use the BGK format of Tian et al. to update the entire flow field from t ⁿ to t ⁿ⁺¹ (if this is not the 0th iteration, you only need to update the boundary around the impact of acceleration, The unit that is not 0, that is, the oblique line unit shown in Figure 1);

6.利用式(4)，从时间tⁿ⁺¹的欧拉点速度插值出拉格朗日点速度 6. Using formula (4), the Euler point velocity from time t ⁿ⁺¹ Interpolate Lagrangian point velocities

7.利用式(6)计算拉格朗日点上加速度的修正值 7. Use formula (6) to calculate the acceleration on the Lagrangian point Correction value of

8.判定修正值是否满足收敛判据式(7)，若式(7)得到满足，则迭代结束，为最终满足一定精度的拉格朗日加速度，得到的tⁿ⁺¹的流场为本时间步最终流场，反之，利用式(8)计算新的拉格朗日点加速度值回到步骤4进行第k+1步迭代；8. Determine the correction value Whether the convergence criterion formula (7) is satisfied, if the formula (7) is satisfied, the iteration ends, In order to finally meet the Lagrangian acceleration with a certain accuracy, the obtained flow field at t ^{n+1 is} the final flow field at this time step, otherwise, use formula (8) to calculate the new Lagrangian point acceleration value Go back to step 4 for the k+1th iteration;

9.将边界每个点的位置从tⁿ更新到tⁿ⁺¹，回到步骤3进行下一时间步的计算。9. Update the position of each point on the boundary from t ⁿ to t ⁿ⁺¹ , and return to step 3 to calculate the next time step.

Claims

1. An IB method, a uniform Cartesian grid is used around the boundary, and a series of Lagrangian points are used to represent the boundary of the object. The center of the unit is an Euler point, and the effect of the boundary is controlled by the force of the Lagrange point on the surrounding fluid. It is characterized by: the method is applied in the BGK format, the force on the boundary will affect the interface flux of the unit, and the force on the Lagrangian point is solved implicitly by iteration to satisfy the infinite force on the Lagrangian point. Slip boundary conditions.

2. The IB method according to claim 1, characterized in that: the boundary force adopts the form of acceleration, the acceleration on the Lagrange point is distributed to the Euler point by a one-dimensional discrete Delta function, and the Euler The velocities at the points are interpolated to the Lagrangian points by a two-dimensional discrete Delta function.

3. according to the described IB method of claim 1 and 2, it is characterized in that, the described iterative method for solving the acceleration on the Lagrangian point is: the acceleration is distributed on the Euler point and acts on fluid, update After the flow field, the velocity is interpolated to the Lagrangian point and compared with the actual velocity of the object surface, and then the acceleration on the Lagrangian point is fed back to adjust, so that after the final evolution of each time step, the interpolation to the Lagrangian point is The velocity at the Grangian point is equal to the actual velocity of the object's surface.

4. the IB method according to claim 3 is characterized in that: after obtaining the acceleration on the Euler point, the described flow field update calculation adopts the BGK format containing the external force, and the acceleration of the Euler point passes through the BGK format containing the external force item. The BGK solver has an effect on element interface fluxes.

5. The IB method according to claim 3, characterized in that: in each Lagrangian point acceleration iteration of each time step, the described mode of feedback adjustment to the acceleration on the Lagrangian point Yes, assuming that the influence of the Euler point acceleration on the unit interface flux is ignored, the Lagrangian point acceleration required to make the Lagrangian point velocity equal to the actual moving velocity of the object surface can be directly algebraically solved under this assumption The correction amount, which is used as a feedback to adjust the Lagrangian point acceleration of each iteration.