CN104036126A - Numerical simulation method and system for particulate matter motion in contracting and expanding process of pulmonary alveoli - Google Patents

Abstract

The invention provides a numerical simulation method for particulate matter motion in the contracting and expanding process of pulmonary alveoli. The numerical simulation method comprises the following steps: establishing a geometrical model of the pulmonary alveoli according to preset geometrical parameters of the pulmonary alveoli and carrying out meshing on the geometrical model of the pulmonary alveoli to generate a dynamic mesh model of the pulmonary alveoli; acquiring particle data and a respiration parameter and a fluid parameter which correspond to the pulmonary alveoli; according to the dynamic mesh model, the geometrical parameters, the respiration parameter and the fluid parameter of the pulmonary alveoli and the particle data, simulating a flow field model in the pulmonary alveoli when the pulmonary alveoli contracts and expands; monitoring the motion process of particulate matter in the flow field model and acquiring motion data of the particulate matter. The invention also provides a corresponding system. The measurement period can be shortened; an operation load is reduced; the motion data of the particulate matter in the pulmonary alveoli is accurately acquired.

Description

Method for numerical simulation and the system of particle motion in alveolar contraction and expansion process

Technical field

The present invention relates to alveolar analog detection technical field, particularly relate to the method for numerical simulation of particle motion in a kind of alveolar contraction and expansion process, and the numerical simulation system of particle motion in the contraction of a kind of alveolar and expansion process.

Background technology

The major function of human body respiration be for each tissue of health provide oxygen and get rid of CO 2 waste gas, the respiratory of human body can be divided into two stages: by external environment to blood transport gas; Gas enters into each tissue via blood.Along with social, economic continuous progress, the mankind improve constantly the quality requirements of living environment, and the protection consciousness of living environment is strengthened gradually.The Particulate Pollution that commercial production and ecological deterioration bring has become one of important indicator evaluating quality of life and air quality.In atmospheric aerosol Particulate Pollution, the small aerosol particle of part wherein, especially inhalable particles is far-reaching especially on the mankind's healthy impact, entering after human respiratory tract, be not deposited on the conduction tracheae of respiratory tract, but being deep into the pulmonary deposition of the whole end of human respiratory tract, a lot of research shows the harm maximum of these particles to health.Research points out, mankind's numerous disease all pollutes and has direct or indirect contact with pellet.Alveolar is the main position of pulmonary gases exchange, and the kinetic characteristic of research pellet in alveolar, for helping pathogenesis and the gasoloid treatment of understanding inhalable particles to have very important meaning.Because the flow performance yardstick of human respiratory tract's alveolar region is the magnitude that is less than micron, gas phase in alveolar flows need to consider gas flow characteristic from miniature scale, the mode of traditional laboratory facilities monitoring alveolar endoparticle thing motion, its cycle is long, calculated amount is high and precision is low, is difficult to clearly describe the motion process of IA particle.

Summary of the invention

Based on this, the invention provides method for numerical simulation and the system of particle motion in a kind of alveolar contraction and expansion process, can shorten measuring period, reduce computational burden, obtain exactly the exercise data of IA particle.

In alveolar contraction and expansion process, a method for numerical simulation for particle motion, comprises the steps:

Set up the geometric model of alveolar according to the geometric parameter of default alveolar, the geometric model of described alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar;

Obtain particle data, and corresponding respiration parameter and the fluid parameter of described alveolar;

According to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, IA flow field model while simulating described alveolar contraction and expansion;

In described flow field model, monitor the motion process of described particle, obtain the exercise data of described particle.

A numerical simulation system for particle motion in alveolar contraction and expansion process, comprising:

Set up module, for set up the geometric model of alveolar according to the geometric parameter of default alveolar, the geometric model of described alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar;

Acquisition module, for obtaining particle data, and corresponding respiration parameter and the fluid parameter of described alveolar;

Analog module, for according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, simulates that described alveolar shrinks and IA flow field model when expansion;

Monitoring modular, for monitor the motion process of described particle at described flow field model, obtains the exercise data of described particle.

Method for numerical simulation and the system of particle motion in alveolar contraction of the present invention and expansion process, by alveolar is set up to geometric model, again the geometric model of alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar, in order to simulate the motion of alveolar in the time breathing; Then obtain particle data, geometric parameter, respiration parameter and fluid parameter that described alveolar is corresponding, according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter and fluid parameter, IA flow field model can simulate described alveolar contraction and expansion time, in described flow field model, monitor the motion process of described particle, obtain the exercise data of described particle, can shorten measuring period, reduce computational burden, obtain exactly IA Field Characteristics data.

Brief description of the drawings

Fig. 1 be alveolar of the present invention shrink and expansion process in the method for numerical simulation that moves of particle schematic flow sheet in one embodiment.

Fig. 2 is the method for numerical simulation geometric model schematic diagram of alveolar in one embodiment of particle motion in alveolar contraction of the present invention and expansion process.

Fig. 3 is the method for numerical simulation dynamic mesh model schematic diagram of alveolar in one embodiment of particle motion in alveolar contraction of the present invention and expansion process.

Fig. 4 be alveolar of the present invention shrink and expansion process in the method for numerical simulation that moves of the particle model schematic diagram of alveolar and tracheae in one embodiment.

Fig. 5 be alveolar of the present invention shrink and expansion process in the numerical simulation system that moves of particle structural representation in one embodiment.

Embodiment

Below in conjunction with embodiment and accompanying drawing, the present invention is described in further detail, but embodiments of the present invention are not limited to this.

As shown in Figure 1, be alveolar of the present invention shrink and expansion process in the method for numerical simulation that moves of particle schematic flow sheet in one embodiment, comprise the steps:

S11, set up the geometric model of alveolar according to default alveolar geometric parameter, the geometric model of described alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar;

S12, obtain particle data, and corresponding respiration parameter and the fluid parameter of described alveolar;

S13, according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, simulate that described alveolar shrinks and IA flow field model when expansion;

S14, in described flow field model, monitor the motion process of described particle, obtain the exercise data of described particle;

In the present embodiment, by alveolar is set up to geometric model, then the geometric model of alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar, in order to simulate the motion of alveolar in the time breathing; Then obtain geometric parameter, respiration parameter and fluid parameter that described alveolar is corresponding, according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter and fluid parameter, IA flow field model can simulate described alveolar contraction and expansion time, can measure described alveolar shrinking and Field Characteristics when expansion process according to the flow field model generating, can shorten measuring period, reduce computational burden, obtain exactly IA Field Characteristics data.

In the present embodiment, consider that natural flowing can be divided into Laminar Flow (laminar) and turbulent flow (turbulence).From the angle of experiment, in Laminar Flow, fluid does not have phase mutual interference between layers, between layer, does not have the transmission of quality there is no the transmission of momentum yet; And in turbulent flow process, fluid particle between layers has constantly mixed phenomenon of oozing mutually, the physical quantitys such as speed and pressure have the pulsating quantity of random nature on time and space.

Turbulent flow is that nature is ubiquitous mobile, and laminar flow belongs to other situation.Common laminar flow has flowing in kapillary or porous medium, flowing in flow-disturbing body surface boundary layer, and laminar flow generally appears in the characteristic length of the density of fluid, characteristic velocity and object less, or the situation such as the viscosity of fluid is very large.With Laminar Flow comparison, the research more complicated of turbulent flow, some basic problems solve not yet completely so far.

According to human respiratory tract's geometric shape, therein mobile of air can be considered as in pipe, flow, can define Reynold number (Re number) for flowing in pipe and be:

$\mathrm{Re}=\frac{\mathrm{ud}}{v}$

Wherein, the flow velocity that u is fluid, d is caliber, the dynamic viscosity that v is fluid.

It is generally acknowledged Re≤2300 o'clock, flowing for Laminar Flow in pipe; In the time of Re >=8000, flowing for turbulent flow in pipe; In the time of 2300<Re<8000, the zone of transition to turbulent flow in laminar flow flows.Calculate according to normal person's breathing state, gas is in the time that the upper respiratory tract flows, and because pharyngeal and throat's shape play the effect of jet pipe, local Re number may exceed 8000, flows and changes turbulence state into by laminar flow.And human respiratory tract alveolar region in the present embodiment, due to the continuous bifurcated of respiratory tract, the xsect of gas flow constantly increases, in respiratory tract end, flow velocity is generally cm/s, the diameter of tracheae is millimeter, and local Re number is far smaller than 100, belongs to typical Laminar Flow.

Be whether that constant can be divided into flow regimes and compressiblely flows and can not baric flow is moving flow according to the density of fluid in flow process.In the time that density is constant, fluid is incompressible fluid; When density is variable, fluid is compressible flowing.Although air is compressible fluid, in the laminar flow low speed of the present embodiment flows, can be used as incompressible fluid and process.

It is that whether temporal evolution is distinguished according to physical quantitys such as flow process speed medium velocity, pressure, temperature that stable-state flow and unstable state are flowed.When each physical quantity not when time to time change, flows for stable state Steady Flow in flow process; When physical quantity changes in time in flow process, for unstable state is flowed.

Because the flow performance yardstick of human respiratory tract's alveolar region is the magnitude that is less than micron, therefore need the impact of considering that small yardstick is analyzed gas flow.An important hypothesis studying flow field problem under macro-scale is exactly continuous medium hypothesis, in the time that the characteristic dimension of studied problem is far longer than the mean free path of fluid molecule, fluid can be considered as to a kind of non-individual body and treat, think that the molecular structure of fluid and molecular thermalmotion can only carry out by affecting the thermodynamic behaviour of fluid itself macroscopic motion of remote effect fluid.But along with the geometric scale of flow field problem constantly dwindles, while reaching micron or even nanometer scale, flow process just can not continue to be described with traditional macroscopic theory.

When different in research small size passage and macro-size passage between flow phenomenon, what first need to determine is the partition problem of channel size, generally the yardstick that is greater than 1mm is called to macro-scale, and the yardstick between 1um～1mm is called microscale.Also ununified standard of concrete partition of the scale, has plenty of according to physical dimension and divides, and has plenty of according to the scope of different acting forces in flowing and carries out partition of the scale.

In Micro-flows, some factors of usually ignoring in macroscopic view flows may be no longer negligible.Less characteristic dimension affects major embodiment in following aspect for mobile:

Under macroscopic conditions, the viscosity of fluid is only relevant with the character of fluid itself.And under the condition of microcosmic, the viscosity of fluid is subject to the impact of the many factors such as line size, shape of cross section, temperature, pressure, and can't use at present the mode of quantification to express exactly viscosity with the relation between these factors.In Micro-flows, just may not simply fluid viscosity be considered as to constant and calculate.

Along with reducing of geometric scale, in Micro-flows, surface force is main by plaing a part.Although fluid can not present polarity on the whole,, the fluid that contains polarity ion and non-polar fluid on flow performance, there is significant difference.Along with reducing of the yardstick that flows, roughness increases with respect to the ratio of the diameter of pipeline, and the impact that causes surfaceness is significantly increased.

Knudson number (Knudson number) is defined as:

$\mathrm{Kn}=\frac{\mathrm{\&lambda;}}{L}$

Wherein, the mean free path that λ is molecule, L is studied a question characteristic dimension.

Between the mean free path of molecule refers to that molecule collides for twice continuously, a free-moving mean path of molecule.The calculating formula of mean free path is:

$\mathrm{\&lambda;}=\frac{\mathrm{kT}}{\sqrt{2}{\mathrm{\&pi;d}}^{2}P}$

According to the size of Kn number in actual flow process, can be by gas the flow regimes division in pipeline as follows:

(1) Kn≤0.001 o'clock, flow process is observed continuous medium hypothesis, collision frequency between gas molecule is far above the collision frequency between gas molecule and body surface, in this case, newtonian viscous stress law and Fourier heat conduction law can be suitable for, and can solve the governing equation having without slip boundary condition and describe flow process;

(2) 0.001≤Kn≤0.1 o'clock, flow in slip region, collision between gas molecule is still far above the collision of gas molecule and body surface, but the intersection of gas and solid surface exists velocity-slip and temperature to jump, can solve the equation with slip boundary condition and describe flow process;

(3) 0.1≤Kn≤10 o'clock, flow in zone of transition, in this case, collision between collision between molecule and molecule and body surface is in same magnitude, need to consider to be positioned at flowing of this region very high to the requirement of computation model, DMSC (direct simulation Monte Carlo), the Molecular Dynamics method (MD) that can adopt the Boltzmann equation that comprises integration item or flow to carry out direct modeling simultaneously;

(4), when Kn > 10, flow for molecule Free-flow.Compare with the collision between gas molecule and body surface, intermolecular collision is negligible, can adopt the Boltzmann equation of collisionless integration item to describe flow process.

The effective diameter of air molecule is about d=3.1 × 10 ^-10m, can calculate in temperature according to the definition of kn number be 310K, pressure while being 101325Pa,, under the airflow condition in the present embodiment, molecule mean free path is approximately 8.7 × 10 ^-8m.

Be that 0.5mm calculates according to the mean diameter of the alveolar region main flow tracheae of human body, the passable Kn number to flow field is approximately Kn ≈ 1.74 × 10 ^-4.Therefore according to the division of Kn number, in the present embodiment, alveolar flow field need be set up the governing equation being set with without slip boundary condition and measures.

The gas exchange area geometric configuration of lower respiratory tract is very complicated, and as shown in Figure 2, geometric model is carried out the dynamic mesh model generating after grid is divided by Fig. 3 to the geometric model of the alveolar that the present embodiment is set up;

The foundation of model can be set up according to the geometric parameter of alveolar, and geometric parameter comprises centre distance and the alveolar angular aperture between alveolar tracheae pipe range, intracavity diameter, alveolar diameter, two alveolars; The geometric parameter of the present embodiment can be: pipe range (duct length) 740 μ m, intracavity diameter (lumen Diameter) 320 μ m, alveolar diameter (alveolar diameter) 320 μ m, centre distance between two alveolars is 340 μ m, 120 ° of alveolar angular apertures (opening angle), Re=ud/v=0.283.

In model, alveolar entrance is speed entrance, exports as pressure export border; Average alveolar flow when the inlet velocity of air-flow is pressed adult's eupnea is given, the parameters such as respiration parameter atmospheric density, respiratory cycle and the pressure coefficient under adult normal sitting state;

Wherein, when stable state, Re=0.283, the parameter of air is: ρ=1.225kg/m ³, Re=UD/ ν, U=0.0129m/s, μ=1.7894e-05kg/m.s, ν=μ/ρ=1.46e-05m ²/ s.When unstable state, IA unstable state air flow rate becomes sinusoidal curve to distribute in time, and respiratory cycle T is 4 seconds, is calculated as time step 0.02 second, and each time interval iterations is limited so that residual values trend is stable, $\stackrel{\&OverBar;}{u}=0.129m/s,$ Try to achieve $u_{\max }=\stackrel{\&OverBar;}{u}\&times;\mathrm{\&pi;}/2.$

In a preferred embodiment, described according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter and fluid parameter, while simulating described alveolar contraction and expansion, the step of IA flow field model comprises:

In the dynamic mesh model of described alveolar, simulate described IA flow field according to following formula:

$\left\{\begin{array}{c}\frac{{\&PartialD;u}_i}{{\&PartialD;x}_i}=0\\ \frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{{\&PartialD;x}_j}=-\frac{\&PartialD;P}{{\&PartialD;x}_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{{\&PartialD;x}_j}\\ \frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_{j}T\right)}{{\&PartialD;x}_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P{\&PartialD;x}_j}\right)}{{\&PartialD;x}_j}+S_T\end{array}\right.$

Wherein, ${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})-\frac{{2\mathrm{\&mu;\&delta;}}_{\mathrm{ij}}{\&PartialD;u}_i}{{3\&PartialD;x}_j};$

I, j=1,2,3; u _iand u _jfor fluid velocity, x _ifor flow field variable, t is the time, and ρ is fluid density, and P is the pressure on fluid infinitesimal, τ _ijfor viscous stress, c _pfor specific heat capacity, T is temperature, the heat transfer coefficient that k is fluid, S _tfor default viscous dissipation item, μ is kinetic viscosity.

Described according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter and fluid parameter, while simulating described alveolar contraction and expansion, the step of IA flow field model also can comprise:

In the dynamic mesh model of described alveolar, simulate contraction and the expansion process of described alveolar according to following formula:

$\begin{array}{c}\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\end{array}$

$\frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}$

$\frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l$

Wherein, $V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t},{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{{\mathrm{\&delta;V}}_l}{\mathrm{\&Delta;t}};$

ρ is fluid density, for fluid velocity, for the speed of motion grid node in described dynamic mesh model, Γ is coefficient of diffusion, for the border of control volume, for the area vector on border; V is the unit volume of motion grid, V ⁿfor the volume of motion grid described in the n moment, dV/dt controls the derivative of volume to the time, and n is the number of controlling volume face, S _φfor default broad sense source item, A _lbe the area vector of l face, δ V _lrepresent l face inswept volume in time Δ t, D is the diameter of alveolar place pipeline, and S is control area.

Described motion process of monitoring described particle in described flow field model, the step that obtains the exercise data of described particle comprises:

To motion control equation carry out the velocity of integration acquisition particle;

Wherein, the speed that u is air flow field, u _pfor the speed of particle, ρ _pfor the density of particle, the density that ρ is fluid, F _dfor acting on the power of dragging on particle, g _ifor gravitational vector, F _bthe Blang's power acting on particle;

$F_D=\frac{18\widetilde{\mathrm{\&mu;}}}{d_p^2{\mathrm{\&rho;}}_p{C}_C},C_C=1+\frac{2\mathrm{\&lambda;}}{d_p}[1.257+0.4\exp ({-1.1d}_p/2\mathrm{\&lambda;})],$ for the viscosity of fluid, C _cfor Cunningham slip correction coefficient, the mean free path that λ is molecule, d _pfor the diameter of particle;

ζ _ibe a Gaussian random vector with zero mean and unit variance, S ₀for the amplitude of white-noise process, the kinetic viscosity that v is air, k _b=1.38 × 10 ^-16for Boltzmann constant, the temperature that T is air flow field;

To governing equation carry out integration, obtain the position of particle in described flow field model described in each moment;

Wherein, u _pfor described velocity.

Any flow field problem all must meet law of conservation of mass, can obtain the mass-conservation equation (mass conservation equation) of fluid:

$\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{{\&PartialD;x}_i}=0,(i,j=\mathrm{1,2,3})$

In the present embodiment, fluid can be considered as incompressible fluid, and the density of air is considered as to constant, therefore above formula can be reduced to

$\frac{{\&PartialD;u}_i}{{\&PartialD;x}_i}=0$

The law of conservation of momentum and law of conservation of energy are also all necessary satisfied philosophys of any flow system, obtain energy conservation equation (momentum conservation equation) according to the law of conservation of momentum:

$\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{{\&PartialD;x}_j}=-\frac{\&PartialD;P}{{\&PartialD;x}_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{{\&PartialD;x}_j}$

In above formula, P is the pressure on fluid infinitesimal, τ _ijfor viscous stress, can be expressed as:

${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})-\frac{{2\mathrm{\&mu;\&delta;}}_{\mathrm{ij}}{\&PartialD;u}_i}{{3\&PartialD;x}_j};$

Can obtain energy conservation equation (energy conservation equation) according to law of conservation of energy:

$\frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_{j}T\right)}{{\&PartialD;x}_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P{\&PartialD;x}_j}\right)}{{\&PartialD;x}_j}+S_T$

Wherein cP is specific heat capacity, and T is temperature, the heat transfer coefficient that k is fluid, S _tfor viscous dissipation item.For moving by baric flow, if can not consider energy equation when heat exchange amount is very little, and in the present embodiment in order to calculate the effect of Blang's power, also need to solve energy equation.

Above-mentioned equation is the necessary satisfied Instantaneous Control equation of institute that flows, and can be:

$\frac{\&PartialD;\left(\mathrm{\&rho;\&phi;}\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_j\mathrm{\&phi;}\right)}{{\&PartialD;x}_j}=\frac{\&PartialD;\left(\frac{{\mathrm{\&Gamma;}}_{\mathrm{\&phi;}}\&PartialD;\mathrm{\&phi;}}{{\&PartialD;x}_j}\right)}{{\&PartialD;x}_j}+S_{\mathrm{\&phi;}}$

Wherein Γ _φfor generalized diffusion process coefficient, S _φfor broad sense source item, every time variation item, convective term, diffusion term and the source item of being followed successively by above formula.

For the partial differential equation of describing fluid motion, only the in the situation that of only a few, can obtain analytical solution.And cannot obtain the Fluid Control Equation of Accurate Analysis solution for those, can only use the method for numerical discretization, partial differential equation is discrete on time and space, be converted into algebraic equation and carry out computer solving, the precision of numerical method depends on the selection of discrete method.

Numerical simulation method for solving can comprise separation solving method and coupling solving method.In separation solving method, continuity equation, the equation of momentum and energy equation solve respectively; In coupling solving method, continuity equation, the equation of momentum and energy equation solve simultaneously, after above-mentioned governing equation is solved, then solve the equation such as turbulent flow, radiation.Solving liking of these two kinds of derivation algorithms is identical, be continuity equation, the equation of momentum and the energy equation that governing equation that they solve is describing mass conservation, momentum conservation and energy conservation, and they all use finite volume method as calculating object being carried out to the discrete basic methods solving.

Finite volume method (Finite Volume Method) is called again control volumetric method.It mainly comprises the steps:

(1) zoning is divided into a series of unduplicated discrete control volumes, and makes each net point have a control volume around;

(2) differential governing equation to be solved is controlled to volume integral to each, draw one group of discrete equation.

(3), by the linearize of discrete equation group, then obtain the iterative solution of variable by solving lienarized equation group.

Numerical simulation in the present embodiment is separation solving method.In the time solving low speed incompressible problem, the degree of coupling between speed, pressure, temperature is not very high, separates solving method and has very high numerical stability.The calculation procedure that separates solving method may be summarized to be:

(1) flow field variable update.In the time calculating for the first time, variable is upgraded by initialization procedure.In calculating subsequently, every iteration had once both obtained the solution of a renewal.

(2) with solving the equation of momentum when the value of fore pressure and mass flux, to obtain new velocity field.

(3) because the numerical solution of the velocity field obtaining in (2) cannot meet continuity equation completely, so solve again pressure update equation.Pressure update equation is the Poisson type equation of being derived by continuity equation, solves this equation and can obtain the correction to pressure field, velocity field and mass flux, and then continuity equation is met.

(4) utilize the solution of obtaining above, solve Equations of Turbulence, energy equation, constituent element equation and energy equation.

(5) if consider interphase interference in polyphasic flow calculates, need to calculate the source item solution in external phase equation by solving disperse phase track.

(6) whether the inspection condition of convergence is satisfied.If the condition of convergence is satisfied, stop calculating.If calculate not convergence, continue iterative process.

In separation derivation algorithm, because continuity equation and the equation of momentum do not need coupling, the flow field velocity value being calculated by the equation of momentum may not meet continuity equation simultaneously, one of solution is pressure to be revised by " Poisson class " equation obtaining of deriving of the equation of momentum after continuity equation and linearization with one, to be met pressure and the speed of continuity equation.Separate the pressure speed coupling process adopting in derivation algorithm, comprise three kinds of SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) method, SIMPLEC method and PISO (Pressure-Implicit with Splitting of Operators) methods.

The basic thought of SIMPLE is: first solve the equation of momentum with the pressure field of supposition and obtain the flux on frontier point, if the pressure field of supposing in calculating is inaccurate, the flux tried to achieve so just can not meet continuity equation.So add correction term on flux, make the flux of gained after revising can meet continuity equation.And add to such an extent that flux correction term is the function of pressure correction term, and therefore corrected flux is brought into continuity equation, just can obtain an equation about pressure correction term, obtain the solution of pressure correction term by solving this equation.Pressure correction term be multiplied by above sub-relaxation factor, then be added and just obtain a new pressure field with former pressure.Repeat said process taking this new pressure field as starting point, just form the iterative process that alternately solves pressure field, velocity field, to the last obtain convergence solution.

The main thought of PISO algorithm is that the double counting needing in the time solving pressure correction equation in SIMPLE and SIMPLEC algorithm is removed.After one or more extra PISO iteration step increasing, revised speed more closely meets continuity equation and the equation of momentum.In each iteration step, PISO algorithm increases a small amount of CPU time, but can reduce to a great extent iterations.PISO algorithm is highly suitable for the calculating of transient problem.

SIMPLEC algorithm is one of improvement algorithm of SIMPLE algorithm, and its algorithm is consistent with the basic ideas of SIMPLE algorithm, only on flux modification method, improves to some extent, and accelerates the speed of convergence of calculating.The stability of SIMPLEC form is better, can be by amplification suitable sub-relaxation factor in computation process.If do not use the subsidiary equations such as radiation model in laminar flow calculates, use SIMPLEC algorithm can accelerate to calculate speed of convergence.Because the mobile of the present embodiment research is typical Laminar Flow, therefore calculate air flow field with SIMPLEC algorithm, to accelerate to solve speed.

Lagrangian method and Euler's method are two kinds of basic skills describing fluid motion: under Lagrange is described, grid moves with fluid, and the speed of net point is identical with the speed of local fluid micellar; Under Euler describes, the locus of grid is fixed, and the speed of net point is zero, the control volume that fluid micellar forms through grid cell.These two kinds of basic skills are united in Arbitrary Lagrangian Euler's method (arbitary lagrangian-eulerian, ALE), allow grid with speed motion arbitrarily.In the time that Grid Velocity is zero, be Euler's method, be Lagrangian method in the time that the speed of mesh motion equals the speed of fluid.

For scalar φ arbitrarily, can obtain ALE form limited bulk governing equation:

$\begin{array}{c}\frac{d}{\mathrm{dt}}{\&Integral;}_V\mathrm{\&rho;\&phi;dV}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\end{array}$

The density that wherein ρ is fluid, for the speed of fluid, the speed of motion grid node, Γ is coefficient of diffusion, represent the border of control volume, represent the area vector on border.Above formula left end Section 1 time-derivative item is discrete with poor form after single order:

$\frac{d}{\mathrm{dt}}{\&Integral;}_V\mathrm{\&rho;\&phi;dV}=\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}$

The volume of grid cell is no longer constant but time dependent variable, and the volume Vn+1 in n+1 moment can be write as:

$V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t}$

In formula, dV/dt controls the derivative of volume to the time.

In the time of application motion grid numerical method, introduce extra error for fear of the motion of grid, how much conservation laws (GCL) were numerically met with the mass conservation, momentum conservation, the same needs of energy conservation.According to how much laws of conservation:

$\frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}$

$\frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{j}{\overset{n_j}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,j}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_j$

Wherein, n _jfor controlling the number of volume face, A _jit is the area vector of j face.Wherein

${\stackrel{\&RightArrow;}{u}}_{g,j}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_j=\frac{{\mathrm{\&delta;V}}_j}{\mathrm{\&Delta;t}}$

Wherein, δ V _jrepresent j face inswept volume in time Δ t.After grid upgrades, according to the variation of mesh volume, first solve the speed of grid, then solve the speed of fluid and other variable, therefore can in the dynamic mesh model of described alveolar, simulate described IA flow field.

The present embodiment is set up the dynamic mesh model of alveolar and is simulated, and grid is divided into structure and the large class of non-structured grid two in Fluid Mechanics Computation.

Structured grid is the Meshing Method generally adopting in the past.In structured grid, node is arranged in order, and the relation between adjacent node is clear and definite.In the conventional generation method of structured grid, there are algebraically generation method, conformal mapping method, Partial Differential Equation method, variational principle method etc.The advantage of structured grid is to process exactly the boundary condition in calculating, and computational accuracy is high, and can use a lot of implicit algorithm and multi-grid methods efficiently, and counting yield is high.But for the zoning with complex appearance, structured grid generates and is conventionally difficult to realize, even and if generate multi-block structured grid, it is also very complicated that boundary condition between piece and piece is processed, therefore the use of structured grid is very restricted.

In non-structured grid, grid cell and node each other unfixing rule can follow, and the distribution of node is arbitrarily completely.The hypothesis of its basic thought based on such: any area of space can be filled up by tetrahedron or triangular element.For the zoning with complicated shape, also can process very easily.Adopt the dividing mode of non-structured grid to be also convenient to adopt adaptive method refinement and alligatoring grid, make the division of grid more convenient and flexible, reduce the time of grid division in computation process, change violent, gradient flowing very greatly for parameter, be also convenient to catch the physical characteristics in flow field.Isostructure grid is compared, and non-structured grid also exists that calculated amount is large, to resolve device development be not also the shortcoming such as very ripe to high precision.

At present, in the time that pack processing contains the flow field problem of moving boundaries, generally process the distortion of grid based on non-structured grid.Conventional dynamic mesh deformation method has:

Spring analogy method (Spring Analog): spring analogy method solves aerofoil profile compulsive vibration flow-disturbing for the non-structured grid of triangle the earliest.The limit between any two grid nodes in system is idealized as a spring by elasticity smoothing method, and the initial mesh of system is regarded as a spring net case system that keeps balance.In system, the displacement of any one grid node all can cause being attached thereto in the spring connecing producing elastic force, and then causes the dynamic balance on proximity network lattice point to be broken.Involve away thus, in the time that the iteration through repeatedly to whole grid system reaches new balance, obtain a grid system after distortion.The spring analogy method of standard is not owing to considering the effect of shear stress, and deformability is not very strong, in two and three dimensions grid system, can not avoid crossing one another of Grid Edge.In the time there is larger moving displacement, grid cell serious distortion, it is non-negative situation that control volume easily occurs.For above-mentioned defect, can use improved spring analogy method, the deformability of grid is improved greatly and kept higher efficiency.

Spring analogy method can be divided into again seamed edge spring (segment spring) and two kinds of description forms of summit spring (vertex spring).In seamed edge spring is described, balance length equals the former length on limit, and net point is zero in suffered the making a concerted effort of original state.According to Hooke's law, after net point moves in spring system, act on making a concerted effort on certain grid node i wherein and be:

${\stackrel{\&RightArrow;}{F}}_i=\underset{j}{\overset{n_i}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ij}}(\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_j-\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_i)$

Wherein, with be respectively node i and and the vector shift of the adjacent node j of i.N _ifor the number of nodes being connected with node i, k _ijit is the elasticity coefficient of spring between node i, j.

In the time that spring system reaches equilibrium state, what each node was subject to above makes a concerted effort is zero.Carry out Jacobi iterative by the balance equation group to system node:

${\mathrm{\&Delta;}}_{x_i}^{\&RightArrow;m+1}=\frac{{\mathrm{\&Sigma;}}_i^{n_i}{k}_{\mathrm{ij}}\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_j^m}{{\mathrm{\&Sigma;}}_j^{n_i}{k}_{\mathrm{ij}}}$

Wherein, on border the displacement of node known, after borderline node updates, can carry out iterative to above formula.In above formula, be the shift value of node j in the time of the m time iteration, calculate and restrain the rear each node post exercise displacement vector in grid system that just can obtain.The position of node i after grid upgrades is:

${\stackrel{\&RightArrow;}{x}}_i^{n+1}={\stackrel{\&RightArrow;}{x}}_i^n+{\mathrm{\&Delta;}\stackrel{\&RightArrow;}{x}}_i^{\mathrm{converged}}$

Different with seamed edge spring describing method, the balance length of describing medi-spring at summit spring is zero.In fact, the description of these two kinds of spring analogy methods is consistent, and its concrete proof procedure can be with reference to other document.In order to make grid still can keep rational density to distribute after distortion, the coefficient of stiffiness of spring can have multiple choices, and wherein most widely used method is that certain power of the length of side between the coefficient of stiffiness of hypothesis spring and two nodes is inversely proportional to.

Elastomer method (Elasticity Method): in elastomer method, whole net region is seen as an elastic body, thereby the movement of all net points is by one group of fundamental equation control of Elasticity, by solving this prescription journey, just can obtain the displacement of grid system internal point.Elastomer method can be suitable in any trellis-type and type of sports, and there is stronger deformability, be more suitable in the situation that causes grid rotation due to object rotation, but the amount of calculation of elastomer method is large, counting yield is not as spring analogy method.

Local grid reconstructing method (Local Remeshing): Local grid reconstructing method develops at first in the input numerical simulation of two-dimentional store Combinations.In the problem of large displacement, in the time that larger distortion occurs grid, can, to the regenerating of Local grid, then try to achieve the flow field property value of newly-generated grid by the method for interpolation.Local grid reconstructing method generally comprises and calls automatic mesh generator and carry out interpolation two parts by old grid to new grid.Local grid reconstructing method has following shortcoming: the one, and the reconstructing method of Local grid need to regenerate grid in each time step, has increased widely the burden of calculating; The 2nd, the restructuring procedure of grid is returned numerical evaluation and has been brought interpolation error.In the non-structure dynamic mesh of three-dimensional problem, Local grid reconstruct is a very hard work.

Can adopt the control method of three kinds of dynamic meshes: elasticity theory of adjustment (Spring-based Smoothing), Local grid are heavily drawn method (Local Remeshing) and dynamic layer method (Dynamic Layering).

Elasticity theory of adjustment is described elasticity approximation method.In seamed edge spring describing method, the elasticity coefficient value of seamed edge is:

$k_{\mathrm{ij}}=\frac{1}{\sqrt{{\stackrel{\&RightArrow;}{x}}_j-{\stackrel{\&RightArrow;}{x}}_i}}$

The Local grid heavily method of drawing is supplementing elasticity smoothing method.In the time that grid system forms with triangle or tetrahedral grid, if the movement and distortion on border is far longer than size of mesh opening, may cause Local grid that serious distortion occurs, even occur that volume is negative situation, or the distortion of grid ambassador calculate and cannot restrain.Grid excessive aberration rate or that change in size is too violent can be concentrated in together and carries out repartitioning of Local grid, if the grid after repartitioning can meet the requirement of aberration rate and size, replace original grid with new grid, if new grid still cannot meet the demands, abandon the result of repartitioning.The Local grid heavily method of drawing can be applied to two-dimentional triangle and three-dimensional tetrahedral grid unit.

Dynamic layering method is the technology that dynamically increases or reduce clathrum on border according to the amount of movement on border, therefore dynamic layering method be applicable to hexahedral mesh, wedge shape grid etc. can be at the grid system of border higher slice.Dynamic layer technology is supposed the clathrum height value of an optimization on border, when being moved, being out of shape on border, if close on the height of one deck grid on border with optimizing height value while comparing greatly to a certain extent, just increase one deck grid between boundary surface and adjacent net compartment.Otherwise, if when closing on the grid on border and being compressed to a certain degree, closing on one deck grid again can be deleted, make borderline grid remain on certain density by this way.

The method for numerical simulation of Multiphase Flow can be divided into two kinds of Euler-Lagrangian method and Euler-Euler's methods.In Euler-Lagrangian method, gas phase is considered as external phase, by time equal Navier-Stokes equation solve.Follow the tracks of the motion of particle, bubble or the drop of some in flow field by the external phase flow field calculating, when calculating, can consider that Particle Phase is with the momentum occurring between gas phase, quality and energy exchange.This model is applied in the corresponding enough sparse situations of discrete particle in airflow field, and discrete phase only occupies very little volume fraction, and the volume fraction of the effect between particle and particle, particle can be ignored to the effect of external phase like this.

In Euler-Euler's method, different being regarded as mutually can interpenetrative external phase, introduce therein the concept of phase volume fraction, the volume fraction sum of each phase is 1, the volume fraction of each phase is a function continuous on time and space, according to conservation equation, each being on good terms obtained having the equation of analog structure, can make system of equations sealing solve by introducing the information of other experience.

In the present embodiment, particle only occupies volume fraction very little in diphasic flow, adopts Euler-Lagrangian method to carry out the motion of count particles phase.

Stressed in Gas-solid Two-phase Flow process of particle is quite complicated, and particle is relevant to the interaction between fluid, particle and particle, particle and wall.Generally speaking,, according to the mechanism of action between particle and fluid, two alternate acting forces comprise air resistance, Basset power, Saffman lift, Magnus lift, false mass force, pressure gradient power, thermophoretic forces, Bu Langli.Wherein, Basset power is that particle acceleration or deceleration in fluid moves and the power of generation, Saffman power is due to the power that exists velocity gradient to cause in flow field, Magnus power is that particle rotates the power producing in flow field, and Blang's power is the power that molecule is subject to the irregular collision effect generation of fluid molecule in flow field.

Although stressed very complicated in diphasic flow of particle, not all power all plays a part of equal importance.For the lower respiratory tract alveolar region of studying herein, the flow velocity of gas is very little, neglects the suffered Basset power of particle, Saffman lift, Magnus lift, false mass force in the calculating of the present embodiment.For the particle that can enter respiratory tract alveolar region, its temperature has approached the temperature of tracheae or alveolar wall, also the effect of the power such as thermophoretic forces, electrophoretic force can be ignored.Molecule in Laminar Flow, the resistance that Blang's power that it is subject to is subject to gravity and particle has identical magnitude, can not ignore.

In table 2.1, listed under the atmospheric pressure of standard, temperature is in the still air of 26 degree, the particle of different-grain diameter is the displacement in 1s under Blang's diffusion and gravity settling effect, can find out at the particle that is greater than 1 micron for particle diameter, the effect of Blang's power is negligible with respect to gravity, and for the particle that is less than 0.5um, the effect of Blang's power is greater than gravity.

The relative effect to movement of particles of table 2.1 Blang power and gravity

Under Lagrangian framework, the equation of motion of particle, by the stress balance in count particles, can calculate by following formula:

$\frac{{\mathrm{du}}_p}{\mathrm{dt}}=F_D(u-u_P)+\frac{g_i({\mathrm{\&rho;}}_p-\mathrm{\&rho;})}{{\mathrm{\&rho;}}_p}+F_B$

Wherein, the speed that u is air flow field, u _pfor the speed of particle, ρ _pfor the density of particle, the density that ρ is fluid, F _dfor acting on the power of dragging on particle, g _ifor gravitational vector, F _bthe Blang's power acting on particle.

For the molecule of studying for the present embodiment, intermolecular phorogenesis may play a big part.Therefore, the present embodiment calculates resistance F wherein by following formula _d:

$F_D=\frac{18\mathrm{\&mu;}}{d_p^2{\mathrm{\&rho;}}_p{C}_C}$

$C_C=1+\frac{2\mathrm{\&lambda;}}{d_p}[1.257+0.4\exp ({-1.1d}_p/2\mathrm{\&lambda;})]$

μ is the viscosity of fluid, C _cfor Cunningham slip correction coefficient, the mean free path that λ is molecule, d _pfor the diameter of particle.

Blang's power F in movement of particles equation _b, can regard a random vector as, solve time step with the physical property of particle and fluid and calculating and define:

$F_B={\mathrm{\&zeta;}}_i\sqrt{\frac{\mathrm{\&pi;}{S}_0}{\mathrm{\&Delta;t}}}$

$S_0=\frac{216v{k}_{B}T}{{\mathrm{\&pi;}}^2\mathrm{\&rho;}{d}_p^5{\left(\frac{{\mathrm{\&rho;}}_p}{\mathrm{\&rho;}}\right)}^2{C}_C}$

In above formula, ζ _ibe a Gaussian random vector with zero mean and unit variance, S ₀for the amplitude of white-noise process, the kinetic viscosity that v is air, k _b=1.38 × 10 ^-16for Boltzmann constant, the temperature that T is air flow field.

Obtain the velocity up of particle by the equation of motion integration to particle, try to achieve according to the integration of following formula the position in zoning of particle:

$\frac{{\mathrm{dX}}_p}{\mathrm{dt}}=u_P$

Above formula is carried out to integration, just can obtain each position of moment particle in flow field.

In simulation process, need the air-flow of a complete cycle of calculating in stable state and the transient of IA mobile 3D, the calculating of particle to be fled from alveolar completely or is deposited as end with particle, the maximum analog time is T.The computing method of the deposition (DE) of particle are: DE=Nd/Nt × 100%.Nd represents the precipitation number of particle at wall, and Nt represents total particle number.

As shown in Figure 4, be the model schematic diagram of the present embodiment alveolar and tracheae, the tracheae before and after using in the present embodiment in the cylindrical tube simulated respiration road of opening, describes the alveolar in respiratory tract with being connected to spherical cap on cylinder wall.R in the present embodiment _aand R _dbe respectively the radius of alveolar and tracheae, the angular aperture that γ is alveolar.

The radius R of tracheae in the present embodiment _dcan be 250um, the radius of alveolar is 200um, and the angular aperture of alveolar on main tracheae is made as 60 degree.The length of tracheae is 1000um, and in the computation model of the present embodiment, alveolar is positioned at the center position of this section of tracheae.

Gentle human body alveolar pipe can be considered as to an enlargement and contraction motion that maintenance geometric configuration is similar at the deformation process in respiratory.Each characteristic dimension of alveolar and tracheae is time dependent sine function:

$L\left(t\right)=L_0[1+\frac{\mathrm{\&beta;}}{2}+\frac{\mathrm{\&beta;}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\&pi;}}{2})]=L_0\mathrm{\&lambda;}$

Wherein L ₀for the value of while finishing (when air-breathing beginning or exhale) when the moment t=0 of each characteristic dimension in geometric model, f=2 π/T is the frequency of breathing, and T is the cycle length of breathing, the expansion amplitude that β is each characteristic dimension.

The feature bulking factor (excursion factor) that can define alveolar is:

C＝(V _max-V _min)/V _min

Wherein V _maxand V _minfor the minimum and maximum volume of geometric model, respectively corresponding expiration in actual respiratory ended and air-breathing volume at the end.Can obtain thus the characteristic dimension flare factor of alveolar:

β＝(C+1) ^1/3-1

In t=T/2 moment (when air-breathing end or while exhaling beginning), the characteristic length in model reaches maximal value: L _max=L ₀(1+ β);

In respiratory, the gas velocity in local Re number and tracheae on every first-order model is the function of time:

Re(t)＝U _d(t)D _d(t)/υ

$U_d\left(t\right)=4Q_d\left(t\right)/\left(\mathrm{\&pi;}{D}_d^2\left(t\right)\right)$

Can characterize alveolar residing position on respiratory tract by a dimensionless characteristic parameter, be defined as the gas flow Q that passes in and out alveolar in respiratory _awith the gas flow Q that flows through alveolar tracheae of living in _dratio.

On every one-level tracheae of alveolar region, according to the incompressible hypothesis of gas, the gas flow of turnover alveolar is:

$Q_a\left(t\right)=\frac{{\mathrm{dV}}_a}{\mathrm{dt}}---(3-7)$

V _afor the volume of unit alveolar on tracheae.Can obtain according to geometric model:

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\&pi;\&beta;}{\mathrm{fR}}_a^2{\mathrm{\&lambda;}}^2\cos (\mathrm{ft}-\frac{\mathrm{\&pi;}}{2})---(3-8)$

The gas flow that flows through this grade of tracheae can expression formula below be tried to achieve:

$Q_d\left(t\right)=\frac{{\mathrm{dV}}_t}{\mathrm{dt}}---(3-9)$

V _tby alveolate volume sum on this section of retrotracheal all tracheaes of process and these tracheaes, Q _dit is the flow that this part spatial expansion and contraction change produce.

In the present embodiment, also can be by get different Q in each geometric model _a/ Q _dvalue, represents this model residing position on whole alveolar region tracheae tree.In the respiratory of human body, maximum Re number occurs in the first order of alveolar region, expiration or air-breathing peak value moment, calculates the value that shows this maximum Re greatly about 20 left and right, and the Re number of most applications current downflow is less than 1.Under such flow field condition, can ignore the impact of inflow point's velocity distribution on whole flow field.For the pipe inner laminar flow of low reynolds number, the even speed from cross section distributes and develops into parabolic type velocity distribution, the pipe range L needing _ecan be by obtaining according to following expression formula:

$\frac{L_e}{D}=0.06\mathrm{Re}$

Wherein, the diameter that D is pipeline.With respect to the tracheae length before alveolar in the model of the present embodiment, this segment length can be disregarded.Therefore, set in the present embodiment one in inlet surface equally distributed speed be more rational, arrive the endotracheal fully development of having flowed when alveolar place.

As shown in Figure 5, be a kind of alveolar of the present invention shrink and expansion process in the numerical simulation system that moves of particle structural representation in one embodiment, comprising:

Set up module 51, for set up the geometric model of alveolar according to the geometric parameter of default alveolar, the geometric model of described alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar;

Acquisition module 52, for obtaining particle data, and corresponding respiration parameter and the fluid parameter of described alveolar;

Analog module 53, for according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, simulates that described alveolar shrinks and IA flow field model when expansion;

Monitoring modular 54, for monitor the motion process of described particle at described flow field model, obtains the exercise data of described particle.

In a preferred embodiment, described analog module 53 also for:

In the dynamic mesh model of described alveolar, simulate described IA flow field according to following formula:

$\left\{\begin{array}{c}\frac{{\&PartialD;u}_i}{{\&PartialD;x}_i}=0\\ \frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{{\&PartialD;x}_j}=-\frac{\&PartialD;P}{{\&PartialD;x}_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{{\&PartialD;x}_j}\\ \frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_{j}T\right)}{{\&PartialD;x}_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P{\&PartialD;x}_j}\right)}{{\&PartialD;x}_j}+S_T\end{array}\right.$

Wherein, ${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})-\frac{{2\mathrm{\&mu;\&delta;}}_{\mathrm{ij}}{\&PartialD;u}_i}{{3\&PartialD;x}_j};$

I, j=1,2,3; u _iand u _jfor fluid velocity, x _ifor flow field variable, t is the time, and ρ is fluid density, and P is the pressure on fluid infinitesimal, τ _ijfor viscous stress, c _pfor specific heat capacity, T is temperature, the heat transfer coefficient that k is fluid, S _tfor default viscous dissipation item, μ is kinetic viscosity;

And/or for

In the dynamic mesh model of described alveolar, simulate contraction and the expansion process of described alveolar according to following formula:

$\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;\&Delta;\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}$

$\frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}$

$\frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l$

Wherein, $V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t},{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{{\mathrm{\&delta;V}}_l}{\mathrm{\&Delta;t}};$

ρ is fluid density, for fluid velocity, for the speed of motion grid node in described dynamic mesh model, Γ is coefficient of diffusion, for the border of control volume, for the area vector on border; V is the unit volume of motion grid, V ⁿfor the volume of motion grid described in the n moment, dV/dt controls the derivative of volume to the time, and n is the number of controlling volume face, S _φfor default broad sense source item, A _lbe the area vector of l face, δ V _lrepresent l face inswept volume in time Δ t, D is the diameter of alveolar place pipeline, and S is control area.

In a preferred embodiment, described monitoring modular 94 also for:

To motion control equation carry out the velocity of integration acquisition particle;

Wherein, the speed that u is air flow field, u _pfor the speed of particle, ρ _pfor the density of particle, the density that ρ is fluid, F _dfor acting on the power of dragging on particle, g _ifor gravitational vector, F _bthe Blang's power acting on particle;

$F_D=\frac{18\widetilde{\mathrm{\&mu;}}}{d_p^2{\mathrm{\&rho;}}_p{C}_C},C_C=1+\frac{2\mathrm{\&lambda;}}{d_p}[1.257+0.4\exp ({-1.1d}_p/2\mathrm{\&lambda;})],$ for the viscosity of fluid, C _cfor Cunningham slip correction coefficient, the mean free path that λ is molecule, d _pfor the diameter of particle;

ζ _ibe a Gaussian random vector with zero mean and unit variance, S ₀for the amplitude of white-noise process, the kinetic viscosity that v is air, k _b=1.38 × 10 ^-16for Boltzmann constant, the temperature that T is air flow field;

To governing equation carry out integration, obtain the position of particle in described flow field model described in each moment;

Wherein, u _pfor described velocity.

Method for numerical simulation and the system of particle motion in alveolar contraction of the present invention and expansion process, by alveolar is set up to geometric model, again the geometric model of alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar, in order to simulate the motion of alveolar in the time breathing; Then obtain particle data, geometric parameter, respiration parameter and fluid parameter that described alveolar is corresponding, according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter and fluid parameter, IA flow field model can simulate described alveolar contraction and expansion time, in described flow field model, monitor the motion process of described particle, obtain the exercise data of described particle, can shorten measuring period, reduce computational burden, obtain exactly IA Field Characteristics data.

The above embodiment has only expressed several embodiment of the present invention, and it describes comparatively concrete and detailed, but can not therefore be interpreted as the restriction to the scope of the claims of the present invention.It should be pointed out that for the person of ordinary skill of the art, without departing from the inventive concept of the premise, can also make some distortion and improvement, these all belong to protection scope of the present invention.Therefore, the protection domain of patent of the present invention should be as the criterion with claims.

Claims (10)

1. a method for numerical simulation for particle motion in alveolar contraction and expansion process, is characterized in that, comprises the steps:

Set up the geometric model of alveolar according to the geometric parameter of default alveolar, the geometric model of described alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar;

Obtain particle data, and corresponding respiration parameter and the fluid parameter of described alveolar;

According to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, IA flow field model while simulating described alveolar contraction and expansion;

In described flow field model, monitor the motion process of described particle, obtain the exercise data of described particle.

2. the method for numerical simulation of particle motion in alveolar contraction according to claim 1 and expansion process, it is characterized in that, the geometric parameter of described alveolar comprises centre distance and the alveolar angular aperture between alveolar tracheae pipe range, intracavity diameter, alveolar diameter, two alveolars.

3. the method for numerical simulation of particle motion in alveolar contraction according to claim 1 and expansion process, is characterized in that, described respiration parameter comprises atmospheric density, respiratory cycle and pressure coefficient.

4. the method for numerical simulation of particle motion in alveolar contraction according to claim 1 and expansion process, is characterized in that, fluid parameter comprises Reynolds number, fluid density, fluid velocity and the fluid viscosity of air-flow.

5. the method for numerical simulation of particle motion in alveolar contraction according to claim 1 and expansion process, it is characterized in that, described according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, while simulating described alveolar contraction and expansion, the step of IA flow field model comprises:

In the dynamic mesh model of described alveolar, simulate described IA flow field according to following formula:

$\left\{\begin{array}{c}\frac{{\&PartialD;u}_i}{{\&PartialD;x}_i}=0\\ \frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{{\&PartialD;x}_j}=-\frac{\&PartialD;P}{{\&PartialD;x}_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{{\&PartialD;x}_j}\\ \frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_{j}T\right)}{{\&PartialD;x}_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P{\&PartialD;x}_j}\right)}{{\&PartialD;x}_j}+S_T\end{array}\right.$

Wherein, ${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})-\frac{{2\mathrm{\&mu;\&delta;}}_{\mathrm{ij}}{\&PartialD;u}_i}{{3\&PartialD;x}_j};$

I, j=1,2,3; u _iand u _jfor fluid velocity, x _ifor flow field variable, t is the time, and ρ is fluid density, and P is the pressure on fluid infinitesimal, τ _ijfor viscous stress, c _pfor specific heat capacity, T is temperature, the heat transfer coefficient that k is fluid, S _tfor default viscous dissipation item, μ is kinetic viscosity.

6. the method for numerical simulation of particle motion in alveolar contraction according to claim 5 and expansion process, it is characterized in that, described according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, while simulating described alveolar contraction and expansion, the step of IA flow field model comprises:

In the dynamic mesh model of described alveolar, simulate contraction and the expansion process of described alveolar according to following formula:

$\left\{\begin{array}{c}\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\\ \frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}\\ \frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l\end{array}\right.$

Wherein, $V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t},{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{{\mathrm{\&delta;V}}_l}{\mathrm{\&Delta;t}};$

ρ is fluid density, for fluid velocity, for the speed of motion grid node in described dynamic mesh model, Γ is coefficient of diffusion, for the border of control volume, for the area vector on border; V is the unit volume of motion grid, V ⁿfor the volume of motion grid described in the n moment, dV/dt controls the derivative of volume to the time, and n is the number of controlling volume face, S _φfor default broad sense source item, A _lbe the area vector of l face, δ V _lrepresent l face inswept volume in time Δ t, D is the diameter of alveolar place pipeline, and S is control area.

7. the method for numerical simulation of particle motion in alveolar contraction according to claim 5 and expansion process, is characterized in that, described motion process of monitoring described particle in described flow field model, and the step that obtains the exercise data of described particle comprises:

To motion control equation carry out the velocity of integration acquisition particle;

Wherein, the speed that u is air flow field, u _pfor the speed of particle, ρ _pfor the density of particle, the density that ρ is fluid, F _dfor acting on the power of dragging on particle, g _ifor gravitational vector, F _bthe Blang's power acting on particle;

$F_D=\frac{18\widetilde{\mathrm{\&mu;}}}{d_p^2{\mathrm{\&rho;}}_p{C}_C},C_C=1+\frac{2\mathrm{\&lambda;}}{d_p}[1.257+0.4\exp ({-1.1d}_p/2\mathrm{\&lambda;})],$ for the viscosity of fluid, C _cfor Cunningham slip correction coefficient, the mean free path that λ is molecule, d _pfor the diameter of particle;

ζ _ibe a Gaussian random vector with zero mean and unit variance, S ₀for the amplitude of white-noise process, the kinetic viscosity that v is air, k _b=1.38 × 10 ^-16for Boltzmann constant, the temperature that T is air flow field;

To governing equation carry out integration, obtain the position of particle in described flow field model described in each moment;

Wherein, u _pfor described velocity.

8. a numerical simulation system for particle motion in alveolar contraction and expansion process, is characterized in that, comprising:

Set up module, for set up the geometric model of alveolar according to the geometric parameter of default alveolar, the geometric model of described alveolar is carried out to grid division, generate the dynamic mesh model of described alveolar;

Acquisition module, for obtaining particle data, and corresponding respiration parameter and the fluid parameter of described alveolar;

Analog module, for according to the dynamic mesh model of described alveolar and geometric parameter, respiration parameter, fluid parameter and particle data, simulates that described alveolar shrinks and IA flow field model when expansion;

Monitoring modular, for monitor the motion process of described particle at described flow field model, obtains the exercise data of described particle.

Alveolar according to claim 8 shrink and expansion process in the numerical simulation system that moves of particle, it is characterized in that, described analog module also for:

In the dynamic mesh model of described alveolar, simulate described IA flow field according to following formula:

$\left\{\begin{array}{c}\frac{{\&PartialD;u}_i}{{\&PartialD;x}_i}=0\\ \frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{{\&PartialD;x}_j}=-\frac{\&PartialD;P}{{\&PartialD;x}_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{{\&PartialD;x}_j}\\ \frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left({\mathrm{\&rho;u}}_{j}T\right)}{{\&PartialD;x}_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P{\&PartialD;x}_j}\right)}{{\&PartialD;x}_j}+S_T\end{array}\right.$

Wherein, ${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})-\frac{{2\mathrm{\&mu;\&delta;}}_{\mathrm{ij}}{\&PartialD;u}_i}{{3\&PartialD;x}_j};$

I, j=1,2,3; u _iand u _jfor fluid velocity, x _ifor flow field variable, t is the time, and ρ is fluid density, and P is the pressure on fluid infinitesimal, τ _ijfor viscous stress, c _pfor specific heat capacity, T is temperature, the heat transfer coefficient that k is fluid, S _tfor default viscous dissipation item, μ is kinetic viscosity;

And/or for

In the dynamic mesh model of described alveolar, simulate contraction and the expansion process of described alveolar according to following formula:

$\left\{\begin{array}{c}\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\\ \frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}\\ \frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l\end{array}\right.$

Wherein, $V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t},{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{{\mathrm{\&delta;V}}_l}{\mathrm{\&Delta;t}};$

ρ is fluid density, for fluid velocity, for the speed of motion grid node in described dynamic mesh model, Γ is coefficient of diffusion, for the border of control volume, for the area vector on border; V is the unit volume of motion grid, V ⁿfor the volume of motion grid described in the n moment, dV/dt controls the derivative of volume to the time, and n is the number of controlling volume face, S _φfor default broad sense source item, A _lbe the area vector of l face, δ V _lrepresent l face inswept volume in time Δ t, D is the diameter of alveolar place pipeline, and S is control area.

Alveolar according to claim 8 shrink and expansion process in the numerical simulation system that moves of particle, it is characterized in that, described monitoring modular also for:

To motion control equation carry out the velocity of integration acquisition particle;

Wherein, the speed that u is air flow field, u _pfor the speed of particle, ρ _pfor the density of particle, the density that ρ is fluid, F _dfor acting on the power of dragging on particle, g _ifor gravitational vector, F _bthe Blang's power acting on particle;

$F_D=\frac{18\widetilde{\mathrm{\&mu;}}}{d_p^2{\mathrm{\&rho;}}_p{C}_C},C_C=1+\frac{2\mathrm{\&lambda;}}{d_p}[1.257+0.4\exp ({-1.1d}_p/2\mathrm{\&lambda;})],$ for the viscosity of fluid, C _cfor Cunningham slip correction coefficient, the mean free path that λ is molecule, the diameter that dp is particle;

ζ _ibe a Gaussian random vector with zero mean and unit variance, S ₀for the amplitude of white-noise process, the kinetic viscosity that v is air, k _b=1.38 × 10 ^-16for Boltzmann constant, the temperature that T is air flow field;

To governing equation carry out integration, obtain the position of particle in described flow field model described in each moment;

Wherein, u _pfor described velocity.