JP2011040055A - Simulation device and program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a fluid simulation device which easily copes with grid generation and shape changes and reduces computation time while maintaining accuracy. <P>SOLUTION: The simulation device specifies a range corresponding to a viscous bottom layer and generates a prism grid with respect to the specified range corresponding to the viscous bottom layer, generates a tetrahedral grid on the outside of the range corresponding to the viscous bottom layer, defines a composite grid including the generated prism grid and tetrahedral grid, as a grid for computation, and executes numerical fluid computation to compute at least one of a flow velocity in each grid and shearing stress to a flow passage wall surface adjacent to the grid. <P>COPYRIGHT: (C)2011,JPO&INPIT

Description

  The present invention relates to a simulation apparatus that simulates blood flow and the like by numerical fluid dynamics calculation.

  Cardiovascular diseases account for 30% of the deaths of the Japanese people, and are expected to increase further in the future. Information on the blood flow of a patient can evaluate vascular loads such as shear stress and pressure, and is important in diagnosis and treatment. Patent Document 1 discloses an example in which a blood vessel shape extracted from a medical image is reproduced with silicone rubber, and a blood flow is predicted by flowing a glycerin aqueous solution. Patent Document 2 uses a simple voxel-like calculation grid, and corrects the error between the partial flow velocity measurement result and the calculation result by the ultrasonic Doppler method by assuming a virtual body force to improve accuracy. Techniques for making them disclosed are disclosed. However, in practice, there are many cases where the ultrasonic Doppler method cannot be used in the head or in the deep part of the body, and a highly reliable blood flow simulation method is desired.

  Recent advances in computational fluid dynamics have been remarkable due to improvements in numerical computing technology and computer capabilities, and have reached the stage of direct analysis of fine vortices and boundary layers. The range of application is expanding not only to industrial fields such as automobiles, but also to natural phenomena such as weather, as well as biological phenomena related to medicine. Blood flow simulation can evaluate physical quantities related to vascular load that cannot be measured directly, and the development of a diagnosis / treatment support system using this as shown in FIG. 1 will bring a big breakthrough in the diagnosis / treatment of cardiovascular diseases. Expected.

  The vascular network of the human body is a complex structure that repeats the branching and merging of many blood vessels with various diameters, and the grid division corresponding to the discretization in the computational fluid dynamics analysis plays an extremely important role. There are basically two methods for dividing a three-dimensional volume of a blood vessel channel into a grid.

  One is a structured lattice of hexahedral lattices, and it is possible to perform an ordered lattice division in which lattice boundary surfaces responsible for transport between adjacent lattices correspond well to three independent dimensions in three-dimensional space. On the other hand, it is impossible to apply an orderly structure to complicated flow channel shapes at the branching part and the merging part, and it is necessary to use a hexahedral lattice having a small distorted shape. Such determination of the lattice form requires a complicated procedure and a great amount of calculation time. Moreover, it is known that such a calculation grid often causes divergence due to numerical calculation errors and instability.

  The other is an unstructured grid with a tetrahedral grid, which can easily divide a complicatedly shaped flow path with a uniform tetrahedral grid and has a relatively fast fluid calculation. On the other hand, the direction of the boundary surface between adjacent lattices is not orderly, and the distance from the wall surface varies depending on each element, so it is not suitable for resolving a phenomenon with a gradient in the vertical direction of the wall surface like a boundary layer. It can be said that.

  Cardiovascular diseases are often caused by vascular wall lesions such as occlusion of blood vessels due to arteriosclerosis and bleeding due to rupture of aneurysms. The load applied to the blood vessel wall by the blood flow is pressure and shear stress, and it is known that the shear stress particularly affects the cells of the blood vessel wall. The shear stress is dominated by a flow state in a thin layer having a large velocity gradient called a velocity boundary layer in contact with the blood vessel wall.

JP 2005-40299 A JP 2004-121735 A

  Since blood flow is basically laminar, the velocity boundary layer can be analyzed with high accuracy using a hexahedral grid without the need for a special model if there is an appropriate calculation resolution. When a hexahedral element is used for the shape, enormous time is required to create a calculation grid. Although the tetrahedral lattice is easy to create a calculation lattice, it is difficult to accurately analyze the velocity boundary layer from the above-described features even if the same number of calculation lattices as the hexahedral lattice are used. Calculations with a large number of grids with an increased calculation resolution are difficult to execute due to limitations in calculation time and computer capacity.

  On the other hand, since the velocity gradient is gentle at the center of the blood vessel cross section far from the blood vessel wall, there is a situation where numerical diffusion does not become an obstacle even if a tetrahedral lattice is used. In the field of computational engineering, a method is often used in which the vicinity of a wall surface is created by a prism lattice and the other part is created by a tetrahedral lattice when performing flow analysis of a complex shape. However, since the flow is different between the flow analysis and the blood flow analysis in the engineering field, it is difficult to use the calculation grid creation method for the flow analysis in the engineering field as it is for creating the calculation grid for the blood flow analysis. No optimal grid partitioning guidelines have been established.

  The present invention has been made in view of the above circumstances, and provides a fluid simulation apparatus that can easily cope with grid generation and shape change and can reduce calculation time while maintaining accuracy. One of its purposes.

  The present invention for solving the problems of the conventional example described above is a program for creating a calculation grid for causing a computer to obtain a three-dimensional flow of blood inside a blood vessel, and a velocity boundary having a large velocity gradient near the blood vessel wall. The prism extrusion grid is the layer and the range where the influence is strong, the tetrahedral grid is the range where the velocity gradient at the center of the blood vessel cross section is gentle, and the distance between the first grid point of the prism extrusion layer that touches the blood vessel wall is the friction velocity. And the distance corresponding to the dimensionless distance 1 using the kinematic viscosity coefficient are in the range of ± 20%. Here, in the lattice constituting the prism push-out layer, the ratio of the distance between adjacent lattice points after the second lattice point from the blood vessel wall may be between 1.1 and 1.25. The height of the entire grid in the height direction from the blood vessel wall forming the layer may be the minimum number of grids exceeding a distance corresponding to 11.6 in a dimensionless distance using the friction velocity and the kinematic viscosity coefficient. Further, the size of the tetrahedral lattice at the center of the blood vessel cross section may be between 8% and 10% of the blood vessel cross section diameter.

Further, the simulation apparatus according to one aspect of the present invention includes a unit that acquires information for specifying a flow path to be analyzed, and a range corresponding to the viscous bottom layer in the flow path specified by the acquired information. Means for calculating, means for generating a prism grating for a range corresponding to the specified viscous bottom layer, means for generating a tetrahedral grating outside the range corresponding to the viscous bottom layer, and the generated prism grating And calculating a flow velocity in each lattice or a shear stress to a flow path wall adjacent to the lattice by performing a numerical fluid calculation. And display means for displaying information on the calculated flow velocity or shear stress. Here means for generating the prism grating calculates the first grating position y ₀ representing the position of the first grid from the channel wall surface, with the first grating position y _0, the viscous bottom layer height from the flow path wall in the range between _{y max, (a) 1.1 <} y i / y i-1 <1.25 (i = 0,1,2 ...), satisfying the _{(b) Σy i <y max} y _i (i = 1, 2,...) may be determined, and a prism grating having a height of each grating position y _i (i = 0, 1, 2,...) may be generated.

  According to the present invention, it is possible to easily cope with grid generation and shape change, and it is possible to reduce calculation time while maintaining accuracy.

It is a flowchart figure showing an example of processing which a simulation device concerning an embodiment of the invention performs. It is explanatory drawing showing the example of the flow path used as the object of the process of the simulation apparatus which concerns on embodiment of this invention. It is explanatory drawing showing the comparative example of the grating | lattice which the simulation apparatus which concerns on embodiment of this invention produces | generates, and the conventional grating | lattice. It is explanatory drawing showing the comparison with the grating | lattice which the simulation apparatus which concerns on embodiment of this invention produces | generates, and the conventional grating | lattice. It is explanatory drawing showing the example of the analysis result by the grating | lattice which the simulation apparatus which concerns on embodiment of this invention produces | generates, and the conventional grating | lattice. It is explanatory drawing showing the example of the grating | lattice by hexahedral lattice division. It is explanatory drawing showing another example of the flow path analyzed with the simulation apparatus which concerns on embodiment of this invention. It is explanatory drawing showing the example of the composite grating | lattice produced | generated by the simulation apparatus which concerns on embodiment of this invention. It is explanatory drawing showing the example of the analysis result by the grating | lattice which the simulation apparatus which concerns on embodiment of this invention produces | generates, and the conventional grating | lattice. It is a block diagram showing the structural example of the simulation apparatus which concerns on embodiment of this invention. It is explanatory drawing showing another example of the flow path analyzed by the simulation apparatus which concerns on embodiment of this invention. It is explanatory drawing showing the example of the analysis result by the grating | lattice which the simulation apparatus which concerns on embodiment of this invention produces | generates, and the conventional grating | lattice.

[Overview]
In creating a blood vessel shape calculation grid, it is necessary to be able to accurately resolve the wall shear stress (velocity boundary layer near the wall surface) and to accurately resolve the flow velocity distribution in the tube. In addition, in order to reduce the calculation load, it is necessary to keep the number of calculation grids to a minimum. In the analysis of blood flow, the flow in the tube is almost always laminar and the diameter of the blood vessel that can be analyzed is within a certain range. Can be created.

One of the problems to be solved in the present embodiment to be described later is that a prism lattice and a tetrahedron can be achieved in order to achieve both of ensuring the calculation accuracy of the velocity boundary layer governing the shear stress and improving the calculation processing speed. A criterion for determining how much and in what size each grid should be arranged is determined based on known information before blood flow analysis.
At that time, (1) the first lattice point distance of the prism layer, (2) the expansion ratio of the prism layer, (3) the total height of the prism layer, and (4) the length of one side of the tetrahedral lattice are required. Become.

ここで

であり、y_minは壁面からの第１格子点の距離（第１格子点距離）、u^*は壁面摩擦速度、νは動粘性係数である。これら（１）、（２）式と血管内の流れが層流であるという仮定から、血流解析用の計算格子の第１格子点距離の基準を以下のように導く。 When creating a computational grid that can resolve the velocity boundary layer with sufficient accuracy, the first grid point from the wall surface (closest to the wall surface) is such that the dimensionless distance y ⁺ expressed as Taking a grid point) is one criterion.

here

Y _min is the distance of the first lattice point from the wall surface (first lattice point distance), u ^* is the wall surface friction velocity, and ν is the kinematic viscosity coefficient. Based on these equations (1) and (2) and the assumption that the flow in the blood vessel is a laminar flow, the reference for the first grid point distance of the calculation grid for blood flow analysis is derived as follows.

ここで

である。y⁺=1とおき、血液の平均的な密度、粘性、及び最大流速を（３），（４）式に与えると、第１格子点の距離y_minが見積もられる。見積もられた値を中心にプラスマイナス２０％程度の範囲を最適範囲とする。この最適範囲を超えてy_minを大きく取り過ぎると壁面近傍での境界層を十分な精度で解析出来なくなる。また、最適範囲を下回ってy_minを小さく取り過ぎると、プリズム要素を構成する３角形の辺とプリズム格子高さのアスペクト比が大きくなり過ぎるため、壁面近傍の計算格子の品質が低下する。 That is, when the relationship between the blood vessel radius and the wall friction velocity is derived from the equation of the flow velocity of the Hagen-Poiseuille flow, the distance y _min of the first grid point is as follows.

here

It is. When y ⁺ = 1 is set and the average density, viscosity, and maximum flow rate of blood are given in the equations (3) and (4), the distance y _min of the first lattice point can be estimated. A range of about plus or minus 20% around the estimated value is set as the optimum range. If y _min is too large beyond this optimum range, the boundary layer near the wall cannot be analyzed with sufficient accuracy. On the other hand, if y _min is made too small below the optimum range, the aspect ratio between the triangle side constituting the prism element and the prism grating height becomes too large, so that the quality of the calculation grating near the wall surface deteriorates.

  Generally, it is desirable that the reference value of the adjacent lattice spacing ratio in the vertical direction of the wall be within about 1.2. When creating a calculation grid for blood flow analysis, the expansion ratio range is 1.1 to 1.25, and the expansion ratio is adjusted within this range to satisfy the distance of the first grid point and the overall height of the prism layer described later. To create a prism layer.

In order to resolve the boundary layer with sufficient accuracy, the total height of the prism layer is determined so that the prism layer can be disposed in a region called a viscous bottom layer in the boundary layer near the wall surface where the velocity gradient is large. The viscous bottom layer is generally defined as y ⁺ <11.6. If y ⁺ = 11.6 in the equation (3), the thickness of the viscous bottom layer can be obtained. Like the Willis artery ring, which is a site where an aneurysm is common, the blood vessel shape to be analyzed often has a blood vessel diameter of about 2 mm to 20 mm and an inflow / outflow diameter ratio of 3 or less. As a result of trial manufacture considering the connection with the tetrahedral lattice in the central part in addition to these things, when the total prism layer height is about 25% (± 5%) of the minimum cross-sectional diameter of the blood vessel shape Therefore, it is possible to create a calculation grid that has good overall quality and can efficiently resolve the viscous bottom layer over the entire blood vessel.

  The length of one side of the tetrahedron lattice is about 10% of the diameter in the inflow / outflow cross section of the blood vessel shape as a result of trial manufacture in consideration of the connection with the prism layer described above and the overall balance. On the wall surface, a calculation grid having a good overall balance could be created when the maximum cross-sectional diameter was about 8% to 10%. If the length of one side of the tetrahedral lattice is made too small, the calculation result obtained is almost the same as the result of the optimum value, although the calculation load becomes enormous. If it is too large, the flow in the blood vessel cannot be resolved with sufficient accuracy.

  In this way, in the blood flow simulation, the type and arrangement of grids can be allocated according to the state of blood flow in the blood vessel cross section and the degree of numerical diffusion failure, and the time required for grid formation and calculation execution can be saved. Moreover, since the type, arrangement, and dimensions of the calculation grid are numerically designated, it is easy to automate the calculation grid generation program from the medical image.

[hardware]
A simulation apparatus according to an embodiment of the present invention will be described with reference to the drawings. Note that the simulation device of the present embodiment is applicable not only to simulating blood flow (blood flow) in blood vessels but also to all fluids that flow through channels defined in a three-dimensional space. In the following example, a case where blood flow is simulated will be described as an example.

FIG. 10 is a configuration block diagram illustrating an example of the simulation apparatus 1 according to the present embodiment. As illustrated in FIG. 10, the simulation apparatus 1 according to the present embodiment includes a control unit 11, a storage unit 12, an operation unit 13, and a display unit 14.
Here, the control unit 11 is a program control device such as a CPU, and operates according to a program stored in the storage unit 12. The control unit 11 of the present embodiment first acquires information for specifying the fluid flow path to be simulated. Then, a lattice for numerical fluid calculation is generated in the flow path. The control unit 11 executes the numerical fluid calculation using this lattice and outputs the calculation result. One characteristic of the present embodiment is that the grating generated by the control unit 11 is a combination of a prism grating and a tetrahedral grating. Here, the prism grating is a grating having a substantially triangular prism shape, and depending on the shape of the flow path, in order to arrange the grating without gaps, one of the opposing triangular surfaces is tapered so as to be smaller than the other. Put on. Details of the processing of the control unit 11 will be described later.

  The storage unit 12 is realized by a memory device, a disk device, or the like. The storage unit 12 holds a program executed by the control unit 11. In the present embodiment, the program may be provided by being stored in a computer-readable recording medium such as a DVD-ROM and installed in the storage unit 12. The storage unit 12 also operates as a work memory for the control unit 11.

  The operation unit 13 is a keyboard, a mouse, or the like, and accepts a user's instruction operation and outputs it to the control unit 11. The display unit 14 is a display or the like, and displays and outputs information according to an instruction input from the control unit 11. The simulation apparatus 1 may include a network interface. When a network interface is provided, the simulation apparatus 1 may acquire information specifying the fluid flow path via a network.

[Processing of control unit 11]
Here, specific processing by the control unit 11 will be described. As illustrated in FIG. 1, the control unit 11 accepts input of a CT image (Computed Tomography image) or an MRI (Magnetic Resonance Imager) image (S1). Then, the control unit 11 analyzes these images and extracts three-dimensional flow path information (blood vessel shape information) (S2). Since a method for obtaining three-dimensional blood vessel shape information based on CT images and MRI images is widely known, detailed description thereof is omitted here. Here, the three-dimensional flow path information may be expressed as, for example, a flow path represented by a mesh and represented as a three-dimensional coordinate of a lattice point of the mesh.

  The controller 11 refers to the extracted flow path information and generates a grid for numerical fluid calculation in the flow path (S3). In the immediate vicinity of the inner wall surface of the channel (blood vessel wall, hereinafter referred to as the wall surface), a viscous bottom layer in which the molecular viscosity action dominates occurs, and the outside of the viscous bottom layer (away from the wall surface from the thickness of the viscous bottom layer) On the center side), there is a layer where the flow velocity follows the smooth logarithmic distribution law and the velocity gradient is gentle.

  The control unit 11 according to the present embodiment performs calculation in which a triangular prism-shaped prism lattice is arranged in the viscous bottom layer without gaps, and a tetrahedral lattice is arranged without gaps in a range closer to the center of the blood vessel than the viscous bottom layer. Generate a grid for. An outline of this lattice is shown in FIG. That is, when the tube surface is viewed in plan, the shape is such that triangles are arranged (FIG. 8 (a)), and when the tube surface is cut by a plane perpendicular to the flow path, the viscous bottom layer has a substantially triangular prism shape. As seen from the side, a square shape is arranged. In addition, on the inner side (vessel center side) from the viscous bottom layer, a state in which tetrahedron shapes are arranged is seen (FIG. 8B).

In each of these layers, the formulas representing the flow velocities are different, but it is a matter of course that the velocities of the two coincide at the boundary. Therefore, an equation is generated on the assumption that these flow velocity equations are equal to each other, and is solved for the dimensionless distance y ⁺ in the equation (1), a value y ⁺ = 11.6 is obtained. That is, the thickness of the viscous bottom layer can be expressed as y ⁺ = 11.6.

  A widely known method can be adopted as a method for generating the lattice by the control unit 11, and an example thereof will be described below. The control unit 11 arranges one bottom surface (a larger triangular surface in the case of applying a taper) of the plurality of prism gratings 21 with no gap so as to contact the wall surface 20. As described above, a widely known method can be adopted as a method of arranging triangles on an arbitrary surface without a gap.

Here, the control unit 11 calculates the first grid point distance y _min as follows. In the following, it is assumed that the local blood flow in the lattice is a Hagen-Poiseuile flow. One of the vertices of a triangle arranged on the wall surface 20 is regarded as an attention point, and an average radius of the inner surface (wall surface) of the blood vessel in a surface perpendicular to the axis of the blood vessel (flow direction) is defined as a. In addition, when the cross section of the blood vessel on this surface is approximated by a circle, the distance from the center of the circle (the radius of this circle is the average radius) is r, and the wall friction velocity u ^* is

で定義する。

Define in.

  Note that ν is a kinematic viscosity coefficient of the target fluid (here, blood). Then, the control unit 11 determines the distance of the first lattice point from the wall surface,

の無次元距離y⁺を１とするようなy_minとして定める。

_Is defined as y _min such that the dimensionless distance y ⁺ is 1.

Specifically, using the equation of the flow velocity of the Hagen-Poiseuille flow and deriving the relationship between the blood vessel radius and the wall friction velocity, the distance y _min of the first lattice point is as follows from the equation (1).

ここで

である。制御部１１は、この（３），（４）式を用い、y⁺=1とおいて、血液の平均的な密度、粘性、及び最大流速（既知である）を代入して、第１格子点の距離y_minを得る。

here

It is. The control unit 11 uses the equations (3) and (4), sets y ⁺ = 1, substitutes the average density, viscosity, and maximum flow rate (known) of the blood to obtain the first lattice point. Get the distance y _min .

Here, the control unit 11 performs the following process using the distance y _min of the first grid point calculated in this way, but experimentally, this value is set to plus or minus 20%. It is known that it may be within a range.

That is, due to other factors, the control unit 11 sets the distance y _min of the first grid point calculated here as a temporary distance, and based on this temporary distance, the distance y ′ of the first grid point actually used for the calculation. You may calculate _min . At this time, y ′ _min / y _min may be between 0.8 and 1.2.

As already described, if y _min is set too large beyond the optimum range, the boundary layer near the wall surface cannot be analyzed with sufficient accuracy. On the other hand, if y _min is made too small below the optimum range, the aspect ratio between the triangle side constituting the prism element and the prism grating height becomes too large, so that the quality of the calculation grating near the wall surface deteriorates.

Next, the control part 11 determines the range (distance from a wall surface) which arrange | positions a prism grating. As already described, since the viscous bottom layer is in the range of y ⁺ <11.6, the control unit 11 corresponds to (3).

を用いてy⁺=11.6を代入し、粘性底層の厚さy_maxを求める。

Substituting y ⁺ = 11.6 using, find the thickness y _max of the viscous bottom layer.

The control unit 11 determines the first lattice point based on a predetermined adjacent lattice spacing ratio (empirically 1.1 or more and 1.25 or less when blood is used as in this example). The distance y ₀ = y _min and the distances y ₁ , y ₂ ,... To the following grid points (second, third,.
(A) 1.1 <y _i / y _i-1 <1.25 (i = 0, 1, 2,...),
(B) Σy _i <y _max
It is determined to satisfy each of the conditions. Here, in (b), i is summed up.

For example, the distance y ₀ = y _{min of} the first lattice point is set, and y ₁ , y ₂ ,... Are sequentially repeated while satisfying the condition (b) so as to satisfy the condition (a). Just ask for it.

Note that the ratio of y _i / y _i-1 may be obtained as a value at which Σy _i is closest to y _max . Since there is a widely known method for obtaining such a value by a CPU such as the control unit 11, detailed description thereof is omitted here.

The thickness of the viscous bottom layer is also affected by the shape of the blood vessel, the diameter of the blood vessel, the diameter ratio of the inflow / outflow part, and the like. Therefore, the range (distance from the wall surface) in which the prism grating is arranged is found in advance as the diameter r _min where the cross-sectional diameter is minimum in the blood vessel shape obtained in step S2, and 25% (± 5) of the r _min The dimensionless distance y ⁺ corresponding to%) may be set to _ymax .

The control unit 11 creates virtual wall surfaces from the wall surfaces at positions i _i (i = 0, 1, 2,...) Of the i-th lattice point obtained as described above. The control part 11 divides | segments a wall surface into the mesh which consists of a mutually adjacent triangle (triangular division). Then, vertex coordinates T ₁ , T ₂ ... Of each triangle obtained by this triangulation are obtained.

Next, the control unit 11 also performs triangulation on the virtual wall surface generated at the distance of the first grid point by the same method, and apex coordinates T ₁ and T _{2 of} each triangle obtained by triangulating the wall surface. ... vertex coordinates of each triangle corresponding to T _{0_1,} T _{0_2} ... obtained. Similarly, on the virtual wall surface generated at the position of the distance between the second, third,... Lattice points (i-th lattice point), the vertex coordinates T ₁ , T _{3 of} each triangle obtained by triangulating the wall surface are similarly obtained. The vertex coordinates T _{i_1} , T _{i_2} ... Of each triangle corresponding to T ₂ .

  Then, corresponding vertex coordinates are connected to each other on the wall surface and each virtual wall surface, and a prism layer made of a substantially triangular prism lattice is formed densely.

  The control unit 11 forms a tetrahedral lattice 22 including a tetrahedron having a triangle formed by triangulation as one of the faces on a virtual wall closest to the center of the blood vessel, and forms a prism layer on the center side of the blood vessel. Distribute densely in areas not covered. Since this method can employ a well-known division process into a tetrahedral lattice, a detailed description thereof is omitted here.

  Note that the control unit 11 may determine the size of one side of the triangle generated on the wall surface 20 when performing the triangulation on the wall surface 20 that is the inner surface of the blood vessel as follows. In other words, the length of one side of the tetrahedral lattice 22 is about 10% of the diameter in the inflow / outflow cross section of the blood vessel shape by an experiment conducted while considering the connection with the prism layer described above and the overall balance. In addition, it has been found that it is appropriate to set the maximum cross-sectional diameter to about 8% to 10% on the blood vessel wall surface. This is because if the length of one side of the tetrahedral lattice 22 is made smaller than the above level, the calculation accuracy cannot be improved in spite of the huge calculation load. Based on the fact that the flow in the blood vessel cannot be resolved with sufficient accuracy.

  Therefore, for example, if the length of one side of the tetrahedral lattice 22 is set to 10% of the cross section of the blood vessel, the length of one side of the triangle formed by triangulation on the virtual wall surface closest to the center of the blood vessel is also obtained. Since the blood vessel cross-section is 10%, the wall surface (based on the distance from the center of the virtual wall surface closest to the center of the blood vessel and the distance from the center to the wall surface (inner surface of the blood vessel) The size of one side of the triangle generated on the inner surface 20 of the blood vessel can be determined.

  As already described, the control unit 11 uses other methods to arrange the triangular prism prisms 21 in the viscous bottom layer without any gaps, and the area on the center side of the blood vessel from the viscous bottom layer is divided into four surfaces. A grid for calculation in which the body grid 22 is arranged without a gap may be generated.

  When the control unit 11 generates a calculation grid for numerical fluid calculation by the above processing, the control unit 11 executes a fluid simulation process using the calculation grid (S4). For this simulation process, since a well-known numerical dropout calculation method can be used, a detailed description is omitted, but with this simulation process, the time change of the flow velocity of the fluid in each calculation grid is calculated.

  When the controller 11 obtains the time change of the fluid flow velocity in each calculation grid, it calculates the stress on the wall surface (blood vessel inner surface) adjacent to each calculation grid (S5). Then, the control unit 11 generates the three-dimensional graphics generated based on the three-dimensional flow path information obtained in the process S2 with respect to the display unit 14 (the generation of this graphic is also widely known, and thus detailed In the three-dimensional graphics, a predetermined color information is obtained in association with the stress value calculated for each wall surface adjacent to the calculation grid. The wall portion is colored in the acquired color (S6). An example of the image thus obtained is shown in FIG.

  The user inputs an instruction from the operation unit 13 and performs an operation such as confirming a position where stress is concentrated while rotating the image three-dimensionally.

[Example 1]
FIG. 2 shows a simulation of a bifurcation that is a basic blood vessel form. The case where the T-shaped branch pipe shown in FIG. 2 is divided into a hexahedral lattice, a case where the T-shaped branch tube is divided into a tetrahedral lattice, and a combination of a prism lattice and a tetrahedral lattice by the control unit 11 of the present embodiment (hereinafter referred to as a “tetrahedral lattice”) FIG. 3 shows respective processing times (computation grid creation times) when divided into composite grids). Here, the resolution is substantially the same. As a result, the generation time of the composite lattice according to the present embodiment is one hundredth of that of the hexahedral lattice, and is almost the same as that of the tetrahedral lattice.

  Further, FIG. 4 shows a comparison of calculation times when a numerical fluid calculation is performed using these calculation grids. As will be described later, in the composite lattice of the present embodiment, the time step can be increased (the number of steps can be reduced). Therefore, even though the number of lattices is increased and the resolution is increased compared to the hexahedral lattice, the calculation time is hexahedral. It was saved to 1/10 or less of the lattice. In FIG. 4, an example in which the time step is not increased (the number of steps is decreased) in the composite grid is also included for comparison.

  FIG. 5 compares the calculation results of the velocity distribution near the wall surface that is the basis for the shear stress evaluation. The composite lattice of this embodiment was able to obtain a boundary layer velocity distribution that was in good agreement with the high-precision calculation using a hexahedral lattice. In the hexahedral lattice division, a small lattice having a large distortion is generated at the branch pipe connecting portion as shown in FIG. 6, and the time step is likely to be restricted.

  Specifically, the Courant number C (flow rate u, time step width Δt, lattice width Δx, C = uΔt / Δx) of the hexahedral lattice is 2.298, whereas the composite lattice according to the present embodiment The Courant number C of 1.442 is 1.442, and if the flow velocity is the same, Δt / Δx can be made smaller than that of the hexahedral lattice.

[Example 2]
As yet another embodiment, FIG. 8 shows analysis conditions for blood flow simulation using the carotid artery shape obtained from the medical image shown in FIG. Compared with the conventional tetrahedral lattice, the composite lattice of the present embodiment was able to calculate a smooth and reasonable shear stress distribution as shown in FIG. Since tetrahedral lattices are not ordered lattices, distribution calculations tend to be uneven.

In the calculation of FIG. 8, as a discretization method,
Space: Finite volume method (convection term: second-order accuracy center difference)
Time: Completely implicit method, time step size Δ t = 1.0 × 10 ^-3 seconds, boundary condition as
Inflow: 0.4m / s (uniform),
Wall: no slip, rigid wall,
Outflow: free outflow boundary condition,
It was. The calculation grid (analysis grid) is the composite grid of the present embodiment. At this time, the number of elements was 225,114.

  In the hexahedral lattice and the composite lattice of the present embodiment, since the lattice is orthogonal to the wall surface vertical direction, the continuity of the shear stress applied to the wall surface is easily maintained.

[Example 3]
Next, when the flow path shape is the middle cerebral artery (MCA: see FIG. 11), the inflow condition is Poiseuille flow with an average flow velocity of 0.28 m / s, the outflow condition is free outflow, and the Reynolds number is 575. The results of the numerical fluid calculation will be described.

  Here, the calculation scheme of the numerical fluid calculation employs a second-order center difference, and the crank-Nicholson method is employed for time integration. Moreover, when outputting the results, those averaged at intervals of 0.02 seconds were used.

  The results of the numerical fluid calculation performed on each of the composite lattice, hexahedral lattice, and tetrahedral lattice of the present embodiment are as shown in FIGS. 12 (a) to 12 (c). FIG. 12A shows the result of using the composite lattice of the present embodiment, and FIG. 12B shows the result of using a hexahedral lattice. The result of using the composite lattice of the present embodiment and the result of using the hexahedral lattice (display result of the distribution of shear stress) are almost the same. On the other hand, in the result using the tetrahedral lattice shown in FIG. 12C, there is no continuity in the distribution of shear stress, unevenness occurs, and the stress concentration state does not appear clearly at the position where the stress should be concentrated. Such a phenomenon was seen.

  Further, the calculation time (hour: minute) when the calculation step (time step Δt) is changed to 1.0d-4, 5.0d-4, 1.0d-3 is shown in the following table. Show.

  The hexahedral lattice diverged without converging the calculation result when the time step width Δt = 1.0 d−3. In addition, regarding the tetrahedral lattice having a relatively coarse Δx (coarse), the flow could not be reproduced, and an appropriate result could not be obtained.

  As shown in this table, in the composite lattice of the present embodiment, even when the time interval Δt = 1.0d−3 is not diverged, the calculation is completed in 28 hours and 18 minutes, and more than in the other lattices. Reasonable results can be obtained in a short time.

  In addition, the time for generating the calculation grid is compared as shown in the following table.

  As described above, in the simulation apparatus according to the present embodiment, the portion corresponding to the viscous bottom layer is a substantially triangular prism-shaped prism lattice, and the portion corresponding to a layer having a gentle velocity gradient is a tetrahedral lattice. Since a composite grid is generated and used for computational fluid calculations, the grid creation time and calculation time are relatively close to a tetrahedral grid, making it easier to handle grid generation and shape change, compared to a hexahedral grid. The time is reduced, and the result of the analysis is close to a hexahedral lattice, and continuous reasonable results are obtained. That is, the calculation time can be shortened while maintaining accuracy.

  Note that the simulation apparatus of the present embodiment can be applied not only to diagnosis and treatment of circulatory system diseases, design of medical equipment, etc., but also to general fluid calculations.

DESCRIPTION OF SYMBOLS 1 Simulation apparatus, 11 Control part, 12 Storage part, 13 Operation part, 14 Display part, 20 Wall surface, 21 Prism grating | lattice, 22 Tetrahedral grating | lattice.

Claims (6)

  This is a program that creates a calculation grid that allows a computer to determine the three-dimensional flow of blood inside a blood vessel. The velocity boundary layer near the blood vessel wall has a large velocity boundary layer and the range of its influence is a prism extrusion grid. The range in which the velocity gradient at the center of the cross section is gentle is a tetrahedral lattice, and the distance of the first lattice point in contact with the blood vessel wall of the prism extruded layer corresponds to the dimensionless distance 1 using the friction velocity and the kinematic viscosity coefficient. The program is characterized in that the difference between is in the range of ± 20%.   2. The program according to claim 1, wherein the ratio of the distance between adjacent lattice points after the second lattice point from the blood vessel wall is between 1.1 and 1.25 in the lattice constituting the prism push layer.   The height of the entire lattice in the height direction from the blood vessel wall forming the prism pushing layer is a minimum number of lattices exceeding a distance corresponding to 11.6 in a dimensionless distance using the friction velocity and the kinematic viscosity coefficient. The program according to claim 2.   The program according to claim 3, wherein the size of the tetrahedral lattice at the center of the blood vessel cross section is between 8% and 10% of the blood vessel cross section diameter. Means for acquiring information for identifying a flow path to be analyzed;
Means for specifying the range corresponding to the viscous bottom layer by calculation in the flow path specified by the acquired information;
Means for generating a prism grating for a range corresponding to the identified viscous bottom layer;
Outside the range corresponding to the viscous bottom layer, means for generating a tetrahedral lattice;
A composite lattice including the generated prism lattice and tetrahedral lattice is used as a lattice for calculation, and numerical fluid calculation is performed to perform at least flow velocity in each lattice or shear stress on the channel wall surface adjacent to the lattice. Means for computing one;
Display means for displaying information of the calculated flow velocity or shear stress;
Including a simulation apparatus. The simulation device according to claim 5,
The means for generating the prism grating calculates a first grating position y0 representing the position of the first grating from the channel wall surface, and between the first grating position y0 and the viscous bottom layer height ymax from the channel wall surface. In the range of
(A) 1.1 <yi / yi-1 <1.25 (i = 0, 1, 2,...)
(B) Σyi <ymax
Y i (i = 1, 2,...) That satisfies the above conditions is determined, and a prism grating having a height of each grating position y i (i = 0, 1, 2,...) Is generated.

Images