JP2015073858A - Aneurysm growth determination method, and device used for the same - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a method for making it possible to accurately determine whether or not an aneurysm grows.SOLUTION: An aneurysm growth determination method is for determining whether or not an aneurysm grows, on the basis of inflow energy E of blood flowing into the aneurysm. In the method, the inflow energy E is calculated on the basis of formula of E=(1/2)×ρ×Qin×|v'in|+Pin×Qin or E=Pin×Qin, where ρ is a density of the blood flowing into the aneurysm, |v'in| is a flow velocity of the blood flowing into the aneurysm in a direction orthogonal to a boundary plane between the aneurysm and a vessel connecting with the aneurysm, Qin is a volume flow rate of the blood flowing into the aneurysm, and Pin is a pressure of the blood flowing into the aneurysm.

Description

  The present invention relates to a method for determining whether an aneurysm grows and a device used in the method.

  One of the conventionally known methods for determining the presence or absence of an aneurysm is a part that is different from the part where the aneurysm is expected and the pulse wave signal of the vascular route through the part where the aneurysm is expected to occur The presence or absence of an aneurysm is determined by detecting the pulse wave signal of the blood vessel path that has passed through and comparing the frequency characteristics of those pulse wave signals (see, for example, Patent Document 1). In addition, several methods for determining the presence or absence of an aneurysm have been proposed (see, for example, Patent Documents 2 and 3).

JP 2013-94264 A International Publication No. WO2009 / 050304 International Publication No. WO00 / 053081

  The above conventional technique is a method for determining the presence or absence of an aneurysm. However, few methods are known for determining whether an aneurysm grows. Furthermore, since various factors are associated with aneurysm growth, it is difficult to accurately determine whether or not the aneurysm grows. Therefore, there is a need for a method that can accurately determine whether or not an aneurysm grows.

The present invention provides a method that can accurately determine whether an aneurysm grows.
More specifically, the aneurysm growth determination method of the present invention is:
As a first step, the inflow energy E of the blood flowing into the aneurysm is determined based on the flow velocity and pressure of the blood flowing into the aneurysm,
As a second step, it is determined whether or not the aneurysm grows based on the inflow energy E.
A method including steps.

In this case, the step of obtaining the inflow energy E is as follows:
The inflow energy E is based on the pressure energy which is “the product Pin · Qin of the volume flow rate Qin obtained based on the flow velocity of blood flowing into the aneurysm and the pressure Pin of blood flowing into the aneurysm”. It may be a step to obtain.

Further, the inflow energy E is, for example,
E = (1/2) · ρ · Qin · | v'in | ² + Pin · Qin
Based on the following formula.
here,
ρ is the density of blood flowing into the aneurysm,
| V′in | is the flow velocity of blood flowing into the aneurysm in the direction perpendicular to the boundary plane between the blood vessel connected to the aneurysm and the aneurysm,
Qin is the volume flow of blood flowing into the aneurysm,
Pin is the pressure of blood flowing into the aneurysm,
It is.

That is, the inflow energy E is expressed as “kinetic energy ((1/2) · ρ · Qin · | v'in | ² ) and pressure energy (Pin · Qin) of blood flowing into the aneurysm through the boundary plane. Can be calculated as However, the inflow energy E may be energy including only the pressure energy (Pin · Qin).

  The inventor investigated the relationship between the inflow energy E and the aneurysm surface area increase for each of a plurality of patients. The results are shown in FIGS. FIG. 9 is a graph that employs the aforementioned “sum of pressure energy and kinetic energy” of blood as the inflow energy E. FIG. 10 is a graph in which the aforementioned “pressure energy” of blood is adopted as the inflow energy E. The surface area increase is an increase in the surface area of the aneurysm per year for each patient.

  From FIG. 9 and FIG. 10, it is understood that the amount of increase in the surface area of the aneurysm increases as the inflow energy E increases. That is, it was found that the larger the inflow energy E, the more the aneurysm grows. Therefore, it can be accurately determined whether or not the aneurysm grows based on the inflow energy E.

  In a more specific example, it is determined whether or not the inflow energy E is equal to or greater than a predetermined threshold value, and it is determined that an aneurysm grows when it is determined that the inflow energy E is equal to or greater than the threshold value. To do. According to this, it can be determined easily and accurately whether or not the aneurysm grows. Of course, the present invention is not limited to such a specific example, and any method can be used as long as it is an aneurysm growth determination method using the fact that the inflow energy E is a parameter strongly correlated with the growth of the aneurysm. Is included in the present invention.

  Furthermore, another aspect of the present invention is an apparatus used for the aneurysm growth determination method. The apparatus treats the aneurysm and the blood vessel connected to the aneurysm as a structure through which “blood as an incompressible fluid” flows, and applies the finite volume method to the blood flowing inside the blood vessel. In this apparatus, the blood flow velocity and pressure are obtained by fluid analysis, and the inflow energy E is obtained using the obtained flow velocity and pressure. According to this apparatus, the inflow energy E can be easily acquired.

  Other objects, other features and attendant advantages of the present invention will be readily understood from the description of the embodiments of the present invention described with reference to the following drawings.

It is the figure which showed typically the three-dimensional image of "aneurysm and its peripheral part" imaged using CT scan (or MRI). It is a figure for demonstrating the boundary conditions for calculating inflow energy. It is a figure for demonstrating the calculation method of inflow energy. It is a graph which shows a secular change of the surface area and inflow energy of an aneurysm of a certain patient. It is a graph which shows the secular change of the surface area and inflow energy of another patient's aneurysm. It is a graph which shows the secular change of the surface area and inflow energy of another patient's aneurysm. It is a graph which shows the secular change of the surface area and inflow energy of another patient's aneurysm. It is a graph which shows the secular change of the surface area and inflow energy of another patient's aneurysm. It is a graph which shows "the relationship between the surface area increase amount of an aneurysm and inflow energy (the sum of pressure energy and kinetic energy)" about each of several patients. It is a graph which shows "the relationship between the surface area increase amount of an aneurysm and inflow energy (pressure energy)" about each of a some patient. It is a perspective view of an element used for fluid analysis concerning an embodiment of the present invention. It is a perspective view of an element used for fluid analysis concerning an embodiment of the present invention. It is a flowchart which shows the process sequence by the SIMPLE method employ | adopted in the fluid analysis which concerns on embodiment of this invention. It is a perspective view of "the element which touches a boundary surface" used for the fluid analysis which concerns on embodiment of this invention. It is sectional drawing of "the element which touches a boundary surface" used for the fluid analysis which concerns on embodiment of this invention. It is a flowchart which shows the process sequence in the fluid analysis which concerns on embodiment of this invention.

  Hereinafter, embodiments of “an aneurysm growth determination method and an apparatus used in the method” according to the present invention will be described with reference to the drawings.

<Image acquisition>
First, a three-dimensional image of “an aneurysm (in this example, a cerebral aneurysm) and its peripheral part” is acquired by CT scan (computer tomography). FIG. 1 schematically shows an example of the acquired three-dimensional image. In this example, the aneurysm 10 is generated at a portion where a plurality of blood vessels 11 to 15 are gathered. The blood vessel 13 is a blood vessel through which blood flowing into the aneurysm 10 flows. The other blood vessels (11, 12, 14 and 15) are blood vessels through which blood flowing out of the aneurysm 10 flows. Note that the shape shown in FIG. 1 is an example, and it goes without saying that the shape of the site where the aneurysm 10 occurs and the blood vessels around the aneurysm 10 are various. The image is not limited to the CT scan, and may be acquired using MRI or the like.

<Acquisition of blood flow rate and pressure>
Next, the “blood flow rate vin and pressure Pin” at the blood inlet 10 a of the aneurysm 10 is acquired. In the present embodiment, these values are obtained by fluid analysis (unsteady fluid analysis) described later. The basics of such analysis methods are well known (for example, “Computational Fluid Dynamics: Versteeg & Malalasekera: Morikita Publishing”, “Computational Fluid Dynamics: Seiichi Koshizuka: Bafukan” and “Computational Fluid Engineering: Taichi Arakawa”). Author: The University of Tokyo Press ”.

  In this analysis, the aneurysm and a blood vessel connected to the aneurysm are treated as a structure through which “blood as an incompressible fluid” flows, and the blood is applied by applying a finite volume method to the blood flowing inside the blood vessel. Obtain the flow rate and pressure. That is, the finite volume method is employed to numerically solve the basic equation of the incompressible fluid. As a mesh used in the finite volume method, a pentahedral element is used near the blood vessel (wall), and a tetrahedral element is used at the center of the blood vessel.

  In addition, as shown in FIG. 2, “pressure as a boundary condition” at each of the end portions 11 a, 12 a, 14 a and 15 a of the blood vessels 11, 12, 14 and 15 through which blood flowing out from the aneurysm 10 flows is “ 0 "is set. Note that the pressure at the end 13a of the blood vessel 13 through which the blood flowing into the aneurysm 10 flows is obtained by calculation, which will be described later.

  Further, the “flow velocity as a boundary condition” at the end portion (inlet end) 13a of the blood vessel 13 through which blood flowing into the aneurysm 10 flows is set to “predetermined flow velocity”. This flow velocity is a vector, and the direction is a direction orthogonal to the end portion 13a (a direction parallel to the central axis of the blood vessel 13).

This “predetermined flow velocity (speed)” is determined as follows.
(1) For various blood vessels having different “thickness and position on the body (for example, the head)”, the flow velocity of blood flowing through the blood vessels is measured in advance and converted into data.
(2) The thickness of the blood vessel 13 and the position of the end portion 13a in the body are specified from the “image of the aneurysm 10” acquired by the CT scan described above.
(3) The flow velocity of the blood flowing in the blood vessel corresponding to the specified thickness and position (having a correlation) is determined based on “the previously measured data”. This determined flow velocity is adopted as the “predetermined flow velocity” described above.

  For example, the graph of FIG. 2 is a blood vessel having a thickness equivalent to that of the blood vessel 13, and measures “the blood flow velocity in one pulsation (pulsation)” for the blood vessel passing through the vicinity of the end portion 13a. Results are shown. Such a flow velocity is given as a boundary condition.

  Further, in this analysis, the flow velocity of the blood surface (ie, the portion in contact with the blood vessel wall) is set to “0”.

  By performing fluid analysis (unsteady fluid analysis) using such boundary conditions, “blood flow velocity (vector) vin and pressure Pin, etc.” at the blood inlet (blood inflow / outflow section) 10a of the aneurysm 10 is obtained. Is acquired. In addition, the blood inflow port 10a means the boundary plane (plane area | region) between the blood vessel connected to the aneurysm 10 and the aneurysm 10, as shown in FIG.1 and FIG.2.

<Acquisition of inflow energy E of blood flowing into the aneurysm>
Next, as shown in FIG. 3, a component (flow velocity surface direct component) v′in (vector) of the flow velocity (vector) vin that is orthogonal to the blood inlet (boundary plane) 10 a of blood flowing into the aneurysm 10. ) Next, the volumetric flow rate Qin of the blood flowing through the blood inlet (boundary plane) 10a and flowing into the aneurysm 10 is obtained based on the following equation (1). In the formula (1), Ain is an area of the boundary plane 10a into which “blood of the flow velocity surface direct component v′in” flows.

Next, the inflow energy E flowing into the aneurysm 10 is obtained based on the following equation (2). In the equation (2), ρ is the blood density. The first term on the right side of the equation (2) is the kinetic energy of blood flowing into the aneurysm 10 through the blood inlet (boundary plane) 10a, and the second term on the right side is the pressure energy of the blood.

  In actual calculation, inflow energy Ei is calculated for each minute surface area along the blood inlet (boundary plane) 10a of each element of the finite volume method, and the inflow energy Ei is the blood in the boundary plane 10a. Is accumulated over the entire region flowing into the aneurysm 10, whereby the inflow energy E is obtained. Further, as described later, the pressure energy (Pin · Qin) of blood flowing into the aneurysm 10 may be employed as the inflow energy E.

<Relationship between inflow energy E and aneurysm growth>
Next, the inventor investigated the relationship between the inflow energy E and the growth of the aneurysm (cerebral aneurysm) 10. The results are shown in FIGS. FIG. 4 to FIG. 8 are graphs in which “points composed of the inflow energy and the surface area of the aneurysm” of the patient in each year are plotted. However, the inflow energy in FIGS. 4 to 8 is the sum of kinetic energy and pressure energy. The surface area of the aneurysm does not include the area of blood inlet (boundary plane) 10a.

  4 and 5 are data on patients (N1, N2) whose inflow energy E is relatively large. For example, the inflow energy E of these patients is 1.5 (N · m / s) or more. For these patients, it was found that the surface area of the aneurysm was gradually increasing (ie, the aneurysm was growing).

  On the other hand, FIGS. 6 to 8 are data on patients (N3, N4, and N5) whose inflow energy E is relatively small. For example, the inflow energy E of these patients is less than 1.5 (N · m / s). For these patients, the aneurysm surface area was found to be substantially constant and hardly increased (ie, the aneurysm was not growing and stable).

  FIG. 9 is a graph showing the relationship between the increase in the surface area of the aneurysm and the inflow energy E. However, the inflow energy E in FIG. 9 is “the sum of pressure energy and kinetic energy” of blood flowing into the aneurysm described above. FIG. 10 is a graph showing the relationship between the increase in the surface area of the aneurysm and the inflow energy E. However, inflow energy E in FIG. 10 is “pressure energy” of blood flowing into the aneurysm described above. In addition, the surface area increase amount in these figures is “an increase amount of the surface area of the aneurysm” per one year (fixed time).

  From FIG. 9 and FIG. 10, it was concluded that the greater the inflow energy E, the greater the amount of increase in the surface area of the aneurysm. Therefore, in this embodiment, it is determined whether or not the aneurysm grows based on the inflow energy E. That is, it is determined that the greater the inflow energy E, the higher the possibility of aneurysm growth. Thereby, the possibility that the aneurysm grows can be accurately estimated.

  In another method, it is determined whether or not the inflow energy E is equal to or greater than a predetermined threshold Eth (for example, Eth = 1.5 (N · m / s)), and the inflow energy E is equal to or greater than the predetermined threshold Eth. Determine that the aneurysm grows. On the other hand, when the inflow energy E is less than the predetermined threshold Eth, it is determined that the aneurysm is stable. This also makes it possible to accurately determine whether or not the aneurysm grows.

<Fluid analysis for obtaining blood flow velocity vin and pressure Pin>
By the way, as described above, the “flow velocity vin and pressure Pin” used for calculating the inflow energy E in the above formulas (1) and (2) are calculated based on the fluid analysis (unsteady fluid analysis). It is calculated by performing it with a computer. This analysis method will be described below. However, the following method is an example, and the present invention is not limited to this.

1. In this analysis, blood is treated as an incompressible fluid. Further, it is assumed that the incompressible fluid satisfies the “continuity equation and momentum conservation equation” as a basic equation. These equations are well known. In the momentum conservation formula, blood viscosity is also taken into account. On the other hand, since it may be considered that there is no temperature distribution, the energy conservation formula is ignored. Furthermore, since the blood flow has a low Reynolds number, it is treated as a laminar flow.

The continuous equation is the following equation (1-1).

The above equations for conservation of momentum are the following equations (1-2), (1-3), and (1-4). The equation of conservation of momentum is also called “Navier-Stokes equation”.

here,
{U} = (u, v, w) is the speed ([m / s])
P is the pressure ([Pa] = [N / m ² ] = [Kg / (m · s ² )]),
ρ is the density ([Kg / m ³ ]),
μ is a viscosity coefficient ([Kg / (m · s)]).

2. Numerical method: Finite volume method 2-1: Theoretical formula of the finite volume method As described above, in order to solve these basic equations, the finite volume method was adopted in this example. For the purpose of ensuring calculation accuracy and automatically creating a mesh, modeling was performed using a pentahedral element (prism element) for the vicinity of the blood vessel and a tetrahedral element (tetra element) for the central part of the blood vessel.

More specifically, the following equation (2-1) is obtained by integrating the above continuous equation (1-1) with a minute space CV inside the fluid region and applying Gauss's divergence theorem. can get. Note that {n} is a normal vector on the surface S of the minute space as shown in FIG.

Similarly, the following equations (2-2 to 4) are obtained from the equations (1-2 to 4) for preserving the momentum.

2-2: Discretization In this example, in order to perform discretization for numerical calculation, the velocity and pressure are assumed to be a collocated grid defined by the center position of each element.

  Further, in this example, the unsteady equation (the following (2-6) equation) in which the pressure and flow velocity for each element are unknown by volume integrating the basic equations for each element (tetrahedral element and pentahedral element). ) And calculate.

More specifically, the flow velocity at the element interface (that is, the flow rate at the element interface) Fs is defined by the following equation (2-5). Here, ΔA _s is the area of the interface, {u _s } is the velocity on the interface, and {n _s } is the normal vector outward from the interface.

Then, the momentum conservation formula (Navier-Stokes formula, formula (2-2-4)) is replaced with the following formula (2-6). Here, ΔV is the volume of the element, P _s is the pressure on the interface, and [τ _s ] is the stress tensor on the interface.

3. Finite volume method for unsteady flow Since blood is treated as an incompressible fluid in this example, the continuous equation (1-1) has no time differential term, but the momentum conservation equation (2-2 ˜4) has a time derivative term. There are explicit methods, a Crank-Nicolson method (central difference method), a frequently used implicit method, and the like as methods for handling the time differential term, and any of them may be adopted. However, in this example, an explicit method that employs easy coding and does not require the solution of simultaneous linear equations is employed.

More specifically, the following equation (3-1) can be obtained by applying the explicit method to the replaced equation (2-6) for conservation of momentum. In this equation, the lower right subscript “P” represents the element of interest, and the upper right subscript “n” represents time. The first term in parentheses in the second term on the right side is a convection term, the second term is a pressure term, and the third term is a diffusion term.

When the explicit method is applied, the time step Δt needs to satisfy the Courant condition expressed by the following equation (3-2). This corresponds to the movement distance of the fluid allowed during the time step Δt being less than 1 mesh. For example, if the maximum speed is 800 [mm / s] and the mesh size is 0.2 [mm], Δt must be 0.25 [ms] or less.

4). Discretization of diffusion term “Difference term in equation (3-1)“ Stress tensor [τ _s ] on interface and vector {γ _PA } connecting the centers between elements sandwiching the interface shown in FIG. The inner product can be calculated by the following equation (4-1), where μ is a viscosity coefficient.

Assuming that the direction of the vector {γ _PA } is equal to the direction of the interface normal vector {n _s }, the diffusion term in the equation (3-1) is expressed by the following equation (4-2). The

The above assumption does not generally hold, but formula (4-2) is adopted for the following reason.
(1) The diffusion term has little influence on the fluid analysis result.
(2) The above assumption holds if the distortion and aspect ratio resulting from the division into meshes are in an appropriate range.
(3) The amount of calculation can be reduced.

By using the equation (4-2), the diffusion term in the equation (3-1) is expressed by the following equation (4-3).

The coefficients of the above equation (4-3) are as follows.

5. Discretization of convection term The convection term in equation (3-1) is given by the following equation.

In order to calculate the flux Fs, it is necessary to calculate the interface center velocity (interface velocity) u _f = {u _s } · {n _s }. In this example, the Rhie and Chow interpolation methods are applied. According to this, the interface velocity _uf between the element P and the element A is given by the following equation.

Here, ∇p | _p is a pressure gradient at the center of the element P, and ∇p | _A is a pressure gradient at the center of the element A. These can be calculated from the pressure of the connecting elements by using a gradient reconstruction (reconstruction) method.

a _p and a _A are values that become diagonal terms related to the speeds of the element P and the element A when the steady-state Navier-Stokes expression is discretized, and the expressions (2-6) and (4-2) ) And the formula. As an example, derivation of _ap is shown below.

この式に（４−２）式を適用すると次式が得られる。

When the equation (2-6) is replaced with a steady state equation, the following equation is obtained.

When the equation (4-2) is applied to this equation, the following equation is obtained.

Further, the following equation is obtained by giving the velocity {u _s } on the interface as the center difference of the elements connected to the interface and multiplying both sides by ρ.

以上により、｛ｕ_ｆ｝が求められるので、（２−５）式から流束Ｆｓを計算することができる。 From this equation, when extracting a coefficient related to the velocity of the element P _{{u p},} the following (5-3) is obtained. The equation (5-3) applies the center difference method instead of the upwind difference method.

As described above, {u _f } is obtained, and thus the flux Fs can be calculated from the equation (2-5).

On the other hand, the upwind difference method, not the center difference method, is applied to the calculation of the velocity {u _s } of the convection term in equation (5-1) in order to improve accuracy. In this example, a TVD (Total Variation Diminishing) method is used as the upwind difference method. In this TVD method, φ _s at the interface between element P and element A is given by the following equation. Here, U represents the leeward element, and D represents the leeward element.

Here, γ _s is defined by the following equation using a vector {γ _PA } between the element P and the element A.

On the other hand, one of a plurality of functions shown below may be adopted as the flow rate limiting function of the formula (5-4).

∇φ _P in the equation (5-5) is obtained by the least square method for gradient reconstruction as described below.

Assuming that the element P is in contact with the elements A, B, C, and D at the interface, the following expression (expression (5-13 to 16)) is obtained. can get.

The least square solution of equation (5-17) is obtained by the following equation. Using the least-squares solution obtained thereby as ∇φ _P.

By the way, when the equation (5-4) is used as it is, the interface velocity appears as a variable. Therefore, the interface speed is replaced with the element center speed by the following procedure.
First, the equation (5-4) is applied to the equation (5-1).

As a result, the equation (5-1) is expressed as the following equation.

The coefficient of the equation (5-22) can be derived as follows. First, when the equation (5-19) is described for each velocity direction component, the following three equations (5-23, 24, 25) are obtained.

When the expression (5-23, 24, 25) is expressed by a matrix, the following expression is obtained.

From this equation, each coefficient of the equation (5-22) is expressed by the following equations (5-27) and (5-28).

6). Discretization of pressure term The pressure term in equation (3-1) is given by the following equation (6-1).

In this example, determined by the pressure section of the next using the pressure P _s on the interface (6-2) below. This expression is an approximate expression.

Here, the pressure P _s on the interface, as shown in the following (6-3) equation, obtained approximately as the average of the pressure of the two elements that share the interface.

From the above, the equation (6-1) is transformed into the following equations (6-4 to 6). These terms determine the pressure term.

7). Update of pressure by SIMPLE method The speed at time n + 1 can be calculated by using the above equation (3-1). However, the pressure itself is calculated as follows using a continuous equation. For calculating the pressure, SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) is applied. The SIMPLE method uses speed and pressure correction values as described below.

First, when the implicit method is applied to the equation (3-1) and further rewritten using the equations (4-3), (5-22) and (6-4), the following (7-1) The formula is obtained.

Here, when the speed predicted by the equation (3-1) is {u _P } ^* , the equation (7-1) is transformed into the following equation (7-2).

Here, it is assumed that there is a relationship of the following equation (7-3) among the true speed u, the predicted speed u ^* , and the speed correction value u ′.
Similarly, it is assumed that there is a relationship of the following equation (7-4) among the true pressure p, the predicted pressure p ^* , and the pressure correction value p ′.

By subtracting (7-2) from (7-1) and applying (7-3) and (7-4), the following (7-5) is obtained.

If the sum regarding the adjacent element velocity in the equation (7-5) is ignored, the following equation is obtained.

The inverse matrix of Expression (7-8) is an inverse matrix of a diagonal matrix (an original matrix that is a diagonal matrix). Therefore, the inverse matrix of the expression (7-8) can be obtained by making the diagonal term of the original matrix the reciprocal (that is, the reciprocal may be multiplied for each vector component).
In addition, the meaning of the symbol which appears afterwards is represented by the following (7-9) Formula.

  The above equation (7-7) indicates that “the speed correction value matches the pressure correction value of the adjacent element”.

Further, an expression relating to the correction value of the pressure is derived by applying the relationship of the expression (7-7) to the continuous expression.
First, equation (2-1), which is a continuous equation for element P, is expressed by the following equation.

Here, as shown in the following equation, it is assumed that the interface center velocity {u _s } ^{n + 1} is expressed by the average of the velocity of adjacent elements.

Here, the following relational expression is established between the surface area of the closed space and the normal vector.

When this is applied to the equation (7-12), the following equation (7-14) is obtained.

When the equation (7-7) is applied to this, the following equation (7-15) is obtained.

If the formula (7-15) is rearranged, the following formula (7-16) is obtained. The coefficient of the expression (7-16) is expressed by the expressions (7-17) and (7-18).

If the equation (7-16) is made simultaneous for each element, a simultaneous linear equation relating to the pressure correction value p ′ is obtained. By solving this, a pressure correction value can be obtained. This simultaneous linear equation is considered to correspond to the solution of the Poisson equation.
Further, by substituting the pressure correction value p ′ obtained by solving the equation (7-16) into the equation (7-7), the element center speed correction value {u _P } ′ can be obtained. it can.

The pressure and speed are updated using the correction value obtained as described above. However, a relaxation method is applied here to prevent the pressure correction value from diverging. When the relaxation coefficient α _p is a value of (0 <α _p <1), the pressure update equation is expressed as the following equation (7-19).

The speed is also updated using the relaxation factor. However, the following equation (7-20) is used.
Here, {u _P } ^old is a value obtained by the previous iterative calculation in the SIMPLE method.

Calculating the flux Fs ^{new new} by applying the rate has been updated by the (7-20) expression (7-21) below. If the continuous equation (7-10) is not satisfied, the calculation of the pressure correction value using the equation (7-16) is repeated.

The calculation procedure of the SIMPLE method described above is shown in the flowchart of FIG. 13 as a reference. Whether or not the flow rate satisfies the continuity formula is determined based on the relative value of the maximum flow rate in consideration of numerical errors (see S4). Note that {u _L } is an inflow velocity.

8). Boundary conditions In this example, the following three types of conditions can be set as boundary conditions.
(1) Time change of speed (2) Time change of pressure (3) Wall surface Although the wall surface is movable, the moving speed of the wall surface is negligible compared to the speed, so it is ignored. Therefore, the coupling with the structure side is realized by satisfying the balance of pressure. For the velocity component at the solid wall, the non-slip condition {u} = 0 is an appropriate condition, and this was also applied to this example.

(1) Diffusion term boundary condition (1A) Velocity boundary condition In “4. Discretization of diffusion term” described above, element P is an element connected to the velocity boundary surface. Then, instead of the element A connected to the element P, as shown in FIG. 14, assuming that the position of the leg of the perpendicular line extending from the center of the element P to the boundary surface is L, the equation (4-2) (8-1). {U _L } is the boundary surface velocity. By applying the equation (8-1) to the equation (3-2), the speed boundary condition can be captured.

(1B) Pressure boundary condition When the pressure boundary condition is set, the boundary surface speed {u _L } is considered to be equal to the element center speed {u _p } ^{n + 1} . Therefore, the following equation (8-2) is obtained as the pressure boundary condition.

(1C) Wall boundary condition When the above-mentioned non-slip condition is applied, {u _L } = 0. Therefore, the equation (8-1) becomes the following equation (8-3).

(2) Boundary condition of convection term (2A) Velocity boundary condition When the velocity boundary velocity is {u _L } and the normal vector of the boundary surface viewed from the velocity boundary element is {n _L }, the velocity boundary interface The speed _uf is expressed by the following equation (8-4) instead of the equation (5-2).

The flux F _{s at the} velocity boundary is calculated using this velocity u _f . Furthermore, the velocity {u _S } ^{n + 1 at the} interface in the equation (3-2) can use the velocity {u _L } at the velocity boundary as it is.

By the way, _ap of the boundary element P is necessary for calculation of the interface velocity with the element connected to the boundary element P and creation of simultaneous linear equations by the SIMPLE method. The _{ap of the} boundary element P is also obtained using the above equation (5-3). However, since the boundary surface velocity {u _S } is known ({u _S } = {u _L }), the equation (5-3) for obtaining _{ap of the} boundary element P is expressed by the following (8-5) ) Formula (see FIG. 15 for the variables of formula (8-5)).

(2B) Pressure boundary condition When the pressure boundary condition is set, it is considered that the boundary surface speed {u _L } is equal to the element center speed {u _P }, so that u _{f in} calculating the flux F _L is Calculate from the following equation (8-6).

The interface center velocity {u _P } ^{n + 1} may be used as it is for the interface velocity {u _S } ^{n + 1} in the equation (3-2).
Furthermore, since the velocity {u _S } = {u _P } at the boundary surface, the expression (5-3) is as shown in the following expression (8-7). F _L is the flux of the interface.

(2C) Wall surface boundary condition When the above-mentioned non-slip condition is applied, {u _L } = 0. Therefore, the equation (8-1) becomes the following equation (8-8).

(3) Boundary condition of pressure term When boundary condition is set (that is, when there are not three or more adjacent elements), a minimum of 2 is necessary for reconstructing the pressure gradient ∇p required in equation (6-1). Multiplication cannot be applied. Therefore, in this case, the pressure gradient is zero (∇p = 0).

(4) Handling of SIMPLE method (4A) Velocity boundary condition If the flux (known) of the velocity boundary is Fin, equations (7-4) to (7-7) are as follows. Where inner
surfaces means interfaces other than velocity boundaries.

Further, it may be assumed that the pressure does not change at the velocity boundary. Thus, the interface pressure _{P S} in calculating the pressure gradient at (6-4) equation is calculated by the following (8-14) below.

(4B) Pressure boundary condition When the pressure boundary condition is set, it is considered that the pressure P _{P at} the element center on the boundary is equal to the pressure P _L on the boundary.
(4C) Wall boundary condition Since there is no inflow / outflow to / from the wall surface, the wall surface is excluded from the sum of equations (7-7). Further, it may be assumed that the pressure does not change at the wall boundary. Therefore, the interface pressure P _S when the pressure gradient is calculated by the equation (6-4) is calculated as P _S = P _P.

<Implementation>
In general, the analysis (simulation) using such a finite volume method can be performed by a general-purpose or dedicated computer having an input unit, a storage unit, a calculation unit, and an output unit as hardware. That is, a model is constructed by dividing a “structure in which fluid (blood) flows” into a plurality of elements using an input unit, the model is stored in a storage unit, and a calculation unit uses the model to create a fluid. Perform analysis. The output unit outputs the analysis result. Note that a preprocessor may be used as a process before the analysis, and a post processor may be used as a process after the analysis.

<Processing procedure>
FIG. 16 is a schematic flowchart showing the flow of processing for performing calculations according to the analysis described above. Each process will be briefly described below.
Step 105: The calculation unit sets initial values of blood velocity (flow velocity) {u} and pressure P. For example, the initial value of the velocity {u} may be set to the maximum value in one pulsation (pulsation) of the “predetermined flow rate” determined as described above. May be set to an average value. The initial value of the pressure P may be set to a predetermined constant value (for example, “0” or “a value other than 0” according to the blood pressure of the patient), and the pressure for each part of the blood vessel can be estimated or acquired in advance. In some cases, the pressure may be set for each part.

Step 110: The computing unit advances time t by Δt.
Step 115: The calculation unit calculates the coefficient of the diffusion term according to the above equations (4-1) to (4-5).
Step 120: The computing unit computes the coefficient of the convection term according to the above equations (5-1) to (5-28).
Step 125: The calculation unit calculates the coefficient of the pressure term according to the above formulas (6-1) to (6-6).
Step 130: The computing unit computes the blood velocity (flow velocity) at time n + 1 (ie, t + Δt) according to the above equations (7-1) to (7-9).

The calculation unit performs a predetermined calculation according to the SIMPLE method described in detail with reference to FIG. Here, each process of FIG. 13 is described more simply.
Step 135: The computing unit creates a coefficient matrix and a right-hand side vector of continuous expressions according to the above expressions (7-10) to (7-18) (particularly, expressions (7-17) and (7-18) See).
Step 140: The calculation unit finds a pressure correction value by solving the simultaneous linear equations of the above expression (7-16).
Step 145: The calculation unit corrects the speed and pressure according to the above equations (7-19) to (7-21).
Step 150: The calculation unit determines whether or not the convergence has occurred (see S2 to S4 in FIG. 13). If the calculation unit has not converged, the calculation unit returns to Step 135. If the calculation unit has converged, the calculation unit proceeds to step 155.

  Step 155: The calculation unit determines whether or not the calculation end time has come. The end time is set to an arbitrary time. If the end time has not come, the calculation unit returns to Step 110. If the end time has arrived, the calculation unit proceeds to step 195 and ends the calculation. As a result of these calculations, the “blood flow rate vin and pressure Pin” at the blood inlet 10a is calculated, and the inflow energy E is calculated.

  As described above, according to the embodiment of the present invention, it can be accurately determined whether or not the aneurysm grows based on the energy flowing into the aneurysm. .

  The present invention is not limited to the above embodiment, and various modifications can be employed within the scope of the present invention. For example, “blood flow velocity (vector) vin and pressure Pin” necessary for obtaining the inflow energy may be calculated using an analysis method different from that of the above embodiment.

10 ... aneurysm 10a ... blood inlet (boundary plane)
11-15 ... Blood vessels

Claims (4)

An aneurysm growth determination method for determining whether an aneurysm grows,
Obtaining the inflow energy E of the blood flowing into the aneurysm based on the flow velocity and pressure of the blood flowing into the aneurysm,
Determining whether the aneurysm grows based on the inflow energy E;
Method. The aneurysm growth determination method according to claim 1,
The step of obtaining the inflow energy E includes:
The inflow energy E is obtained based on the pressure energy that is the product Pin · Qin of the volume flow rate Qin obtained based on the flow velocity of blood flowing into the aneurysm and the pressure Pin of blood flowing into the aneurysm. Is a step,
Method. The aneurysm growth determination method according to claim 2,
The step of obtaining the inflow energy E includes:
Ρ, the density of blood flowing into the aneurysm
| V′in |, the magnitude of the flow velocity of blood flowing into the aneurysm in the direction perpendicular to the boundary plane between the blood vessel connected to the aneurysm and the aneurysm
And when
The inflow energy E is
E = (1/2) · ρ · Qin · | v'in | ² + Pin · Qin
Is a step to obtain based on
Method. An apparatus used for the aneurysm growth determination method according to any one of claims 1 to 3,
Treating the aneurysm and the blood vessel connected to the aneurysm as a structure through which blood as an incompressible fluid flows, and applying the finite volume method to the blood flowing inside the blood vessel, the flow velocity and pressure of the blood are reduced. An apparatus that obtains the inflow energy E using a fluid analysis and uses the obtained flow velocity and pressure.

Images