A method for calculating respiratory airflow and particle transport based on a dynamic coupling full airway model

Abstract

The application provides a method for calculating respiratory airflow and particle transport based on a dynamic coupling full airway model. The method comprises: constructing a five-branch full airway three-dimensional model comprising an upper airway segment, a triple bifurcation unit (TBU) segment, and an alveolar segment; establishing a dynamic motion model of the glottis, trachea, and alveolar region to realize the geometric change of the airway structure in the respiratory cycle; adopting a segmented independent simulation and step-by-step boundary condition transmission strategy to complete the transient airflow simulation in each level of the airway; and introducing a particle tracking mechanism to realize the simulation of the transport, deposition, and exhalation process of inhaled particles from the oral cavity to the alveoli. The above method combines real anatomical structure, dynamic physiological motion, and efficient calculation strategy to improve the accuracy and calculation efficiency of airflow and particle transport simulation, and is suitable for scenarios such as inhaled drug delivery, air pollution exposure assessment, and respiratory system disease modeling.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of biomedical engineering and computational fluid dynamics (CFD), and in particular to a method for calculating respiratory airflow and particle transport based on a dynamically coupled full airway model. Background Technology

[0002] With the widespread application of computational fluid dynamics (CFD) in medical engineering, airway airflow simulation and particle transport simulation have become important tools in fields such as drug delivery, disease mechanism research, and environmental exposure assessment. Most existing CFD airway models are based on static, rigid wall assumptions, generally neglecting the dynamic motion characteristics of key airway structures during respiration. For example, in the glottis region, it is often simplified to a fixed cross-section channel, failing to reflect the periodic opening and closing changes caused by inspiratory expansion and expiratory contraction; the trachea is usually modeled as a smooth cylindrical structure, failing to consider the asymmetric geometry of the cartilaginous rings in real anatomy and its deformation response under intrathoracic pressure; and in deep lung regions (such as the terminal bronchi and alveoli), they are often omitted or replaced by symmetrical static branching structures, making it difficult to reflect the actual contraction and expansion processes of alveoli driven by respiration.

[0003] Furthermore, existing models often employ simplified boundary conditions and lack dynamic driving mechanisms coupled with human spontaneous respiratory rhythms (such as tidal volume, respiratory rate, and inspiratory-expiratory ratio). This significantly limits the ability to reproduce changes in airway structure over time during simulations and leads to significant deviations between simulation results such as airflow field distribution, wall shear stress, and particle deposition behavior and actual physiological states.

[0004] These limitations not only affect the accuracy of simulation results but also reduce the application value of CFD models in scenarios such as targeted delivery of inhaled drugs, prediction of pollutant exposure risks, and modeling of chronic respiratory diseases. Therefore, there is an urgent need for a CFD modeling method with higher physiological realism that can dynamically couple changes in airway structure with respiratory behavior to improve the reliability and applicability of simulations and provide a more meaningful simulation data foundation for subsequent clinical or engineering applications. Summary of the Invention

[0005] In view of this, the purpose of this invention is to provide a method for calculating respiratory airflow and particle transport based on a dynamically coupled full airway model, in order to solve the problems of static structure, lack of dynamic coupling, poor continuity of particle transport, and low simulation efficiency in existing airway CFD models, thereby achieving accurate simulation of airway geometric evolution, airflow behavior and particle transport and deposition processes under real physiological rhythm conditions.

[0006] To achieve the above objectives, the present invention provides the following technical solution: In one embodiment of the present invention, a method for calculating respiratory airflow and particle transport based on a dynamically coupled whole airway model is provided, the method comprising the following steps: Step 1: Construct a three-dimensional model of the five-branched airway, which includes the upper airway segment, the triple bifurcation unit (TBU) segment, and the acinar segment: The upper airway segment is reconstructed based on clinical CT scan data and includes a trachea with a cartilage ring structure. The cartilage ring is formed into annular grooves through sweeping excision and rounding treatment, and surface continuity is optimized using a geometric repair module. The TBU segment consists of three sets of triple bifurcation structures, corresponding to the G7–G9, G10–G12 and G13–G15 generation bronchi, respectively. Each TBU is modeled using the regular dichotomy method, and the geometric parameters are determined based on Weibel morphometry data. Each TBU segment is set with a random tilt angle in the range of 30° to 65°. The acinar segment is generated into a polyhedral alveolar structure using a three-dimensional Voronoi diagram. Each acinar unit contains approximately 500 alveoli, with an average equivalent alveolar diameter of 313 μm and an average path distance of 2.5 mm. The upper airway segment, TBU segment, and acinar segment are connected in a series and parallel manner to form five independent pathways. Each pathway extends from the G6 exit to the acinar area, corresponding to the five lobes of the lung.

[0007] Step 2: Establish a dynamic respiratory tract motion model to simulate the periodic deformation behavior of the glottis, trachea, and alveolar regions during inhalation and exhalation. The cross-sectional area change in the glottic region is simulated using a nodal displacement function based on Fourier series, and its displacement is expressed as: ; in: : Represents the horizontal coordinate of the k-th node on the m-th plane at time t; : This represents the initial position of the node; : A time-dependent Fourier series function that describes the trend of cross-sectional changes within the respiratory cycle; f(y) and F(z) represent the spatial distribution functions of the nodal displacements in the y and z directions, respectively. The ratio of the maximum opening to the static opening determines the maximum opening amplitude.

[0008] The tracheal segment was simulated using a dynamic mesh method to model the triaxial anisotropic deformation of the airway wall, with the nodal displacements being: ; in: : Represents the position coordinates (unit: mm) of the i-th node at time t along a certain direction (such as x, y, or z). : This is the reference position of the node in its static state; : The coordinates of the trachea's axis center in this direction, used to define the deformation center (unit: mm); : This represents the deformation rate of the node in the current direction (unit: dimensionless), defined as the ratio of the maximum deformation amplitude to the static distance; for example... = 1, = 0.375; : This is the current time step (in seconds), used for iterative control; : The total duration of a complete inhalation and exhalation cycle (unit: s), usually set to 6s, but can be adapted according to patient parameters; This expression can be set independently in the x, y, and z directions. The value simulates the triaxial anisotropic deformation behavior of the tracheal wall under the action of thoracic cavity pressure.

[0009] The expansion behavior of the acinar region under spontaneous breathing mode is modeled using a sine function:

[0010] in: : The path length at any given time; : Path length in resting state; : is the expansion coefficient, set to 0.053; ω is the angular frequency, and T is the respiratory cycle time; If it is a deep breathing mode, it is described by a combination of a sixth-order polynomial and a sine function; Step 3: Perform transient airflow calculations using a segmented independent simulation and step-by-step data transfer strategy: First, a transient CFD simulation was performed on the upper airway section to obtain the mass flow rate and velocity distribution at each outlet of G6. Set the TBU segment inlet at five selected representative G6 exits and load their speed data as G7 inlet conditions; The calculation results of each TBU segment are passed down from top to bottom as the entry boundary conditions of the next segment, all the way to the acinar region; During the exhalation phase, the inspiratory airflow field data is reversed according to the time node as the initial condition to simulate the reverse exhalation airflow.

[0011] Step 4: Conduct continuous tracking simulation of the particulate matter inhalation-deposition-exhalation process: Particles were injected at a rate of 100,000 particles / second at the upper airway inlet, and the status of particles escaping from each outlet of G6 was recorded. The escaped particles are used as the initial condition for the G7 entry point and are then transmitted to the acinar region step by step through a particle data transfer and reconstruction strategy. The particle reconstruction process takes each "parent" particle as the center and randomly generates multiple "child" particles in the surrounding area to maintain a consistent velocity, thereby enhancing the stability of statistical quantities. During the exhalation phase, the particle trajectory is transmitted from the alveolar region upwards. A UDF program is used to map the particles at the truncated exit point, maintaining the relative position and velocity to ensure particle continuity.

[0012] Furthermore, the cartilage ring modeling steps include: locating the central axis of the CT-reconstructed cartilage ring in SpaceClaim, creating a reference plane perpendicular to the tracheal axis, sweeping and excising along the ring path to form an annular groove, applying rounded corner treatment to the edge of the groove, and finally performing geometric repair through the "merging surface" module.

[0013] Furthermore, in the TBU structure, the TBU segments of the same generation are generated using a completely symmetrical regular dichotomy method. All subbronchial bronchi of each level or generation have the same diameter, length and bifurcation angle, and each bifurcation segment is given a random tilt angle to enhance spatial specificity.

[0014] Preferably, the alveolar region modeling includes: generating an initial point cloud in a 45mm³ space, constructing a uniform Voronoi diagram using a point distance optimization algorithm; generating a secondary point cloud based on the polygon centroid, simulating the iterative path using a slime mold growth model, and extracting alveolar channels; and finally performing topology repair and smoothing to complete the alveolar model construction.

[0015] Furthermore, in the segmented linking strategy, five representative outlets with the largest airflow or the most particle escape are selected at the upper airway G6 outlet, and they are connected to the TBU structure respectively through scaling and rotation rules to construct five independent lobe pathways, while the remaining outlets are truncated.

[0016] Preferably, in the particle reconstruction method, the positions of the "child" particles are randomly distributed within a circular region centered on the "parent" particle, and their velocity direction is consistent with that of the "parent" particle.

[0017] Optionally, the exhalation phase particle mapping method includes: mapping particle information at the alveolar outlet or TBU truncated outlet to all corresponding truncated outlet planes, maintaining consistent relative positions and velocities, and performing step-by-step recursion until the particles complete the exhalation process.

[0018] In one possible implementation, a five-branch full airway simulation system for implementing the above method is provided, comprising: a geometric reconstruction module for constructing a three-dimensional upper airway model with cartilaginous rings based on CT data; a TBU modeling module for generating multi-segment triple-branching structures based on a rule-based bisection method and assigning tilt angles; an alveolar construction module for generating Voronoi structure alveoli and flow channels; a dynamic assignment module for defining time-varying motion equations for the glottis, trachea, and alveolar regions; and a data transfer and simulation module for controlling the airflow field and particle segmentation simulation, boundary condition transitions, and particle reconstruction. In one possible implementation, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the above method. In one possible implementation, an electronic device is provided, comprising a processor, a memory, and a computer program stored in the memory, which, when executed by the processor, implements the above method.

[0019] Based on the above technical solutions, the respiratory airflow and particle transport calculation method based on a dynamically coupled full airway model of the present invention significantly improves the performance of the airway model in terms of physiological realism and calculation accuracy by introducing a segmented construction and series-parallel linkage of a five-branch airway structure, combined with dynamic simulation of the time-varying deformation behavior of the glottis, trachea and alveolar regions. Through independent simulation and a step-by-step transfer strategy, it realizes the calculation of complex multi-level airflow and particle transport tracking, effectively reducing the consumption of simulation resources and ensuring the continuity and reliability of particle deposition path and distribution results. It is suitable for various application scenarios such as inhaled drug evaluation, pollutant exposure research and respiratory disease modeling.

[0020] Specifically, the present invention has achieved substantial technical progress in the following aspects: 1. Improve the physiological restoration of airway geometry By using high-resolution clinical CT data and 3D modeling techniques, a realistic upper airway model with cartilage rings was constructed. Combined with the segmented regular branching structure of TBU and the acinar tissue generated by Voronoi, the geometric distortion problem caused by neglecting soft tissue structure in existing CFD models was effectively compensated, and the ability of the model to match individual anatomical differences was enhanced.

[0021] 2. Dynamic motion coupling improves the accuracy of flow field simulation. This invention is the first to integrate the periodic respiratory movements of the glottis, trachea and alveoli into a unified model. It accurately expresses the dynamic deformation process through various time-varying functions (Fourier series, sine function, polynomial), overcoming the problem that traditional static rigid airway models cannot reflect real respiratory dynamic changes and improving the reliability of results such as transient flow field and shear stress.

[0022] 3. The segmented computation strategy significantly reduces simulation resource overhead. By modeling the bronchial segments of G7–G15 as independently solvable triple bifurcation elements (TBUs) and adopting a stepwise transfer of boundary conditions and particle injection method, the number of global meshes and memory consumption are significantly reduced, enabling high-precision simulations of ultra-large-scale models to be performed on ordinary workstations.

[0023] 4. Particle reconstruction and mapping mechanisms ensure sedimentary continuity By employing a "parent-child" particle strategy and a particle mapping method that cuts off the exit point, the trajectory of particles remains continuous and their mass is conserved throughout the entire process of inhalation, deposition, and exhalation. This method is particularly suitable for simulating the transport and deposition behavior of nanomedicines and air pollution particles in multi-stage airways.

[0024] 5. Possesses good versatility and scalability. The model structure supports modular replacement and expansion. For example, the number of TBUs can be adjusted, the acinar density parameter can be changed, and the particle size can be varied to adapt to the individualized airway simulation needs of different ages, genders, and disease states. It has broad scientific research and industrial promotion value.

[0025] In summary, this invention establishes a whole-airway CFD modeling and simulation method that integrates structural anatomy and dynamic physiological characteristics, systematically solving the problem of insufficient coupling in the spatial, temporal, and functional dimensions of existing models. It has high innovation and engineering practicality, and provides reliable technical support for refined respiratory system analysis and simulation. Attached Figure Description

[0026] Figure 1 Flowchart of a method for calculating respiratory airflow and particle transport based on a dynamically coupled whole airway model. Detailed Implementation

[0027] To make the technical solution of this invention clearer and more complete, the following will describe in detail the method for calculating respiratory airflow and particle transport based on a dynamically coupled whole airway model proposed in this invention, in conjunction with the accompanying drawings and specific embodiments. It should be understood that the described embodiments are only for illustrating this invention and are not intended to limit the scope of protection of this invention. Without departing from the concept of this invention, those skilled in the art can make various modifications and substitutions to the specific embodiments based on the content of this document, and these modifications and substitutions should all fall within the scope of protection of this invention.

[0028] I. Structural Modeling Methods for the Upper Airway Segment like Figure 1 As shown, in a preferred embodiment of the present invention, the geometric modeling of the upper airway segment aims to realistically reproduce the airway anatomy from the oral cavity entrance to the 6th generation bronchus (G6), and in particular, considers incorporating cartilage rings as key physiological features into the trachea model to improve the physiological fidelity and engineering applicability of the flow field simulation.

[0029] The specific steps are as follows: 1. CT Data Acquisition and Preprocessing Clinical CT scan data of the human upper respiratory tract were used to ensure that the data resolution met the accuracy requirements of geometric modeling. The acquired DICOM images were imported into the medical modeling software Geomagic Wrap 2021, and the initial point cloud and surface mesh of the upper airway were obtained through image segmentation and thresholding.

[0030] 2. Extraction of key parameters of cartilage rings In Geomagic software, the built-in measurement tools were used to perform quantitative analysis on each group of tracheal cartilage rings, obtaining the following anatomical features: Interring spacing: the average spacing between each cartilage ring along the longitudinal axis; Inclination angle and torsion angle: the distribution of cartilage rings in the axial rotation and tilt directions; Average radius: the cross-sectional dimension of the ring diameter; Opening angle: The notch angle of the C-shaped non-closed ring.

[0031] This parameter system provides a foundation for subsequent personalized modeling and structural positioning.

[0032] 3. Cartilage ring reconstruction and localization modeling The initial mesh was imported into SpaceClaim 2020 R2, and the outer surface of the trachea was reconstructed based on CT.

[0033] At each cartilage ring location, perform the following procedures: Reference coordinate alignment: Move and rotate the local coordinate system until it aligns with the cartilaginous ring axis; Plane generation: Create a reference plane perpendicular to the local tracheal centerline at each ring location; Sweep resection modeling: Perform a sweep resection operation on the plane using an arc path to form a grooved cartilage ring structure with physiological dimensions.

[0034] 4. Boundary fillet treatment and transition optimization To avoid non-physical shear concentration and mesh distortion issues in airflow simulation, a suitable fillet is applied to the edge of each groove to ensure a natural and smooth transition between adjacent ring structures. The fillet radius is adjusted based on anatomical data to ensure that the overall contour is not disrupted.

[0035] 5. Geometric Repair and Model Improvement Since the process of removing cartilage rings may introduce tiny gaps, redundant boundaries, or topological anomalies, the "Merge Faces" function under the "Repair" module in SpaceClaim is used to automatically repair and optimize the model surface, generating a topologically continuous and smooth three-dimensional tracheal structure.

[0036] 6. Confirmation of the overall structure of the airway. The final three-dimensional geometric model of the upper airway segment includes the complete structure from the oral cavity to G6, in which the tracheal segment integrates multiple cartilage ring grooves. The model has high-fidelity anatomical consistency and meets the requirements for mesh continuity and boundary recognition in subsequent CFD simulations.

[0037] Through the above modeling process, the upper airway section not only has a realistic spatial morphology, but also retains key biomechanical characteristics, which helps to improve the simulation accuracy of wall shear, velocity distribution and particle adhesion behavior in airflow simulation.

[0038] II. Modeling Strategy for Triple Bifurcation Units (TBUs) In a preferred embodiment of the present invention, in order to accurately simulate the geometric structure and branching characteristics of the G7 to G15 generation bronchi, a triple bifurcation unit (TBU) modeling strategy is adopted to construct a middle and lower bronchial model with segmentation, controllability, and statistical anatomical features. This strategy can take into account both the physiological laws of the complex branching structure of the lungs and the control requirements of numerical simulation on computational resources.

[0039] The specific modeling process includes the following steps: 1. Segmentation of TBU structure The bronchi from generation G7 to G15 are divided into three TBU units, with each TBU unit corresponding to a three-generation branching structure: TBU1: G7–G9; TBU2: G10–G12; TBU3: G13–G15.

[0040] This segmentation approach helps to manage modeling complexity hierarchically and is compatible with the step-by-step transfer of computational logic.

[0041] 2. Geometric modeling principle: Regular bisection structure Within each TBU unit, a tree-like structure is constructed using regular dichotomy: Each generation of bronchi is strictly divided into two, with the left and right bronchi maintaining the same diameter, length, and bifurcation angle; All branches at the same level have the same forking angle, forming a highly self-similar structure, which facilitates parametric modeling and simulation accuracy control.

[0042] 3. Structural Parameter Definition The average geometric parameters of the bronchi for each generation were determined with reference to data from the classic Weibel airway model, including: Pipe diameter Pipe length bifurcation angle These parameters support the automatic generation of airways in SpaceClaim 2020 R2.

[0043] 4. Introduction of Spatial Tilt Features To enhance the physiological realism of the geometric structure and break the planar repetitive structure, each TBU unit is generated in an independent spatial plane during construction.

[0044] Relative to the axis of the previous TBU, each TBU inlet direction is set with a random tilt angle; The tilt angle is randomly generated between 30° and 65°. This tilt can introduce reasonable spatial distortion features in the lung lobe, effectively improving the realism of particle distribution simulation.

[0045] 5. Construction and Simplification Strategies of the Five-Path Model Considering that the human lung is composed of five main lobes (upper and lower lobes of the left lung, and upper, middle, and lower lobes of the right lung), this model constructs five independent TBU pathways, each representing a lobe region and consisting of a complete G7–G15 chain.

[0046] Preferably, to reduce computing resource consumption, only one exit point is retained at the exit point of each TBU level to connect to the next TBU segment, and the remaining exit points are truncated.

[0047] Each path thus becomes the only passage extending from the G6 exit to the G15; The truncated path will not participate in subsequent CFD simulations to avoid redundant mesh generation; This strategy significantly reduces the computational scale while maintaining the representativeness and structural integrity of the simulation results.

[0048] 6. Connection between TBU structure and upper airway In the five-branch full airway model, five representative outlets with the largest mass flow rate or the highest particle escape rate are selected from the G6 outlet of the upper airway section and used as the inlets of the five TBU paths.

[0049] The selected exit location will become the entrance face of TBU1 (G7–G9); Ensure the continuity of the boundary between airflow simulation and particle simulation; The upper and lower level TBUs achieve seamless connection by scaling the interface size and rotating the coordinate direction.

[0050] Through the TBU modeling strategy described above, this invention constructs a middle and lower bronchial model that achieves a good balance between geometric structure, anatomical consistency and simulation efficiency, providing an efficient and stable modeling foundation for subsequent numerical simulations of airflow transmission, particle migration and other parameters.

[0051] III. Acinar Region Structural Generation and Flow Channel Construction In a preferred embodiment of the present invention, to simulate the microstructure and airflow behavior of the gas exchange region in the deep lung area, an alveolar region structural model is introduced at the terminal end of the G15 generation bronchus to reproduce the physiological functional characteristics of the bronchioles and alveolar tissue. This alveolar region not only serves as the terminal diffusion unit of airflow but also as an important endpoint region for particle transport and deposition calculations.

[0052] The specific modeling method includes the following steps: 1. Establish the modeling space and initialize the point cloud. At the exit of each TBU path at the end of G15, a cuboid 3D modeling region is defined to represent the terminal space of a single lung lobe.

[0053] The modeling space size is set to 45mm³ (length-width-height ratio of 2:2:5), which is set with reference to the anatomical distribution density and direction of alveoli; A certain number of initial point clouds are randomly generated within this space to represent candidate locations of alveolar unit centers; The number of point clouds is sufficient to meet the distribution requirements for the subsequent generation of approximately 500 alveolar structural units.

[0054] 2. Point cloud optimization processing To improve structural uniformity and avoid generating overlapping or overly dense alveolar models, the initial point cloud is optimized: Use minimum spacing constraints to eliminate points that are too close together; Improve the overall uniformity of point cloud distribution through spatial distribution equalization algorithms; The optimized point cloud serves as the basic input for generating the 3D Voronoi structure.

[0055] 3. Construction of polyhedral alveolar structures The three-dimensional Voronoi diagram algorithm is used to generate spatial polyhedral elements representing alveolar structures. Specific operations include: Voronoi units are constructed based on point clouds to form a topologically closed polyhedral space. Each polyhedral structure corresponds to a standard alveolus, with an average equivalent diameter set at 313 μm; The average alveolar path distance was set to 2.5 mm to reflect the actual diffusion path from the end of the airway to the center of the alveoli.

[0056] 4. Construct an internal flow channel network To simulate the propagation path of airflow within the alveolar region and the particle transport channels, the connecting channels between alveoli were further established: The geometric centroid of the polyhedron is calculated using the 3D Thiessen polygon algorithm. A secondary random point cloud is generated in each polyhedron as the nutrient source input; By using a "slime mold growth model" to simulate pathways, the connection paths are made to follow the shortest energy consumption law, thus forming a reasonable connectivity structure. Finally, the reticular branch structure representing the acinar ventilation pathway was extracted.

[0057] 5. Topology optimization and geometric smoothing To meet the requirements of mesh generation and boundary continuity in subsequent CFD simulations, the preliminary modeling results were subjected to topology and morphology optimization: Eliminate the situation where particle flow channels penetrate the polyhedral wall; Automatically merge boundary nodes at the intersection of polyhedral elements; The alveolar surface is smoothed using quadrilateral reconstruction and subdivision modeling techniques to meet the smooth surface boundary conditions required for CFD simulation.

[0058] 6. Structural connection with TBU path After completing the acinar region modeling, connect each group of acinar structures to the unique outlet face of the G15 segment in the corresponding TBU path via a flow channel: Each five-branch path retains only one G15 exit to avoid redundant distribution of multiple paths; The interface is transitioned using interface scaling and spatial orientation correction. This ensures that gas and particles flow smoothly from the TBU segment into the acinar region, achieving continuity of the physical boundary.

[0059] The acinar structures constructed using the above methods possess realistic spatial diversity, microscopic channel distribution characteristics, and anatomical rationality. They not only meet the scientific requirements of anatomical simulation but also the accuracy and controllability requirements of numerical simulation calculations, providing an ideal simulation basis region for simulating particle deposition efficiency, gas diffusion dynamics, and other phenomena.

[0060] IV. Definition of Dynamic Respiratory Boundary Motion Model In a preferred embodiment of the present invention, to realistically reproduce the airway structural deformation characteristics during respiration, a dynamic boundary motion model integrating glottal movement, tracheal deformation, and alveolar expansion is constructed as the time boundary input for CFD numerical simulation. This model controls the changes of airway inner wall or cross-sectional nodes over time through functional expression, reflecting the physical behavior of structural changes with the rhythm of inhalation and exhalation during respiration.

[0061] 1. Dynamic motion modeling of the glottic region The glottis, located at the end of the upper airway, is a key valve-controlled structure regulating the opening and closing of the ventilation channel during inhalation and exhalation; its movement is characterized by periodic changes in its cross-sectional area. This invention uses a Fourier series-based approach to define the displacement function of the node over time: ; in: : Represents the horizontal coordinate of the k-th node on the m-th plane at time t; : This represents the initial position of the node; : A time-dependent Fourier series function that describes the trend of cross-sectional changes within the respiratory cycle; f(y) and F(z) represent the spatial distribution functions of the nodal displacements in the y and z directions, respectively. The ratio of the maximum opening to the static opening determines the maximum opening amplitude.

[0062] This model can simulate the periodic process of the glottis changing from a resting state to its maximum inspiratory opening and then to its expiratory closure, making it suitable for upper airway simulations that require high sensitivity to flow regulation.

[0063] 2. Simulation of triaxial deformation of the tracheal segment To simulate the periodic compression effect of thoracic pressure on the tracheal wall, this invention defines dynamic motion functions of nodes on the airway wall in the range of G6 to G15, and applies wall deformation of different amplitudes in three dimensions by combining the dynamic mesh method.

[0064] The nodal displacement function is defined as follows: ; in: : represents the position of the i-th node at time t; : This is a reference static position; : Coordinates of the tracheal axis; : Deformation rate in each direction; : The duration of a complete respiratory cycle; : This represents the current time step.

[0065] Preferably, the deformation ratio is set to anisotropically, i.e., the x:y:z ratio is 1:0.375:1, to reflect the characteristics of the actual trachea, which is easier to contract laterally and has stronger axial elasticity. This deformation model is loaded through a User Defined Function (UDF) in Fluent to achieve control of wall deformation over time.

[0066] 3. Dynamic expansion model of the acinar region The acinar region is the core area for deep lung ventilation and gas exchange, and its structure exhibits periodic volume expansion and contraction during respiration. This invention establishes corresponding functional expressions for acinar behavior under "spontaneous breathing mode" and "deep breathing mode," respectively: Spontaneous breathing pattern (sinusoidal expansion): ; in: : The path length at any given time; : Path length in resting state; : is the expansion coefficient, set to 0.053; ω is the angular frequency, and T is the respiratory cycle time (e.g., 6s).

[0067] This model can simulate the slow fluctuation of alveolar length over time under relatively regular natural breathing conditions.

[0068] Deep breathing pattern (polynomial + sine combination): When the respiratory waveform no longer exhibits a single sinusoidal shape, such as during deep breathing or disease simulations involving abnormal respiratory amplitude, the following sixth-order polynomial combination model is used to fit the inspiratory phase, while the expiratory phase continues to use a sinusoidal function: ; ; in: , : These represent the duration of inhalation and exhalation (e.g., T=6s, =3s); polynomial coefficient a i Determined based on physiological data or experimental fit.

[0069] Optionally, the acinar motion function can be adaptively adjusted based on the patient's individual lung function test data to achieve specific lung function modeling.

[0070] Through the dynamic boundary control of the above three types of structures, this invention not only realizes the boundary driving of the respiratory airflow field changing with time, but also enhances the trajectory accuracy of particle dynamics under time-varying geometry, forming the core technical foundation for CFD calculation under real physiological rhythm conditions.

[0071] V. Transient Calculation Method of Airflow Based on Segmented Linkage In a preferred embodiment of the present invention, to address the issues of large grid size and high computational resource requirements in CFD simulations of the entire lung airway, a transient airflow calculation method based on segmented linkage is proposed. This method divides the five-branch full airway model into multiple sub-units according to structural hierarchy, and constructs the temporal evolution process of the overall respiratory flow field through segmented independent solutions and step-by-step result transfer.

[0072] This method includes the following key steps: 1. Model Structure Segmentation Principle The overall realistic airway model is divided into three main regions according to the anatomical pathway: Upper airway segment: from the oral inlet to the G6 outlet, including the cartilaginous ring structure; Triple bifurcation unit (TBU) segment: G7–G15, subdivided into 3 TBU units; Alveolar segment: Connected to the G15 exit, representing the end of alveolar gas exchange.

[0073] The segmentation principles are as follows: Each structural element has closed computational boundary conditions; Adjacent segments are connected only through one or more matching interfaces; Each segment can be independently calculated transiently under the condition that the boundary is known, and the output of its outlet velocity distribution result can be used as the inlet of the next segment.

[0074] 2. Transient CFD simulation of the upper airway segment First, perform transient simulation of the entire respiratory cycle for the upper airway segment: The simulated inlet is set to the inspiratory flow rate boundary (e.g., a periodic function fluctuating in the range of 4–12 L / min). Use dynamic boundaries (such as glottal movement) to drive cavity deformation; Output data on the changes in mass flow rate and velocity at outlet G6 over time.

[0075] Preferably, five representative outlets with the largest flow rate or the highest particle escape rate are selected from the G6 outlets to connect to the downstream TBU path.

[0076] 3. Calculation of the first TBU segment (G7–G9) For the five selected TBU paths, perform transient airflow simulations for segments G7–G9 respectively: Load the velocity distribution time series of its corresponding G6 exit at the G7 entrance; Perform a transient solution involving structural deformation to obtain velocity / pressure data at the G9 outlet for each path.

[0077] 4. Calculation of the second TBU segment (G10–G12) and the third TBU segment (G13–G15) The same strategy is used to solve the subsequent TBU segments in segments: The second TBU (G10–G12) loads the G9 exit data at the G10 inlet; The third TBU (G13–G15) loads the G12 exit data at the G13 inlet; Dynamic meshes are used within each segment to accommodate the deformation of the structure over time, ensuring the accuracy of the flow field simulation.

[0078] The calculation results of the exit of each level of TBU segment will be used as the boundary conditions for the entry of the next level segment, and will be passed down level by level until they reach the acinar region.

[0079] 5. Flow field simulation in the acinar region There are no further branching structures within the acinar region; instead, there are closed polyhedral alveolar channel structures generated based on Voronoi notation. G15 outlet serves as the loading boundary velocity distribution for the acinar inlet; The acinar wall performs dynamic boundary control through periodic expansion / contraction; Simulates the diffusion and reflux of gases in the acinar channels.

[0080] Preferably, each alveolar unit is an independent channel structure, which can be set to a zero-pressure boundary or a uniform reflux pressure boundary condition to reflect the open characteristics of the alveolar respiratory end.

[0081] 6. Reverse simulation strategy during the exhalation phase After the complete inhalation cycle simulation is completed, the exhalation phase begins: The velocity / pressure distribution at the end of each segment of inhalation is taken as the initial state; The inlet and outlet are swapped, with the original outlet set as the new inlet, and the speed direction reversed. Segment-by-segment reverse simulation: from the acinar region → G13–G15 → G10–G12 → G7–G9 → G6 exit; Ensure a complete closed loop of the respiratory cycle, forming a continuous evolution of inspiratory and expiratory airflow.

[0082] Optionally, in the exhalation simulation, if it is necessary to simulate the CO2 or drug expulsion diffusion process, a passively labeled gas or a multi-component model can be set in the acinar region to realize the gas exchange analysis function.

[0083] By employing the above strategy of segmented independent calculation and coupling of stepwise boundary conditions, this invention significantly reduces the computational load and memory requirements of the simulation unit, avoids solving the entire domain containing tens of millions of grids at once, and ensures efficient and time-varying accuracy in constructing airflow fields.

[0084] VI. Simulation Mechanism of Particle Inhalation, Deposition and Exhalation In a preferred embodiment of the present invention, a particle transport simulation mechanism based on a Lagrange particle tracking model is constructed to simulate the spatial migration, wall deposition, and partial expulsion of inhaled particles (such as atomized drug particles, air pollutant particles, etc.) during a complex respiratory cycle. This mechanism covers the entire process of particle transport from oral injection, through multi-stage airway transport to acinar deposition, and then to bottom-up movement during expiration.

[0085] Specifically, the following steps are included: 1. Particle injection and tracking during the inhalation phase During the inhalation phase, a transient simulation method is used to set the particle injection boundary at the oral cavity entrance: The particle injection rate is set to 100,000 particles / second, and continuous injection is maintained throughout the entire inhalation cycle. The initial particle velocity is set to be consistent with the gas inlet velocity, and the direction is towards the depth of the gas passage. Physical parameters such as particle size and density can be customized (e.g., diameter 1μm, density 1000kg / m³) to match different research needs.

[0086] After particles enter the airway, their position, velocity, acceleration, and other state parameters are tracked in real time using the Lagrangian method, and their behavior within the upper airway is recorded, including: Whether it deposits upon collision with the wall; Whether it enters the next calculation area through the G6 outlet with the airflow; Key indicators such as arrival time and length of stay.

[0087] Preferably, at the end of each simulation time step, the particles that successfully pass through the G6 exit are extracted and recorded as "parent particles", and their states include position coordinates, velocity vector and departure time.

[0088] 2. Segmented Particle Data Transfer and Sub-Particle Reconstruction Mechanism To ensure the continuity of particle transmission between TBU segments and to enhance particle statistical stability, the "parent particles" are spatially redistributed before entering a TBU segment (such as the G7 entrance) to generate a group of "child particles".

[0089] The reconstruction method is as follows: Each parent particle serves as the center, and multiple child particles are randomly generated within its normal cross-sectional circular domain. The positions of the child particles are randomly distributed within a certain radius centered on the parent particle; The velocity direction of the child particle is consistent with that of the parent particle, maintaining the continuity of momentum and transport trend.

[0090] This method ensures that even if the particle loss rate is high in the first stage, a sufficient number of particles can still be maintained in the next stage for statistical deposition analysis, avoiding simulation bias caused by an insufficient number of particles.

[0091] The reconstruction and transfer process is used sequentially for: G6 → G7 (First TBU segment) G9 → G10 (Second TBU segment) G12 → G13 (Third TBU segment) G15 → Acinar entrance 3. Simulation of particle deposition in the acinar region Within the acinar region, due to the significant reduction in airflow velocity, particles are deposited on the wall primarily through mechanisms such as Brownian diffusion and inertial collisions. Record the number and proportion of particles deposited on the wall of each alveolus; Undeposited particles will continue to move during the exhalation phase; It supports simulating the deposition efficiency variations of particles with multiple sizes and properties.

[0092] Optionally, the particles can be set to inert or charged particles to achieve adaptability modeling for different application scenarios (such as electrostatic atomization).

[0093] 4. Particle backflow tracking and cutoff outlet mapping during the expiratory phase During the exhalation phase, the simulated process reverses, and undeposited particles in the alveolar region begin to move upward with the expiratory airflow: The initial state of the simulation is the position of the particles in each segment at the end of inhalation; The airflow inlet and outlet are set to reverse direction, and the velocity is the negative value of the inhalation velocity; The original exit boundary becomes the entrance, and the new entrance becomes an open boundary.

[0094] For the truncated exit of the TBU segment, this invention innovatively introduces a particle mapping mechanism to maintain path continuity: For example, when one G15 outlet connects to an acinus, the other 7 G15 branches are truncated; Map the particle output results from the G15 outlet onto these 7 cross-sections; During the mapping process, the relative positions and velocities between particles are kept consistent, and only coordinate affine transformations are performed. The same method is applied to the cut-off exits of downstream TBU sections such as G13 and G10.

[0095] This mechanism ensures the physical rationality and quantity conservation of the particle return path from the deep lung region to the upper airway during exhalation.

[0096] 5. Multi-round respiratory cycle control To more closely resemble physiological conditions, this invention supports repeated simulation of multiple respiratory cycles: If, at the end of exhalation, more than 2% of the particles in a certain area are still in a suspended state, the next inhalation-exhalation cycle will begin automatically. The simulation continues iteratively until the particle suspension rate meets the convergence criterion or reaches the preset number of cycles.

[0097] Through the complete particle transport mechanism described above, this invention achieves continuous control of particle state and dynamic response of deposition behavior throughout the entire process of inhalation, deposition, and exhalation, and is applicable to multiple scenarios such as aerosol drug distribution prediction, lung toxin transport modeling, and basic research on particle dynamics.

[0098] VII. Program Implementation and System Deployment Description In an optional embodiment of the present invention, the above-described method for calculating respiratory airflow and particle transport based on a dynamically coupled full airway model can be implemented by a computer program. This program can be deployed on electronic devices including servers, workstations, or dedicated simulation platforms, and run on a computational fluid dynamics (CFD) software environment, such as ANSYS Fluent or OpenFOAM.

[0099] Preferably, the program includes the following functional modules: The geometric reconstruction module is used to read CT data and construct a three-dimensional model of the upper airway with cartilage rings; The TBU modeling module is used to generate multi-level triple-branching structures in segments; The dynamic assignment module is used to control the time-varying deformation of the airway structure in response to the respiratory rhythm. The segmented simulation control module is used to run the airflow and particle transport simulation segment by segment. The particle reconstruction and mapping module is used to handle cross-segment particle transfer and exit mapping operations.

[0100] The program can be stored in a computer-readable storage medium, such as a local hard drive, solid-state drive, cloud storage, or USB flash drive; the above modules can be called and executed by a processor to implement all the steps of the method of the present invention.

[0101] In summary, this invention constructs a five-branch, three-dimensional airway model integrating the upper airway, TBU segment, and alveolar region, and introduces the dynamic motion mechanism of the glottis, trachea, and alveolar region. This enables the realistic evolution of airway geometry with spontaneous breathing rhythm. Furthermore, by combining segmented simulation and particle reconstruction methods, it accurately simulates the entire process of respiratory airflow and particle transport and deposition. This method significantly reduces computational resource requirements while maintaining simulation accuracy, exhibiting good adaptability, scalability, and engineering practicality. It can be widely applied to various scenarios such as inhaled drug delivery, air pollution exposure assessment, and lung disease modeling.

[0102] It should be noted that the above description is only a preferred embodiment of the present invention. Any equivalent transformations or substitutions made by those skilled in the art without departing from the basic concept of the present invention should be covered within the protection scope of the present invention.

Claims

1. A method for calculating respiratory airflow and particle transport based on a dynamically coupled whole airway model, characterized in that, Includes the following steps: Step 1: Construct a three-dimensional model of the five-branched airway, which includes the upper airway segment, the triple bifurcation unit (TBU) segment, and the acinar segment: The upper airway segment is reconstructed based on clinical CT scan data and includes a trachea with a cartilage ring structure. The cartilage ring is formed into annular grooves through sweeping excision and rounding treatment, and surface continuity is optimized using a geometric repair module. The TBU segment consists of three sets of triple bifurcation structures, corresponding to the G7–G9, G10–G12 and G13–G15 generation bronchi, respectively. Each TBU is modeled using the regular dichotomy method, and the geometric parameters are determined based on Weibel morphometry data. Each TBU segment is set with a random tilt angle in the range of 30° to 65°. The acinar segment is generated into a polyhedral alveolar structure using a three-dimensional Voronoi diagram. Each acinar unit contains approximately 500 alveoli, with an average equivalent alveolar diameter of 313 μm and an average path distance of 2.5 mm. The upper airway segment, TBU segment, and acinar segment are connected in a segmented manner through a combination of series and parallel connections to form five independent pathways. Each pathway extends from the G6 exit to the acinar area, corresponding to the five lobes of the lung. Step 2: Establish a dynamic respiratory tract motion model to simulate the periodic deformation behavior of the glottis, trachea, and alveolar regions during inhalation and exhalation. The cross-sectional area change in the glottic region is simulated using a nodal displacement function based on Fourier series, and its displacement is expressed as: ； in: : Represents the horizontal coordinate of the k-th node on the m-th plane at time t; : This represents the initial position of the node; : A time-dependent Fourier series function that describes the trend of cross-sectional changes within the respiratory cycle; f(y) and F(z) represent the spatial distribution functions of the nodal displacements in the y and z directions, respectively. The ratio of the maximum opening to the static opening determines the maximum opening amplitude. The tracheal segment was simulated using a dynamic mesh method to model the triaxial anisotropic deformation of the airway wall, with the nodal displacements being: ； in: : represents the position of the i-th node at time t; : This is a reference static position; : Coordinates of the tracheal axis; : Deformation rate in each direction; : The duration of a complete respiratory cycle; : This represents the current time step; The expansion behavior of the acinar region under spontaneous breathing mode is modeled using a sine function: ； in: : The path length at any given time; : Path length in resting state; : is the expansion coefficient, set to 0.053; ω is the angular frequency, and T is the respiratory cycle time; If it is a deep breathing mode, it is described by a combination of a sixth-order polynomial and a sine function; Step 3: Perform transient airflow calculations using a segmented independent simulation and step-by-step data transfer strategy: First, a transient CFD simulation was performed on the upper airway section to obtain the mass flow rate and velocity distribution at each outlet of G6. Set the TBU segment inlet at five selected representative G6 exits and load their speed data as G7 inlet conditions; The calculation results of each TBU segment are passed down from top to bottom as the entry boundary conditions of the next segment, all the way to the acinar region; During the exhalation phase, the inspiratory airflow field data is reversed according to the time node as the initial condition to simulate the reverse expiratory airflow. Step 4: Conduct continuous tracking simulation of the particulate matter inhalation-deposition-exhalation process: Particles were injected at a rate of 100,000 particles / second at the upper airway inlet, and the status of particles escaping from each outlet of G6 was recorded. The escaped particles are used as the initial condition for the G7 entry point and are then transmitted to the acinar region step by step through a particle data transfer and reconstruction strategy. The particle reconstruction process takes each "parent" particle as the center and randomly generates multiple "child" particles in the surrounding area to maintain a consistent velocity, thereby enhancing the stability of statistical quantities. During the exhalation phase, the particle trajectory is transmitted from the alveolar region upwards. A UDF program is used to map the particles at the truncated exit point, maintaining the relative position and velocity to ensure particle continuity.

2. The method according to claim 1, characterized in that, The cartilage ring modeling steps include: locating the central axis of the CT-reconstructed cartilage ring in SpaceClaim, creating a reference plane perpendicular to the tracheal axis, sweeping and excising along the ring path to form an annular groove, applying rounded corners to the edge of the groove, and finally performing geometric repair through the "merging surface" module.

3. The method according to claim 1, characterized in that, In the TBU structure, each generation of TBU segments is generated using a completely symmetrical regular bisection method. All subbronchial segments of the same generation or level have the same diameter, length and bifurcation angle, and each bifurcation segment is given a random tilt angle to enhance spatial specificity.

4. The method according to claim 1, characterized in that, The alveolar region modeling includes: at 45mm 3 An initial point cloud is generated in space, and a uniform Voronoi diagram is constructed using a point distance optimization algorithm. A secondary point cloud is generated based on the centroid of the polygon. After simulating the iterative path using a slime mold growth model, alveolar channels are extracted. Finally, topology repair and smoothing are performed to complete the acinar model construction.

5. The method according to claim 1, characterized in that, In the segmented linking strategy, five representative outlets with the largest airflow or the most particle escape are selected at the G6 outlet of the upper airway, and they are connected to the TBU structure respectively through scaling and rotation rules to construct five independent lobar pathways, while the remaining outlets are truncated.

6. The method according to claim 1, characterized in that, In the particle reconstruction method, the positions of the "child" particles are randomly distributed within a circular region centered on the "parent" particle, and their velocity direction is consistent with that of the "parent" particle.

7. The method according to claim 1, characterized in that, The exhalation phase particle mapping method includes: mapping particle information at the alveolar outlet or TBU truncated outlet to all corresponding truncated outlet planes, maintaining consistent relative positions and velocities, and performing step-by-step recursion until the particles complete the exhalation process.

8. A five-branch full airway simulation system for implementing the method of any one of claims 1 to 7, characterized in that, include: Geometric Reconstruction Module: Used to construct a three-dimensional upper airway model with cartilage rings based on CT data; TBU modeling module: used to generate multi-segment triple bifurcation structures based on rule-based bisection and assign tilt angles; Alveolar construction module: used to generate Voronoi structure alveoli and flow channels; Dynamic assignment module: used to define the time-varying motion equations of the glottis, trachea, and alveolar regions; Data transfer and simulation module: used for controlling airflow field and particle segment simulation, boundary condition transition and particle reconstruction.

9. A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method of any one of claims 1 to 7.

10. An electronic device comprising a processor, a memory, and a computer program stored in the memory, wherein the processor, when executing the program, implements the method of any one of claims 1 to 7.

Images