CN104090996A - Method for simulating air flow field in pulmonary alveolus - Google Patents

Abstract

A method for simulating an air flow field in a pulmonary alveolus comprises the following steps: mesh modeling is performed on the pulmonary alveolus according to geometric data of the pulmonary alveolus, a dynamic mesh model of the pulmonary alveolus is built, and initial air flow field pressure and mass flux of the dynamic mesh model of the pulmonary alveolus are set; a momentum model of the air flow field is solved according to the set initialized air flow field pressure and mass flux , and a first speed of the air flow field is determined according to the momentum model; a modified model is obtained according to a continuity equation of the air flow field, and a modified pressure value, a modified speed value and a modified mass flux value are acquired through verification; whether a convergence condition is met is determined according to the preset convergence condition, if yes, an air flow field simulation result is acquired, otherwise, the initialized air flow field pressure and mass flux are replaced with the modified pressure value and the modified mass flux value, and the modified pressure value, the modified speed value and the modified mass flux value are re-calculated until the convergence condition is met. With adoption of the scheme, the air flow field simulation accuracy is increased.

Description

Air flow field analogy method in alveolar

Technical field

The present invention relates to alveolar analog detection technical field, particularly relate to air flow field analogy method in a kind of alveolar.

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 gas exchange area 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 with suck Particulate Pollution and have direct or indirect contact.Therefore; the kinetic characteristic of research pellet in respiratory tract, for helping pathogenesis and the gasoloid treatment of understanding inhalable particles to have very important meaning, for the health of to protect mankind; protection of the environment, improves the quality of living and has positive effect.

In conventional art, adopt the mode of numerical simulation first to simulate IA air flow field, then obtain the kinetic characteristic of particle based on numerical simulation.Be the tracheal bronchus model of simplifying by foundation at present, by computational fluid dynamics, Air Flow wherein simulated, simulation value and actual conditions that this analog form obtains differ larger.

Summary of the invention

Based on this, being necessary, for the low problem of simulation accuracy rate, provides air flow field analogy method in a kind of alveolar.

Air flow field analogy method in a kind of alveolar, comprising:

According to the geometric data of alveolar, alveolar is carried out to mesh modeling, set up alveolar dynamic mesh model, the initialized air flow field pressure of alveolar dynamic mesh model and mass flux are set;

Solve the momentum model of air flow field according to described initialized air flow field pressure and mass flux, and determine the First Speed of air flow field according to described momentum model;

Obtain correction model according to the continuity equation of air flow field, described First Speed, initialized air flow field pressure and mass flux substitution correction model are revised, obtain pressure modified value, speed modified value and mass flux modified value;

According to the default condition of convergence, judge whether described pressure modified value, speed modified value and mass flux modified value meet the condition of convergence, if, solve the momentum model of air flow field according to pressure modified value and mass flux modified value, obtain the second speed of air flow field, pressure modified value, second speed and mass flux modified value are made as to the result of air flow field simulation; If not, described pressure modified value, mass flux modified value are replaced to initialized air flow field pressure and mass flux, recalculate pressure modified value, speed modified value and mass flux modified value, until meet the condition of convergence.

Air flow field analogy method in above-mentioned alveolar, by according to actual alveolar geometric data, alveolar being carried out to mesh modeling, has improved the similarity of model and human body alveolar greatly.Solve the First Speed of the momentum model acquisition air flow field of air flow field by initialized air flow field pressure and mass flux, then revise according to correction model, obtain pressure modified value, speed modified value and mass flux modified value, to meet continuity equation.Judge whether convergence by the default condition of convergence, thereby can obtain the air flow field speed of simulation.Improve the accuracy and the accuracy of obtaining air flow field speed of air flow field simulation, and then improved and obtain the accuracy of alveolar Kinematic Locus, thereby can prepare the good gasoloid of the property of medicine according to movement locus.

Brief description of the drawings

Fig. 1 is the schematic flow sheet of air flow field analogy method embodiment in alveolar of the present invention;

Fig. 2 is the front schematic view of alveolar dynamic mesh model in the embodiment of the present invention;

Fig. 3 is the view of alveolar dynamic mesh model in the embodiment of the present invention;

Fig. 4 is air flow rate schematic diagram when unstable state in respiratory in the embodiment of the present invention;

Fig. 5 is the two-dimensional geometry model of alveolar sac in the embodiment of the present invention;

Fig. 6 is steady-state gas flow contoured velocity YZ cross sectional representation in the embodiment of the present invention;

Fig. 7 is steady-state gas flow contoured velocity XZ schematic cross-section in the embodiment of the present invention;

Fig. 8 is steady-state gas flow contoured velocity XY schematic cross-section in the embodiment of the present invention;

Fig. 9 is steady-state gas flow contoured velocity XZ-XY schematic cross-section in the embodiment of the present invention;

Figure 10 is steady-state gas flow contoured velocity X-Y-Z tri-schematic cross-sections in the embodiment of the present invention;

Figure 11 is the motion pattern of alveolar inside on first order tracheae in the static situation of 1/10T moment alveolar wall in the embodiment of the present invention;

Figure 12 is the motion pattern of alveolar inside on first order tracheae in 1/10T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention;

Figure 13 is the motion pattern of alveolar inside on first order tracheae in the static situation of 1/5T moment alveolar wall in the embodiment of the present invention;

Figure 14 is the motion pattern of alveolar inside on first order tracheae in 1/5T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention;

Figure 15 is the motion pattern of alveolar inside on first order tracheae in the static situation of 3/10T moment alveolar wall in the embodiment of the present invention;

Figure 16 is the motion pattern of alveolar inside on first order tracheae in 3/10T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention;

Figure 17 is the motion pattern of alveolar inside on first order tracheae in the static situation of 2/5T moment alveolar wall in the embodiment of the present invention;

Figure 18 is the motion pattern of alveolar inside on first order tracheae in 2/5T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention;

Figure 19 be in the embodiment of the present invention in the air-breathing peak value moment first order motion pattern of alveolar inside;

Figure 20 be in the embodiment of the present invention on air-breathing peak value moment third level the motion pattern of alveolar inside;

Figure 21 be in the embodiment of the present invention on air-breathing peak value moment level V the motion pattern of alveolar inside;

Figure 22 be in the embodiment of the present invention on the 7th grade of air-breathing peak value moment the motion pattern of alveolar inside;

Figure 23 be in the embodiment of the present invention in the peak value moment first order of exhaling the motion pattern of alveolar inside;

Figure 24 is the motion pattern in the peak value moment of exhaling alveolar inside on the third level in the embodiment of the present invention;

Figure 25 be in the embodiment of the present invention on the peak value moment level V of exhaling the motion pattern of alveolar inside;

Figure 26 be in the embodiment of the present invention on the 7th grade of peak value moment of exhaling the motion pattern of alveolar inside;

Figure 27 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/10T moment alveolar inside;

Figure 28 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/10T moment alveolar inside;

Figure 29 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/5T moment alveolar inside;

Figure 30 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/5T moment alveolar inside;

Figure 31 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 3/10T moment alveolar inside;

Figure 32 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 3/10T moment alveolar inside;

Figure 33 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 2/5T moment alveolar inside;

Figure 34 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 2/5T moment alveolar inside.

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 the schematic flow sheet of air flow field analogy method embodiment in alveolar of the present invention, comprise step:

Step S101: according to the geometric data of alveolar, alveolar is carried out to mesh modeling, set up alveolar dynamic mesh model, the initialized air flow field pressure of alveolar dynamic mesh model and mass flux are set;

The geometric data of alveolar can refer to the geometric data of a certain joint alveolar.Described geometric data can comprise: pipe range parameter, intracavity diameter, alveolar diameter, the centre distance between two alveolars, alveolar angular aperture, respiration parameter.Such as, can be Section 18 alveolar.Geometric data can comprise: pipe range is 740 μ m, and intracavity diameter is 320 μ m, and alveolar diameter is 320 μ m, and the centre distance between two alveolars is 340 μ m, and alveolar angular aperture is 120 °, and respiration parameter is 0.283.

The present embodiment also adopts dynamic mesh model simulation alveolar, improves authenticity and the correlativity of simulation.

Step S102: solve the momentum model of air flow field according to described initialized air flow field pressure and mass flux, and determine the First Speed of air flow field according to described momentum model;

Momentum model can obtain according to the law of conservation of momentum.Can obtain air flow field speed according to momentum model.Such as momentum model can be:

$\frac{\∂\left({\mathrm{\ρu}}_i\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρu}}_i{u}_j\right)}{{\∂x}_j}=-\frac{\∂P}{{\∂x}_i}+\frac{\∂\left({\mathrm{\τ}}_{\mathrm{ij}}\right)}{{\∂x}_j}$

ρ represents fluid density, u _irepresent fluid i direction speed, u _jrepresent fluid j direction speed, x _irepresent fluid i direction position, x _jrepresent fluid j direction position, wherein, i is 1 or 2 or 3, represents respectively disalignment, and j is 1 or 2 or 3, represents respectively disalignment, and P represents the pressure on fluid infinitesimal, τ _ijrepresent viscous stress,

${\mathrm{\τ}}_{\mathrm{ij}}=\mathrm{\μ}(\frac{\∂u_i}{\∂x_j}+\frac{\∂u_j}{\∂x_i})-\frac{2\mathrm{\μ}{\mathrm{\δ}}_{\mathrm{ij}}\∂u_i}{3\∂x_j}$

μ represents the viscosity (claiming again the kinematic viscosity of fluid) of fluid, δ _ijrepresent the kinetic viscosity of fluid.

Step S103: obtain correction model according to the continuity equation of air flow field, described First Speed, initialized air flow field pressure and mass flux substitution correction model are revised, obtain pressure modified value, speed modified value and mass flux modified value;

Correction model can be the Poisson type equation of being derived by continuity equation, is used for the parameters such as pressure to revise, to make revised data meet continuity equation.Wherein, continuity equation can be

Step S104: according to the default condition of convergence, judge whether described pressure modified value, speed modified value and mass flux modified value meet the condition of convergence, if so, enter step S105, if not, enter step S106;

The condition of convergence can be respectively pressure, First Speed, mass flux when fore pressure modified value, speed modified value and mass flux modified value and before revising to be carried out to error ratio, when error is less than corresponding setting value simultaneously, represent to meet the condition of convergence, otherwise, foot with thumb down.

Step S105: solve the momentum model of air flow field according to pressure modified value and mass flux modified value, obtain the second speed of air flow field, pressure modified value, second speed and mass flux modified value are made as to the result of air flow field simulation;

In the time that the condition of convergence meets, can obtain the speed of air flow process.

Step S106: described pressure modified value, mass flux modified value are replaced to initialized air flow field pressure and mass flux, return to step S102.

In the time that the condition of convergence does not meet, the pressure modified value of acquisition, mass flux modified value are replaced to initialized air flow field pressure and mass flux as current data, return to step S102 and calculate, until the condition of convergence meets.

The present embodiment, by according to actual alveolar geometric data, alveolar being carried out to mesh modeling, has improved the similarity of model and human body alveolar greatly.Solve the First Speed of the momentum model acquisition air flow field of air flow field by initialized air flow field pressure and mass flux, then revise according to correction model, obtain pressure modified value, speed modified value and mass flux modified value, to meet continuity equation.Judge whether convergence by the default condition of convergence, thereby can obtain the air flow field speed of simulation.Improve the accuracy of obtaining air flow field speed, and then improved and obtain the accuracy of alveolar Kinematic Locus, thereby can prepare the good gasoloid of the property of medicine according to movement locus.

Set up alveolar dynamic mesh model, can also upgrade dynamic mesh model.The method of upgrading has a variety of, such as elasticity theory of adjustment (Spring-based Smoothing) and Local grid are heavily drawn method (Local Remeshing Method).In elasticity theory of adjustment, can adopt again seamed edge spring.

In seamed edge spring, 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{\→}{F}}_i=\underset{j}{\overset{n_i}{\mathrm{\Σ}}}{k}_{\mathrm{ij}}(\mathrm{\Δ}{\stackrel{\→}{x}}_j-\mathrm{\Δ}{\stackrel{\→}{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 iterative by the balance equation group to system node:

$\mathrm{\Δ}{\stackrel{\→}{x}}_i^{m+1}=\frac{{\mathrm{\Σ}}_i^{n_i}{k}_{\mathrm{ij}}\mathrm{\Δ}{\stackrel{\→}{x}}_j^m}{{\mathrm{\Σ}}_j^{n_i}{k}_{\mathrm{ij}}}$

Wherein on border, the displacement of node is 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:

${\stackrel{\→}{x}}_i^{n+1}={\stackrel{\→}{x}}_i^n+\mathrm{\Δ}{\stackrel{\→}{x}}_i^{\mathrm{converged}}$

The balance length of describing medi-spring at summit spring is zero.In order to make grid still can keep rational density to distribute after distortion, the coefficient of stiffiness of spring can have multiple choices, wherein can suppose that certain power of the length of side between the coefficient of stiffiness of spring and two nodes is inversely proportional to.

In the present embodiment, the elasticity coefficient value of seamed edge is:

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

Two dimension triangular mesh and three-dimensional tetrahedral grid can use elasticity smoothing method.

Local grid is heavily drawn 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 cannot restrain.Disposal route is that grid excessive aberration rate or that change in size is too violent is 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.Local grid heavily the method for drawing can only be applied to two-dimentional triangle and three-dimensional tetrahedral grid unit.Above-mentioned two kinds of methods can realize in Fluent.

The present invention program also provides a kind of preferred update method, and the method is applicable to hexahedral mesh, wedge shape grid etc. can be at the grid system of border higher slice.Specifically comprise:

According to the geometric data of alveolar, alveolar is carried out after mesh modeling, detect the height and the default difference of optimizing height value of adjacent boundary grid, in the time that difference is greater than the first setting value, between this adjacent boundary grid and boundary layer, increase one deck grid, in the time that difference is less than the second setting value, delete this adjacent boundary grid, wherein, described adjacent boundary grid is one deck grid adjacent with boundary layer.

This method is on border, to suppose the clathrum height value of an optimization, 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.

In an embodiment, the present invention program can also adopt the automatic renewal of the grand control dynamic mesh of DEFINE_GRID_MOTION therein.Define the enlargement and contraction motion in time of alveolar wall in respiratory by DEFINE_GRID_MOTION is grand, its basic format is

DEFINE_GRID_MOTION(name，d，dt，time，dtime)，

Name is the title that user specifies, and d and dt are two pointers that point to respectively zoning and dynamic mesh region, and time and dtime are respectively current time and the time step in flow field.

Described alveolar dynamic mesh model is 3D model.In an embodiment, in described alveolar dynamic mesh model, grid can be one or more in tetrahedron, hexahedron, sphenoid etc. therein.Wherein, the trellis-type of alveolar dynamic mesh model can be hybrid grid, and grid number can be 310,000 left and right.As shown in Figures 2 and 3, be front schematic view and the view of alveolar dynamic mesh model in the embodiment of the present invention.

Numerical simulation adopts three-dimensional steady state and unstable state to calculate, and adopts SIMPLE method to solve N-S equation to the calculating of air-flow.The air-flow calculating adopts laminar flow (laminar) computation model, can utilize the DPM model in Fluent to follow the trail of movement of particles track simultaneously.When the inlet velocity of air-flow is pressed adult's eupnea, such as, the flow of Section 18 alveolar is given.The respiration parameter of the Section 18 under adult normal sitting state: Re=0.283, the parameter of air is: ρ=1.225kg/m ³, μ=1.7894e-05kg/m.s, ν=μ/ρ=1.46e-05m ²/ s, Re=UD/ ν, U=0.0129m/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, $Q_{\max }=\stackrel{\‾}{u}\×A\×\mathrm{\π}/2,\stackrel{\‾}{u}=0.129m/s,$ Try to achieve $u_{\max }=\stackrel{\‾}{u}\×\mathrm{\π}/2.$ As shown in Figure 4, be air flow rate schematic diagram when unstable state in respiratory.

In calculating, calculated the air-flow of a complete cycle in stable state and the transient of IA mobile 3D model, the calculating of particle is 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/Nt100%.Nd represents the precipitation number of particle at wall, and Nt represents total particle number.

Therein in an embodiment, setting up in alveolar dynamic mesh model process, also comprise: the two-dimensional geometry model of setting up alveolar sac, be used for the enlargement and contraction of the alveolar of simulated respiration process, wherein, by the tracheae in the cylinder tube shape simulated respiration road of front and back opening, by the alveolar in the round simulation respiratory tract being connected on cylinder wall, alveolar is positioned at the center position of this section of tracheae.Correlation parameter can also be set, and such as the radius of tracheae is 250um, the radius of alveolar is 200um, and the angular aperture of alveolar on main tracheae is made as 60 degree, and the length of tracheae is 1000um.

As shown in Figure 5, be the two-dimensional geometry model of alveolar sac.The present embodiment can utilize parameter model to carry out the enlargement and contraction of alveolar in simulated respiration process.Tracheae before and after using 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.Data in incorporating parametric model, establish the radius R of tracheae in the present embodiment model _dfor 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 computation model herein, alveolar is positioned at the center position of this section of tracheae.

In an embodiment, 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 therein.In the two-dimensional geometry model process of setting up alveolar sac, also comprise that each characteristic dimension that alveolar and tracheae be set is time dependent sine function:

$L\left(t\right)=L_0[1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2}\left)\right]$

Wherein, L (t) represents the value of the characteristic dimension of alveolar or tracheae, L ₀the value of each characteristic dimension in expression geometric model in the time of air-breathing beginning or while exhaling end, f represents the frequency of breathing, f=2 π/T, T represents the cycle length of breathing, β represents the expansion amplitude of each characteristic dimension, β=(C+1) ^1/3-1, C=(V _max-V _min)/V _min, V _maxand V _minrepresent maximum volume and the minimum volume of alveolar dynamic mesh model, described maximum volume and minimum volume respectively corresponding expiration in actual respiratory end and air-breathing volume at the end.

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:

Therein in an embodiment, in the two-dimensional geometry model process of setting up alveolar sac, also comprise: determine alveolar present position on respiratory tract according to the gas flow that passes in and out alveolar in respiratory with the ratio of the gas flow that flows through alveolar tracheae of living in, set up the two-dimensional geometry model of alveolar sac according to described position, wherein:

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\π\βf}{R}_a^2{\mathrm{\λ}}^2\cos (\mathrm{ft}-\frac{\mathrm{\π}}{2}),\mathrm{\λ}=1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2})$

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

Q _a(t) represent to pass in and out in respiratory the gas flow of alveolar, R _arepresent alveolar radius, Q _d(t) gas flow of alveolar tracheae of living in, V are flow through in expression _trepresent the value through the alveolate volume summation acquisition of institute on this section of retrotracheal all tracheae and tracheae.

In the present embodiment 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 present embodiment model, 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.

In the model of the present embodiment, the deflection on the border of alveolar is greater than the height that is positioned at grid on border, can control by elasticity smoothing method and the Local Gravity And method of drawing the amoeboid movement of grid in model, therefore divide and can use triangular mesh for the grid of the computation model of the present embodiment.

This programme is based on above-mentioned analogy method, and the 3D numerical simulation result that carries out unstable state airflow characteristic in the 3D numerical simulation result of alveolar homeostasis airflow characteristic and analysis, alveolar is with analysis, periodically IA flow field analysis is descended in the motion of wall shrinkage expansion.

First, carry out 3D numerical simulation result and the analysis of alveolar homeostasis airflow characteristic.

In Fig. 6 to Figure 10, provide respectively vertical X, Y, the contoured velocity on Z axis interface.Fig. 6 is steady-state gas flow contoured velocity YZ cross sectional representation in the embodiment of the present invention, and Fig. 7 is steady-state gas flow contoured velocity XZ schematic cross-section in the embodiment of the present invention, and Fig. 8 is steady-state gas flow contoured velocity XY schematic cross-section in the embodiment of the present invention; Fig. 9 is steady-state gas flow contoured velocity XZ-XY schematic cross-section in the embodiment of the present invention; Figure 10 is steady-state gas flow contoured velocity X-Y-Z tri-schematic cross-sections in the embodiment of the present invention.

As can be seen from the figure, it is negligible that the intra-acinous speed of lung is compared with the speed in main airway, and in the main airway after lung acinus, maximal value appears in speed, in cylinder main airway, speed is extended with increase tendency along alveolar, and velocity amplitude is along with radius reduces gradually; Meanwhile, XZ cross section is the same with the velocity distribution on XY interface, and visible air current flow is not subject to the impact of gravity.

The IA streamline of 3D is the intra-acinous streamline complex structure of lung particularly.In alveolar, have equally the existence of recirculating zone, but not only there is nonplanar structure recirculating zone, and has the backflow of spatial structure, in alveolar, from outside to inside, recirculating zone is changed spatial structure into from planar structure gradually, refluxes and also becomes and become increasingly complex.In cylinder main airway, it is simple many that air-flow streamline just becomes, and air-flow is substantially parallel with X-axis, in acinus porch, produces gradually to streaming in bubble.Meanwhile, between primary air and backflow, there is not convective exchange, finally still lay out from opening part around the air-flow entering, illustrate between recirculating zone and primary air and have a division surface, enter backflow around entering air-flow.

Then, carry out 3D numerical simulation result and the analysis of unstable state airflow characteristic in alveolar.

For unstable state airflow characteristic is analyzed, can choose 1/8T, 1/4T, 3/8T, 1/2T, 5/8T, 3/4T and the airflow characteristic in T moment compare analysis.Can find out from analyzing, except breathing conversion moment (1/2T and T moment), in other moment, gas velocity is similar in IA distribution, and just the size of speed is not identical.It needs to be noted, the velocity distribution in lateral cross section is the same, on longitudinal section, no matter is expiration or air-breathing, and speed maximal value is all at pressure export place, and the speed of speed porch speed in the time of incoming call and exhalation is all uniform.Meanwhile, the air-flow in main airway is compared with intra-acinous air-flow, and intra-acinous air-flow is negligible.Speed streamline on cross section can be found out, has recirculating zone in acinus, and the position of recirculating zone converts, and has equally separator bar (face) between recirculating zone and primary air simultaneously.In the moment (1/2T and T moment) of expiration and air-breathing conversion, there is obvious change in air-flow, original airflow field is carved at this moment essential variation has been occurred, between this moment recirculate mixing primary air, obvious material convective exchange has occurred, disturbance that now can air-flow is also maximum.

In analyzing according to the 1/4T moment, can find out that each intra-acinous streamline is identical, the many of complexity are wanted in intra-acinous backflow, between acinus backflow simultaneously and primary air, have separation interface.The intra-acinous backflow center in each moment is that variation is slightly occurring, this shows the flow perturbation that unstable state air-flow can generating portion, make, between acinus backflow and primary air, slight convective exchange occurs, thereby produced the variation at the backflow center in expiration and breathing process.But this disturbance effect will be much smaller more than breathing the flow perturbation occurring in transfer process.Meanwhile, due to the complexity of recirculating zone, therefore can not significantly find out whether primary air, around the radian that enters acinus, obvious change has occurred along with respiratory.

Finally, carry out the lower IA flow field analysis of periodicity wall shrinkage expansion motion.

First analyze in flow process the situation that in model, tracheae and respiratory tract are not expanded, calculated the unstable state respiratory without wall motion.In this course, the inlet velocity of gas is still time dependent sinusoidal wave distribution.

Figure 11 is the motion pattern of alveolar inside on first order tracheae in the static situation of 1/10T moment alveolar wall in the embodiment of the present invention, Figure 12 is the motion pattern of alveolar inside on first order tracheae in 1/10T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention, wherein, alveolar ectasia coefficient is 0.05.

Figure 13 is the motion pattern of alveolar inside on first order tracheae in the static situation of 1/5T moment alveolar wall in the embodiment of the present invention, Figure 14 is the motion pattern of alveolar inside on first order tracheae in 1/5T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention, wherein, alveolar ectasia coefficient is 0.05.

Figure 15 is the motion pattern of alveolar inside on first order tracheae in the static situation of 3/10T moment alveolar wall in the embodiment of the present invention, Figure 16 is the motion pattern of alveolar inside on first order tracheae in 3/10T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention, wherein, alveolar ectasia coefficient is 0.05.

Figure 17 is the motion pattern of alveolar inside on first order tracheae in the static situation of 2/5T moment alveolar wall in the embodiment of the present invention, Figure 18 is the motion pattern of alveolar inside on first order tracheae in 2/5T moment alveolar wall expansion contractile motion situation in the embodiment of the present invention, wherein, alveolar ectasia coefficient is 0.05.

Can find, in the situation that not having alveolate wall, the endotracheal alveolar opening part that is flowing in has formed a shear flow, has formed the flow region of a convolution in alveolar region.Whole flow field can be divided into two regions, is respectively endotracheal and flows and IA flowing, and can use approx the line of two jiaos of connection alveolar opening part front and back separately.For the three-dimensional alveolar structure of reality, this separator bar can be regarded a curved surface as.The effect that endotracheal fluid and IA fluid almost do not have convection current to mix.Endotracheal air-flow is made convective motion towards Way out, and the fluid of alveolar inside is done circumnutation to be parallel to alveolar wall, and the region of circumnutation occupies whole alveolar cavity, and this is the outer slightly principal character of cavity shear flow.

The mobile speed of IA convolution is very little, is far smaller than endotracheal flowing velocity, only has the centesimal magnitude of mainstream speed.For axisymmetric cavity, when the speed of shear flow is very little, convolution region should be also axisymmetric.Above the first order because the gas velocity in main flow tracheae is also larger, the partially upstream side of the center in the region of therefore circling round.

The impact of moving on its interior flow field in order to analyze alveolar residing position and alveolar wall on respiratory tract.In the first order, although the gas flow of alveolar is with the gas flow (Q of tracheae _a/ Q _d=0.004%) comparing, is almost a negligible value.But with not having the flow field structure of alveolar inside under alveolar wall motion conditions to compare, can see that the Field Characteristics of alveolar inside has had obvious difference.Shear flow in tracheae has still formed a very large recirculation zone in alveolar inside.But the motion stream field of wall has produced effect: the fluid of endotracheal fraction, with IA fluid, convection current mixing has occurred, this part fluid comes from the position of tracheae near wall.Tangential speed of fluid of doing circumnutation in alveolar has been given in the expansion of alveolar, give radial velocity of endotracheal gas, tracheae is come near fluid and the disengaging of endotracheal main flow of wall, this segment fluid flow enters after alveolar with convective motion, transfer near the locality of outlet at alveolar opening, and move along the wall of alveolar.The form of the convolution flow region of alveolar inside also has a very large change.Do not have in the flow field of alveolar wall motion, can see that it is axisymmetric that convolution is flowed.In alveolar, the mobile center of convolution is also constantly offset, and the offset direction in direction and the breathing process of this skew is contrary, identical with endotracheal gas flow direction, and the area in convolution region is along with skew reduces.

Figure 19 be in the embodiment of the present invention in the air-breathing peak value moment first order motion pattern of alveolar inside, wherein respiratory intensity 10L/min.Figure 20 be in the embodiment of the present invention on air-breathing peak value moment third level the motion pattern of alveolar inside, wherein respiratory intensity 10L/min.Figure 21 be in the embodiment of the present invention on air-breathing peak value moment level V the motion pattern of alveolar inside, wherein respiratory intensity 10L/min.Figure 22 be in the embodiment of the present invention on the 7th grade of air-breathing peak value moment the motion pattern of alveolar inside, wherein respiratory intensity 10L/min.

As can be seen from the figure, the wall of alveolar is in extensional motion, and in alveolar, the mobile off-centring of convolution has arrived the position of close alveolar opening near tracheae upstream.In breathing process, move towards the axis of alveolar at the mobile center of circling round, and mobile and main endotracheal gas flow direction is identical, the area change in convolution region, and in the time of air-breathing end, convolution has relatively approached the axis of alveolar.

Figure 23 be in the embodiment of the present invention in the peak value moment first order of exhaling the motion pattern of alveolar inside, wherein respiratory intensity 10L/min.Figure 24 is the motion pattern in the peak value moment of exhaling alveolar inside on the third level in the embodiment of the present invention, wherein respiratory intensity 10L/min.Figure 25 be in the embodiment of the present invention on the peak value moment level V of exhaling the motion pattern of alveolar inside, wherein respiratory intensity 10L/min.Figure 26 be in the embodiment of the present invention on the 7th grade of peak value moment of exhaling the motion pattern of alveolar inside, wherein respiratory intensity 10L/min.

In the process of exhaling, alveolar is in contractile motion, and at this moment the direction of IA streamline and breathing process synchronization is just in time contrary.Although there is convection current mixing in the fluid of tracheae inside and IA fluid, but only occur on the fraction fluid that approaches tube wall in tracheae, most of fluid in main flow, with IA fluid, convection current does not occur to be mixed, this part fluid, towards tracheae Way out convective motion, enters the next stage of tracheae.Along with the position of alveolar on respiratory tract approaches terminal, in alveolar and tracheae, the convection action of air-flow is just more obvious.

Work as Q _a/ Q _dwhile increasing to a certain value, in whole respiratory, all there will not be convolution region in alveolar, can observe, the 7th grade of alveolar has still occurred therefore in the most alveolars on human respiratory tract, all having circumnutation region in convolution region above.

Under different alveolar ectasia coefficients, the interior slightly difference that flows also of same alveolar.Flare factor is increased to 0.1 by 0.05, the expansion amplitude of alveolar doubles, now endotracheal flow has just in time also increased by one times, it is more obvious that in tracheae, mainstream gas speed increases the effect producing, from the contrast of Figure 27 to Figure 34, can find out, tracheae mainstream speed generation effect increase is more obvious.Figure 27 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/10T moment alveolar inside; Figure 28 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/10T moment alveolar inside.

Figure 29 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/5T moment alveolar inside; Figure 30 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 1/5T moment alveolar inside.

Figure 31 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 3/10T moment alveolar inside; Figure 32 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 3/10T moment alveolar inside.

Figure 33 be when in the embodiment of the present invention, respiratory intensity is 10L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 2/5T moment alveolar inside; Figure 34 be when in the embodiment of the present invention, respiratory intensity is 20L/min the alveolar of the 7th grade at the motion pattern of expiratory phase 2/5T moment alveolar inside.

Various technical characterictics in above embodiment can combine arbitrarily, as long as there is not conflict or contradiction in the combination between feature, but as space is limited, describe one by one, therefore combining arbitrarily of the various technical characterictics in above-mentioned embodiment also belongs to this instructions scope of disclosure.

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 (8)

1. an air flow field analogy method in alveolar, is characterized in that, comprising:

According to the geometric data of alveolar, alveolar is carried out to mesh modeling, set up alveolar dynamic mesh model, the initialized air flow field pressure of alveolar dynamic mesh model and mass flux are set;

Solve the momentum model of air flow field according to described initialized air flow field pressure and mass flux, and determine the First Speed of air flow field according to described momentum model;

Obtain correction model according to the continuity equation of air flow field, described First Speed, initialized air flow field pressure and mass flux substitution correction model are revised, obtain pressure modified value, speed modified value and mass flux modified value;

According to the default condition of convergence, judge whether described pressure modified value, speed modified value and mass flux modified value meet the condition of convergence, if, solve the momentum model of air flow field according to pressure modified value and mass flux modified value, obtain the second speed of air flow field, pressure modified value, second speed and mass flux modified value are made as to the result of air flow field simulation; If not, described pressure modified value, mass flux modified value are replaced to initialized air flow field pressure and mass flux, recalculate pressure modified value, speed modified value and mass flux modified value, until meet the condition of convergence.

2. air flow field analogy method in alveolar according to claim 1, is characterized in that, described geometric data comprises: pipe range parameter, intracavity diameter, alveolar diameter, the centre distance between two alveolars, alveolar angular aperture, respiration parameter.

3. air flow field analogy method in alveolar according to claim 1, is characterized in that, while setting up alveolar dynamic mesh model, also comprises dynamic mesh model is upgraded, and comprising:

According to the geometric data of alveolar, alveolar is carried out after mesh modeling, detect the height and the default difference of optimizing height value of adjacent boundary grid, in the time that difference is greater than the first setting value, between this adjacent boundary grid and boundary layer, increase one deck grid, in the time that difference is less than the second setting value, delete this adjacent boundary grid, wherein, described adjacent boundary grid is one deck grid adjacent with boundary layer.

4. air flow field analogy method in alveolar according to claim 1, is characterized in that, in described alveolar dynamic mesh model, grid comprises one or more in tetrahedron, hexahedron, sphenoid.

5. air flow field analogy method in alveolar according to claim 1, is characterized in that, the trellis-type of described alveolar dynamic mesh model is hybrid grid, and grid number is 310,000.

6. air flow field analogy method in alveolar according to claim 1, is characterized in that, setting up in alveolar dynamic mesh model process, also comprises:

Set up the two-dimensional geometry model of alveolar sac, be used for the enlargement and contraction of the alveolar of simulated respiration process, wherein, by the tracheae in the cylinder tube shape simulated respiration road of front and back opening, by the alveolar in the round simulation respiratory tract being connected on cylinder wall, alveolar is positioned at the center position of this section of tracheae.

7. air flow field analogy method in alveolar according to claim 6, is characterized in that, in the two-dimensional geometry model process of setting up alveolar sac, also comprises that each characteristic dimension that alveolar and tracheae be set is time dependent sine function:

$L\left(t\right)=L_0[1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2}\left)\right]$

Wherein, L (t) represents the value of the characteristic dimension of alveolar or tracheae, L ₀the value of each characteristic dimension in expression geometric model in the time of air-breathing beginning or while exhaling end, f represents the frequency of breathing, f=2 π/T, T represents the cycle length of breathing, β represents the expansion amplitude of each characteristic dimension, β=(C+1) ^1/3-1, C=(V _max-V _min)/V _min, V _maxand V _minrepresent maximum volume and the minimum volume of alveolar dynamic mesh model, described maximum volume and minimum volume respectively corresponding expiration in actual respiratory end and air-breathing volume at the end.

8. air flow field analogy method in alveolar according to claim 7, it is characterized in that, in the two-dimensional geometry model process of setting up alveolar sac, also comprise: determine alveolar present position on respiratory tract according to the gas flow that passes in and out alveolar in respiratory with the ratio of the gas flow that flows through alveolar tracheae of living in, set up the two-dimensional geometry model of alveolar sac according to described position, wherein:

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\π\βf}{R}_a^2{\mathrm{\λ}}^2\cos (\mathrm{ft}-\frac{\mathrm{\π}}{2}),\mathrm{\λ}=1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2})$

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

Q _a(t) represent to pass in and out in respiratory the gas flow of alveolar, R _arepresent alveolar radius, Q _d(t) gas flow of alveolar tracheae of living in, V are flow through in expression _trepresent the value through the alveolate volume summation acquisition of institute on this section of retrotracheal all tracheae and tracheae.