CN104392105A - IB (Immersed Boundary) method - Google Patents

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

A kind of IB method

Art

The present invention relates in a kind of Fluid Mechanics Computation, calculating have sticky can not baric flow move time, process is complicated, the method on moving object border.

Background technology

In modern computational, process is often needed to have the flow field problem of object of complex appearance, motion.To the object with complex curve border, traditional disposal route is the structure or the unstrctured grid that generate laminating object boundary, consider the problem such as mesh topology, mesh quality, generate satisfied grid and need comparatively complicated algorithm, even need experienced personnel repeatedly to adjust, efficiency is lower.The cutter unit method based on cartesian grid that development in recent years is got up, utilizes object boundary to cut out computing grid from background grid, and mess generation is quick, but easily produces irregular, non-convex, small size grid cell etc.To moving object, traditional treating method comprises grid heavily length and mesh motion.Wherein grid heavily long adopt regenerate method that is whole or Local grid process mobile after border, its shortcoming is that efficiency is lower, and also needs flow field interpolation because newly-generated grid does not mate with previous moment grid.Mesh motion comprises moving chimera grids and grid deforming method.To the former, be a problem in nested overlapping region interpolation, and usually need comparatively fine grid to cause lattice number larger in whole possible moving region.To distortion of the mesh, the quality of deformation efficiency, the rear grid of distortion and the more difficult control of validity.

The IB method (immersed boundary method, immersed Boundary Method) that development in recent years is got up is the emerging popular method of process complexity/moving boundaries.The method uses cartesian grid at whole flow field, and grid does not need body surface of fitting, and mess generation is extremely simple.Civilian according to " the Immersed boundary methods " that Mittal and Iaccarino delivered on periodical " Annual Review of Fluid Mechanics " in 2005, current IB method is broadly divided into two classes, one class is calculated the power of border to surrounding fluid, be assigned to net point around equably, be called " continuous forcing approach ", in the method, the calculating of power does not need the relative position relation specifically investigating net point around border, is more easily applied to Problem of Moving Boundaries.Another kind of is the physical quantity of net point around direct interpolation border, is called " direct force method ", and the method needs the spatial relation of net point around comparatively accurately computation bound to carry out interpolation, is more suitable for high reynolds number and calculates.Current IB method, is generally speaking lessly applied to finite volume method, and never a kind of IB method calculates the impact of border on the around boundary flux of limited bulk unit, and this is important to what ensure border without penetrability.

Summary of the invention

The object of the invention is to improve existing IB method, be applied in the BGK form of limited bulk, strengthen border without penetrability, there is degree of precision, conveniently can process and sticky can not move the method for middle complex boundary/moving object by baric flow.

The BGK form that the present invention uses, be loaded in " A three-dimensional multidimensional gas-kinetic scheme for theNavier-Stokes equations under gravitational fields " that the people such as Tian delivered on periodical " Journal of ComputationalPhysics " in 2007 literary composition, the BGK form of people " Tian etc. " that mention hereinafter all refers to this form.This form is finite Volume Scheme, considers the effect of external force simultaneously.Its calculation process is consistent with conventional finite volume format, comprise the evolution in time of initial physical amount that the reconstruct of physical quantity in unit, unit interface reconstruct to obtain numerical flux, projection by spatial averaging in boundary flux computing unit, have 3 steps altogether.Wherein in unit, physical quantity reconstruct adopts simple linear reconstruction, and calculate the BGK flux solver of flux employing containing external force term, can consider the change of boundary flux under External Force Acting, and external force is directly presented as the form of acceleration in the calculating of whole form.Concrete implementation detail repeats no more herein.

Continuous forcing approach in similar conventional I B method, represents object boundary with a series of Lagrangian points in the present invention, and the center of unit is then Euler's point, and border is embodied the power of surrounding fluid by it the impact of surrounding fluid.The general thought of method tries to achieve power on Lagrangian points, be distributed to Euler's point around, the calculating utilizing the BGK form of the people such as Tian to carry out whole flow field upgrades, and then meeting the speed non-slip condition on Lagrangian points, the speed also namely on Lagrangian points needs to equal this actual motion speed of body surface.

In the present invention, the conveniently enforcement of the BGK form of the people such as Tian, power shows as the form of acceleration.By the formula that the acceleration profile on Lagrangian points is put to Euler be

$\stackrel{\→}{a}({\stackrel{\→}{x}}_{i,j},t^{n+1})=\underset{l}{\mathrm{\Σ}}{\stackrel{\→}{a}}_L(s_l,t^{n+1}){\mathrm{\δ}}_h\left(\frac{|{\stackrel{\→}{x}}_{i,j}-{\stackrel{\→}{x}}_L(s_l,t^{n+1})|}{h}\right){\mathrm{\Δs}}_l,---\left(1\right)$

This is a linear relational expression.In formula, for Lagrangian points acceleration, coordinate, for Euler puts acceleration, coordinate, h is spatial mesh size, t ⁿ⁺¹be the (n+1)th step time, s _lbe l Lagrangian points to starting point arc length, Δ s _lbe the corresponding segmental arc length of l Lagrangian points.δ _hfor one-dimensional discrete Delta function, can be expressed as,

${\mathrm{\δ}}_h\left(r\right)=\left\{\begin{array}{cc}\frac{1}{4}(1+\cos \left(\frac{\mathrm{\π}\left|r\right|}{2}\right)),& \left|r\right|\≤2\\ 0,& \left|r\right|>2\end{array}\right.---\left(2\right)$

The effect of visible borders to the power of surrounding cells is only limited to the radius (namely the non-zero scope of formula (2)) of 2 element lengths.If known time t ⁿ⁺¹acceleration (namely the acceleration of all limited bulk unit centers) on Shi Suoyou Euler point, then can utilize the BGK form of the people such as Tian from t by flow field ⁿbe updated to time t ⁿ⁺¹.And acceleration needs the effect reached to be make t ⁿ⁺¹hourly velocity is met without slip boundary condition

${\stackrel{\→}{u}}_L(s_l,t^{n+1})={\stackrel{\→}{U}}_L(s_l,t^{n+1}),---\left(3\right)$

In formula, for Lagrangian spot speed, for the Boundary motion speed of reality.Lagrange spot speed can be obtained to Lagrangian point interpolation by the speed of Euler's point (namely unit center),

${\stackrel{\→}{u}}_L(s_l,t^{n+1})=\underset{i,j}{\mathrm{\Σ}}\stackrel{\→}{u}({\stackrel{\→}{x}}_{i,j},t^{n+1})D_h({\stackrel{\→}{x}}_{i,j}-{\stackrel{\→}{x}}_L(s_l,t^{n+1}))h^2,---\left(4\right)$

This is also a linear relational expression.Wherein, D _hfor discrete two-dimensional Delta function

$D_h({\stackrel{\→}{x}}_{i,j}-{\stackrel{\→}{x}}_L(s_l,t^{n+1}))={\mathrm{\δ}}_h\left(\frac{x_{i,j}-x_L(s_l,t^{n+1})}{h}\right){\mathrm{\δ}}_h\left(\frac{y_{i,j}-y_L(s_l,t^{n+1})}{h}\right),---\left(5\right)$

In above series of computation, it is all linear relationship that Lagrangian points acceleration profile is interpolated into Lagrangian points to Euler's point, Euler's spot speed, but, due to the reason containing the BGK solver of external force term used when this method flux calculates, Euler puts the flux that acceleration can affect its unit interface, so Euler to put acceleration extremely complicated on the impact of Euler's spot speed, cause Lagrangian points acceleration also extremely complicated on the impact of Lagrangian spot speed, be difficult to be write as an analytical expression.So in order to make to be met without slip boundary condition formula (3), this method first supposes the initial value of Lagrangian points acceleration at each time step recycle process of iteration afterwards constantly to revise.From time t ⁿto t ⁿ⁺¹whole time step iterative step following (in describing below, k represents kth step iteration):

1. utilize the Lagrangian acceleration of a upper time step final iterative value as this time step iterative initial value ${\stackrel{\→}{a}}_L^{\left(0\right)}(s_l,t^{n+1}).$

2. with the Lagrangian accekeration of kth step iteration through type (1) calculates Euler and puts acceleration ${\stackrel{\→}{a}}^{\left(k\right)}({\stackrel{\→}{x}}_{i,j},t^{n+1}).$

3. with the BGK form of the people such as Tian by whole flow field from t ⁿto upgrading t ⁿ⁺¹(if this is not the 0th step iteration, then only need upgrade to be subject to around border acceleration action, be not the unit of 0 because remaining element has upgraded in the 0th step iteration and the acceleration perseverance that border applies thereon is 0).

4. utilize formula (4), from time t ⁿ⁺¹euler's spot speed interpolation goes out Lagrangian spot speed

5. the 4th step calculates not necessarily meet without slip boundary condition formula (3), need Lagrangian accekeration revise.As previously mentioned, because the acceleration on unit can affect its boundary flux, so Lagrangian points acceleration is extremely complicated on the impact of Lagrangian spot speed.And if ignore the impact of acceleration on boundary flux, the knots modification of the speed of unit that then acceleration causes is that acceleration is multiplied by time step length, and acceleration is distributed to Euler's point type (1) from Lagrangian points, speed is linear relationship from Euler's point interpolation to Lagrangian point type (4).If so ignore the impact of acceleration on boundary flux, convolution (1) and formula (4) directly can write out the relational expression of acceleration and Lagrangian spot speed knots modification on Lagrangian points, recycle Lagrangian spot speed with actual boundary speed between difference, can directly calculate acceleration on Lagrangian points modified value $\mathrm{\Δ}{\stackrel{\→}{a}}_L^{\left(k\right)}(s_l,t^{n+1}):$

$\begin{array}{c}{\stackrel{\→}{U}}_L(s_l,t^{n+1})-{\stackrel{\→}{u}}_L^{\left(k\right)}(s_l,t^{n+1})=\\ \underset{i,j}{\mathrm{\Σ}}\underset{l}{\mathrm{\Σ}}\mathrm{\Δ}{\stackrel{\→}{a}}_L^{\left(k\right)}(s_l,t^{n+1})\mathrm{\Δt}{\mathrm{\δ}}_h\left(\frac{|{\stackrel{\→}{x}}_{i,j}-{\stackrel{\→}{x}}_L(s_l,t^{n+1})|}{h}\right)D_h({\stackrel{\→}{x}}_{i,j}-{\stackrel{\→}{x}}_L(s_l,t^{n+1}))\mathrm{\Δ}{s}_l{h}^2,\end{array}---\left(6\right)$

Note the hypothesis not affecting boundary flux owing to being used herein acceleration, and this hypothesis is inapplicable when the BGK form of the people such as the actual Tian of utilization solves flow field, still can not ensure to meet completely without slip boundary condition formula (3) through revised Lagrangian points acceleration, so need to carry out this whole iterative program.

6. judge whether meet convergence criterion

$\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{a}}_L^{\left(k\right)}(s_l,t^{n+1})<u_{\mathrm{ref}}\&times;c,---\left(7\right)$

U in formula _reffor the reference velocity of particular problem, c is the parameter of control accuracy, generally more conservatively gets 10 ^-5.If formula (7) is met, then iteration terminates, for finally meeting the Lagrangian acceleration of certain precision, the t obtained ⁿ⁺¹flow field be final flow field.Otherwise, then calculate the Lagrangian points accekeration made new advances:

${\stackrel{\&RightArrow;}{a}}_L^{(k+1)}(s_l,t^{n+1})={\stackrel{\&RightArrow;}{a}}_L^{\left(k\right)}(s_l,t^{n+1})+\mathrm{\&omega;\&Delta;}{\stackrel{\&RightArrow;}{a}}_L^{\left(k\right)}(s_l,t^{n+1}),---\left(8\right)$

In formula, ω is that coefficient of relaxation generally gets 1.2.Then get back to step 2, carry out the iteration of kth+1 step.

In above method, employ process of iteration Implicit Method border to the power of surrounding cells, make to be met well without slip boundary condition.

Generally speaking, beneficial effect of the present invention is:

Relative to the structured grid, the unstrctured grid method that are conventionally used to process complex boundary condition, this method is without the need to generating the grid of laminating object boundary, and only need represent object boundary by series of points, computing grid is then simple cartesian grid, deal with very easy, disposal route is unified.When processing moving boundaries, heavily length is compared with mesh motion method with traditional grid, and this method grid does not need the variation making any complexity, and the Lagrangian points only need expressing object plane border move, all the other calculation process are consistent with when calculating stationary object, very simple and direct.

While there is above advantage, computational accuracy aspect, utilize the present invention to calculate series of numerical examples and compare with other great many of experiments, numerical result, proving that precision of the present invention is excellent, at all not second to, some aspect even exceedes other boundary processing method.Illustrate in table 1 static single cylindrical under the different Reynolds number that the present invention calculates stream in example average resistance coefficient and list of references test, the contrast of numerical result, wherein, the result of calculation quoted in table is by being sequentially respectively from top to bottom:

The numerical result of the people such as 1.Silva in " the Numericalsimulation of two-dimensional flows over a circular cylinder using the immersed boundarymethod " literary composition to deliver on periodical " Journal of Computational Physics " for 2003;

The numerical result of the people such as 2.Park in " the Numerical solutionsof flow past a circular cylinder at Reynolds numbers up to 160 " literary composition to deliver on periodical " KSME International Journal " for 1998;

The numerical result of the people such as 3.Ye in " the An accurateCartesian grid method for viscous incompressible flows with complex immersedboundaries " literary composition to deliver on periodical " Journal of Computational Physics " for 1999;

Experimental result in " Experiments on the flowpast a circular cylinder at low Reynolds numbers " literary composition that 4.Tritton delivered in nineteen fifty-nine on periodical " Journal of Fluid Mechanics ".

Visible, result of calculation of the present invention and list of references coincide good, and all data are all in (result of Tritton is experimental result) between numerous numerical result and experimental result.In addition some are about the visible Figure of description part of example comparison result that is static, moving object simulation.All results all show the present invention and have good precision on static, the moving object circumferential motion problems of simulation.

Reynolds number Re	20	40	46	47	60	80	100
								Result of calculation of the present invention	2.071	1.548	1.465	1.458	1.419	1.376	1.352
The result of the people such as Silva	2.04	1.54	-	1.46	-	1.40	1.39
								The result of the people such as Park	2.01	1.51	-	-	1.39	1.35	1.33
The result of the people such as Ye	2.03	1.52	-	-	-	1.37	-
								The result of Tritton	2.103	1.605	1.527	1.516	1.398	1.316	1.271

Under table 1 different Reynolds number, static single cylindrical streams example average resistance coefficient

In addition, compared with existing IB method, this method is applied in the BGK form of limited bulk, the BGK flux solver containing strong item is adopted when calculating flux, can the effect of the power on border be considered in flux calculates, better ensure that border without penetrability, be never attempted in this former IB method.Consider the effect of the power on border about when calculating flux and do not consider the impact that the effect of power causes result of calculation, the result of calculation contrast of visible accompanying drawing part.Contrast display, if calculate flux to consider the effect of power, closing of the frontier is better, and streamline is fitted border more.

Accompanying drawing explanation

Fig. 1 is when calculating formula (1), and acceleration profile on Lagrangian points is put schematic diagram to Euler.As shown in the figure, object boundary is represented by a series of Lagrangian points, Euler's point is unit center, broken circle is the non-zero coverage of one dimension Delta function, the unit drawing oblique line is the unit affected by boundary force, and acceleration is interpolated into the center of the unit of picture oblique line by Lagrangian points through type (1).

Fig. 2 is when calculating formula (4), by the schematic diagram of Euler's spot speed toward interpolation on Lagrangian points, empty frame is the non-zero coverage of two-dimensional discrete Delta function, the unit drawing oblique line is the unit affected by it, and speed is interpolated on Lagrangian points by the center through type (4) of the unit drawing oblique line.

In Fig. 3, a () is reynolds number Re=20, when 40, the single cylindrical that the present invention and list of references calculate streams surface pressure distribution results comparison diagram, in figure, arrow represents Re augment direction, and list of references is wherein Fornberg was published on periodical " Journal of Fluid Mechanics " paper " A numerical study of steady viscous flow past a circularcylinder " in 1980.B () is Re=60,80, when 100, the single cylindrical that the present invention and list of references calculate streams surface pressure distribution results comparison diagram, in figure, arrow represents Re augment direction, and list of references is here the paper " Numerical solutions of flow past a circular cylinder at Reynolds numbersup to 160 " that Park delivered on periodical " KSME InternationalJournal " in 1998.Show in figure, each group result is coincide good.

Fig. 4 is under Re=100, KC=5 condition in oscillating circular cylinder example, during 180 ° of phase places, and dimensionless x coordinate near cylinder vibration balancing position the nondimensional velocity of four position sections is along longitudinal y directions distribution plan.Wherein (a) is x direction nondimensional velocity component distribution, (b) is y direction nondimensional velocity component distribution.Curve/open symbols/filled symbols in figure, the respectively numerical result/list of references experimental result of corresponding numerical result/list of references of the present invention.Here, list of references is the paper " Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenternumbers " that the people such as D ü tsch delivered on periodical " Journal of Fluid Mechanics " in 1998.Show in figure, each group result is coincide good, partial data (as place distribution) result of calculation of the present invention is even closer to experimental result.

Fig. 5 is under Re=100, KC=5 condition in oscillating circular cylinder example, within a cylinder vibration period, and the variation diagram in time of axial force suffered by cylinder.In figure, the longitudinal axis is dimensionless axial force transverse axis is dimensionless time t/T, and wherein T is the cylinder vibration period.Three groups of results in figure are respectively result of calculation of the present invention, list of references numerical result, Morison formula fitting result.Here, list of references is the paper " Low-Reynolds-number flow around an oscillating circular cylinder at lowKeulegan-Carpenter numbers " that the people such as D ü tsch delivered on periodical " Journal of Fluid Mechanics " in 1998.Morison formula is widely used in engineer applied, estimates the semiempirical formula of object suffered axial force in oscillating fluid.Simulate in Morison formula after two undetermined parameters by numerical result, axial force changes in time can be calculated.As seen from Figure 5, result of calculation of the present invention almost overlaps with the result with reference to civilian numerical evaluation, also coincide good with the fitting result of Morison formula.

Fig. 6 is that the impact considering the power on border when calculating flux contrasts with the result of calculation of the impact ignoring power, example is the peripheral flow of Re=40, streamline near a cylinder trailing edge boundary that () is flux calculating consideration power, visible cylindrical inside trailing edge forms whirlpool, and streamline perfection laminating border near border.B () ignores the streamline of power for flux calculates, although do not occur that obvious streamline penetrates, dissipate in cylindrical inside trailing edge whirlpool, and near border, streamline laminating degree obviously declines.This illustrates that flux calculates consideration power and better can ensure closing of the frontier.

Embodiment

If given a certain specific static/Problem of Moving Boundaries, more excellent, the more representational embodiment of one of the present invention is:

1. generate the cartesian grid of whole computational fields, the initial physical amount of given each unit;

2. object boundary is expressed as the discrete point of a series of known coordinate, dot spacing is 1.4 ~ 2 Gridding lengths;

3. set this time step t ⁿ⁺¹-t ⁿlagrangian acceleration initial value, if first time step t ¹-t ⁰, desirable arbitrary value, if not very first time step, then gets a time step Lagrange acceleration final iterative value as this time step iterative initial value

4. with the Lagrangian accekeration of this time step kth step iteration through type (1) calculates Euler and puts acceleration ${\stackrel{\&RightArrow;}{a}}^{\left(k\right)}({\stackrel{\&RightArrow;}{x}}_{i,j},t^{n+1});$

5. with the BGK form of the people such as Tian by whole flow field from t ⁿto upgrading t ⁿ⁺¹(if this is not the 0th step iteration, then only need upgrade to be subject to around border acceleration action, be not the unit of 0, the picture oblique line unit namely shown in Fig. 1);

6. utilize formula (4), from time t ⁿ⁺¹euler's spot speed interpolation goes out Lagrangian spot speed

7. utilize formula (6) to calculate acceleration on Lagrangian points modified value

8. judge modified value whether meet convergence criterion formula (7), if formula (7) is met, then iteration terminates, for finally meeting the Lagrangian acceleration of certain precision, the t obtained ⁿ⁺¹flow field be the final flow field of this time step, otherwise, utilize formula (8) to calculate new Lagrangian points accekeration get back to step 4 and carry out kth+1 step iteration;

9. by the position of each for border point from t ⁿbe updated to t ⁿ⁺¹, get back to the calculating that step 3 carries out future time step.

Claims (5)

1. an IB method, even cartesian grid is adopted around border, object boundary is represented with a series of Lagrangian points, unit center is Euler's point, the effect on border is embodied the power of surrounding fluid by Lagrangian points, it is characterized in that: method is applied in BGK form, the acting force on border can have an impact to the boundary flux of unit, by the power on iteration Implicit Method Lagrangian points meet on Lagrangian points without slip boundary condition.

2. IB method according to claim 1, it is characterized in that: described border acting force adopts the form of acceleration, acceleration on Lagrangian points is distributed on Euler's point by one-dimensional discrete Delta function, and the speed on Euler's point passes through two-dimensional discrete Delta function interpolation on Lagrangian points.

3. IB method according to claim 1 and 2, it is characterized in that, described for solving the alternative manner of the acceleration on Lagrangian points is: degree of will speed up is distributed on Euler's point and acts on fluid, contrast on speed interpolation to Lagrangian points with body surface actual motion speed after upgrading flow field, and then feedback regulation is carried out to the acceleration on Lagrangian points, to make each time step after final evolution, the speed be interpolated on Lagrangian points equals body surface actual motion speed.

4. IB method according to claim 3, it is characterized in that: after trying to achieve the acceleration on Euler's point, described flow field upgrades the BGK form calculating and adopt containing external force, and the acceleration of Euler's point is had an impact to unit interface flux by the BGK solver containing external force term.

5. IB method according to claim 3, it is characterized in that: in each Lagrangian points acceleration iteration of each time step, described to the mode that the acceleration on Lagrangian points carries out feedback regulation be, suppose that ignoring Euler puts the impact of acceleration on unit interface flux, can directly algebraic solution be under this assumption, make Lagrangian spot speed and the correction of the equal required Lagrangian points acceleration of body surface actual motion speed, this correction is regulated the glug Lang point acceleration of each iteration as feedback.