Method of, apparatus for, and program for fluid simulation

Abstract

To carry out easily fluid simulation in a method of, an apparatus for, and a program for fluid simulation, by using a fluid basic equation embedded with a solid boundary.SOLUTION: A method for fluid simulation according to the present invention has the steps of acquiring shape data for a solid state part, setting an optional calculation grid for the shape data, defining a porosity rate corresponding to the shape of the solid state part as an index of smooth fluidness for the calculation grid, and carrying out fluid simulation based on modified Navier-Stokes equation including the porosity rate.SELECTED DRAWING: Figure 2

Description

[Technical Field]

[0001] The present invention relates to a method, an apparatus, and a program for fluid simulation. [Background technology]

[0002] Patent Document 1 discloses a fluid simulation device that uses a three-dimensional computational grid to numerically analyze the behavior of a fluid flowing around an object (solid part). When analyzing this fluid simulation device, it is necessary to create a complex computational grid that matches the solid boundary (interface).

[0003] Non-Patent Document 1 discloses the basic equations for flow with an embedded solid boundary. These equations are proposed as new basic equations for flow in which the porosity (fluid occupancy) corresponding to the solid shape is introduced and the solid boundary is embedded, in contrast to the conventional Navier-Stokes equations which treat the solid shape that forms the boundary of the flow as a boundary condition. [Prior art documents] [Patent documents]

[0004] [Patent Document 1] Japanese Patent Application Publication No. 11-66356 [Non-patent literature]

[0005] [Non-Patent Document 1] Nobuyuki Oshima, "Basic Equations for Flows with Embedded Solid Boundaries," 36th Fluid Mechanics Symposium, December 2022 Summary of the Invention [Problem to be solved by the invention]

[0006] The fluid simulation device of Patent Document 1 requires the creation of a complex computational grid that matches the interface, which requires a large number of steps to implement the simulation. In fact, the creation of a computational grid that matches the interface is a bottleneck process in implementing the simulation, and there is a demand for reducing this process.

[0007] Non-Patent Document 1 does not disclose in detail how to specifically implement the above-mentioned new fundamental equations in a fluid simulation. For example, there is room for improvement in terms of a method for simply defining porosity from shape data of an object (solid part) and in formulating the viscous drag (interfacial force) at an interface.

[0008] The present invention provides a method, apparatus, and program for fluid simulation that utilizes basic equations of flow with embedded solid boundaries to easily perform fluid simulation. [Means for solving the problem]

[0009] A first aspect of the present invention is Prepare the shape data of the solid part, setting an arbitrary computational grid for the shape data; defining a porosity corresponding to the shape of the solid portion as an index of ease of fluid flow for the computational grid; Perform a fluid simulation based on the following modified Navier-Stokes equations including porosity: TIFF2025117536000002.tif4788v: Flow velocity vector p: pressure ρ: Density ε: porosity f w :Interfacial force T: Viscous stress tensor matrix μ: 1st viscosity coefficient λ: Second viscosity coefficient I: identity matrix The present invention provides a method for fluid simulation, comprising:

[0010] This method incorporates porosity into the modified Navier-Stokes equations (basic flow equations with embedded solid boundaries), allowing fluid and solid regions to be represented by porosity. This eliminates the need to create complex computational grids that fit the interfaces, simplifying fluid simulations. While the computational grid itself is still created, it is no longer necessary to create a computational grid tailored to the shape of the solid. This allows for the adoption of a uniform computational grid, potentially making analysis time independent of the complexity of the solid's shape. Here, porosity is defined for the computational grid as a distribution of flowability, ranging from 0 to 1. For example, porosity is defined as 1 in the fluid region, 0 in the solid region, and 0.5 at the interface. This allows the fluid region to follow the conventional Navier-Stokes equations, while the boundary and solid regions follow the modified Navier-Stokes equations. Therefore, boundary and solid regions can be defined by appropriately defining porosity that corresponds to the shape of the solid. In the above equation, ∇p represents the pressure gradient vector, ∇v represents the gradient tensor matrix of the flow velocity vector, and ∇V t denotes the transpose matrix of ∇V, and ∇·V denotes the divergence of the velocity vector.

[0011] The interfacial force may be expressed in the no-slip condition as follows: TIFF2025117536000003.tif1357α:Constant Δ: boundary thickness μ: viscosity coefficient

[0012] This method allows for the specific formulation of modified Navier-Stokes equations that take into account the interfacial forces under no-slip conditions, whereas under slip conditions the interfacial forces are zero.

[0013] The shape data may be a group of coordinate position data of a solid boundary, The coordinate position data group may be converted into a signed distance function to the computational grid to define the porosity by the following formula: TIFF2025117536000004.tif1639ξ:Signed distance function Δ: boundary thickness

[0014] This method makes it easy to define the porosity corresponding to the shape of the solid part. The signed distance function itself is a well-known function. Specifically, the signed distance function is defined to take a positive value in the fluid region and a negative value in the solid region, indicating the distance from the interface. The porosity defined as above can be regarded as a level set function using a sigmoid function.

[0015] The shape data may be black and white data; The pixel values of the monochrome data may be normalized to a range of 0 to 1; The normalized pixel value may be defined as the porosity.

[0016] This method makes it easy to define porosity corresponding to the shape of a solid part for black-and-white data. Specifically, the porosity is normalized so that it takes a value of 1 in a white region (fluid region) with a pixel value of 255 and a value of 0 in a black region (solid region) with a pixel value of 0. Here, black-and-white data includes image data and CAD (Computer Aided Design) data. The black-and-white data may be two-dimensional data or three-dimensional data.

[0017] The shape data may be color data; The color data may be converted to grayscale data; Pixel values of the grayscale data may be normalized to a range of 0 to 1; The normalized pixel value may be defined as the porosity.

[0018] According to this method, by applying a grayscale process to the color data, it is possible to easily define the porosity corresponding to the shape of the solid part, just as with the black-and-white data described above. Here, the color data includes image data and CAD data. The color data may be two-dimensional data or three-dimensional data.

[0019] In the grayscale data, an area where the pixel value is greater than a predetermined value may be corrected to white.

[0020] With this method, grayscale data may contain noise-like portions with erroneous pixel values, such as portions that are actually fluid regions but have pixel values that are not fluid regions, but these portions can be corrected and treated as fluid regions, thereby improving calculation accuracy.

[0021] The shape data may be diffusion filtered.

[0022] According to this method, the contours of solid parts can be smoothed by diffusion filtering, for example, a Gaussian filter.

[0023] The shape data may be noise filtered.

[0024] This method allows for more accurate handling of shape data by removing dot noise through noise filtering. For example, this method is effective when a fluid region is erroneously displayed as dots within a solid region, or when a solid region is erroneously displayed as dots within a fluid region. For example, the noise filter is a median filter (median blur processing).

[0025] For time-dependent porosity, the momentum conservation formulation of the modified Navier-Stokes equations may be rearranged to: TIFF2025117536000005.tif16125The interfacial force may be expressed as the following equation under no-slip conditions for the time-dependent porosity: TIFF2025117536000006.tif1050α:Constant Δ: boundary thickness μ: viscosity coefficient Vr: Relative velocity between the flow velocity vector V and the velocity Vw of the solid part

[0026] This method allows numerical simulations to be performed even when the porosity is time-dependent, i.e., when the solid part is moving.

[0027] When calculating the pressure field, a quadtree method may be used in which the area of the calculation grid is divided into four and parallel calculations are performed.

[0028] This method can accelerate calculations, and in implementation, it has achieved a calculation speed 90 times faster than a general calculation method that does not use parallel calculations.

[0029] A second aspect of the present invention is a porosity definition unit that sets an arbitrary computational grid for shape data of a solid portion and defines a porosity corresponding to the shape of the solid portion for the computational grid as an index of ease of fluid flow; a computing unit that performs a fluid simulation based on the following Navier-Stokes equations including the porosity: TIFF2025117536000007.tif4788v: Flow velocity vector p: pressure ρ: Density ε: porosity f w :Interfacial force T: Viscous stress tensor matrix μ: 1st viscosity coefficient λ: Second viscosity coefficient I: identity matrix The present invention provides an apparatus for fluid simulation, comprising:

[0030] A third aspect of the present invention is Set an arbitrary calculation grid for the shape data of the solid part, defining a porosity corresponding to the shape of the solid portion as an index of ease of fluid flow for the computational grid; Perform a fluid simulation based on the following modified Navier-Stokes equations including porosity: TIFF2025117536000008.tif4788v: Flow velocity vector p: pressure ρ: Density ε: porosity f w :Interfacial force T: Viscous stress tensor matrix μ: 1st viscosity coefficient λ: Second viscosity coefficient I: identity matrix The present invention provides a program for fluid simulation that causes a computer to execute a fluid simulation method including the steps of: [Effects of the Invention]

[0031] According to the present invention, in a method, apparatus, and program for fluid simulation, fluid simulation can be easily performed by utilizing basic equations of flow in which solid boundaries are embedded. [Brief explanation of the drawings]

[0032] [Figure 1] 1 is a block diagram of a fluid simulation apparatus according to a first embodiment of the present invention. [Figure 2] 1 is a flowchart of a fluid simulation method according to the first embodiment. [Figure 3] 4 is a flowchart showing pre-processing in the first embodiment. [Figure 4] 4 is a flowchart showing the definition of porosity in the first embodiment. [Figure 5] Image of solid part shape data converted to binary data. [Figure 6] Contour plot showing the signed distance field of Figure 5. [Figure 7] A contour map of Figure 6 with a computational grid added and downscaled. [Figure 8] Contour plot showing the porosity distribution of Figure 7 . [Figure 9] A diagram of a Koch curve at complexity level 0 converted into a simple closed curve. [Figure 10] A diagram of a Koch curve at complexity level 1 converted into a simple closed curve. [Figure 11] A diagram of a Koch curve at complexity level 2 converted into a simple closed curve. [Figure 12] A diagram of a Koch curve at complexity level 3 converted into a simple closed curve. [Figure 13] FIG. 3 is a diagram showing an example of a fluid simulation result in the first embodiment. [Figure 14] 10 is a flowchart showing pre-processing in the second embodiment. [Figure 15] 10 is a flowchart showing the definition of porosity in the second embodiment. [Figure 16] Black and white image as shape data of solid parts. [Figure 17] A diagram obtained by applying various filter processes to Figure 16. [Figure 18] FIG. 10 is a diagram showing an example of a fluid simulation result in the second embodiment. [Figure 19] 10 is a flowchart showing pre-processing in the third embodiment. [Figure 20] 10 is a flowchart showing the definition of porosity in the third embodiment. [Figure 21] This is a diagram of a color image representing shape data of a solid part that has been converted to grayscale. [Figure 22] Figure 21 after being resized with margins added. [Figure 23] Figure 22 after white correction. [Figure 24] FIG. 10 is a diagram showing an example of a fluid simulation result in the third embodiment. [Figure 25] FIG. 10 is a schematic diagram showing a calculation domain in the fourth embodiment. [Figure 26] FIG. 10 is a contour diagram showing the porosity distribution in the fourth embodiment. [Figure 27] FIG. 13 is a contour diagram showing a fluid simulation result in the first example of the fourth embodiment. [Figure 28] FIG. 13 is a contour diagram showing a fluid simulation result in the first example of the fourth embodiment. [Figure 29] Conceptual diagram showing the two-domain division of the iterative solution method. [Figure 30] A conceptual diagram showing division into four regions using the quadtree method. DETAILED DESCRIPTION OF THE INVENTION

[0033] Hereinafter, an embodiment of the present invention will be described with reference to the accompanying drawings.

[0034] (First embodiment) FIG. 1 is a block diagram of a fluid simulation apparatus 1 according to a first embodiment of the present invention.

[0035] The fluid simulation device 1 performs fluid simulations using modified Navier-Stokes equations (basic flow equations with embedded solid boundaries) that include porosity defined in correspondence with the shape of the solid part as an index of fluid flow ease for the shape data of the solid part.

[0036] The shape data of a solid part broadly includes data relating to the shape of a solid part in two and three dimensions. For example, the shape data of a solid part may be image data or CAD data that displays the solid part. Alternatively, the shape data of a solid part may be data that combines multiple two-dimensional image data, such as medical computed tomography (CT), to display the solid part in three dimensions.

[0037] The fluid simulation apparatus 1 of this embodiment includes an input unit 10, an output unit 20, a storage unit 30, and a control unit 40. The input unit 10, the output unit 20, the storage unit 30, and the control unit 40 are connected to each other via an appropriate bus. Specifically, the fluid simulation apparatus 1 can be configured by an information processing device (computer) such as a desktop personal computer, a notebook computer, a workstation, or a tablet terminal.

[0038] The input unit 10 is a part that generates or receives input data for the fluid simulation apparatus 1. For example, the input unit 10 is configured by a keyboard, a mouse, a touch panel, or the like.

[0039] The output unit 20 is a part that displays the processing results, etc., of the control unit 40. For example, the output unit 20 is configured by a liquid crystal display, an organic EL (organic electro-luminescence) display, a plasma display, or the like.

[0040] The storage unit 30 is a recording medium that records various information. For example, the storage unit 30 is configured by a flash memory, a solid-state device (SSD), a hard disk, other storage devices, or an appropriate combination of these. The storage unit 30 stores the program executed by the control unit 40 and various information for executing the fluid simulation.

[0041] The control unit 40 is a controller that controls the fluid simulation apparatus 1. The control unit 40 includes a general-purpose processor such as a CPU (Central Processing Unit) or an MPU (Micro Processing Unit) that executes programs to achieve predetermined functions. The control unit 40 invokes and executes programs stored in the storage unit 30 to achieve various controls in the fluid simulation apparatus 1. Therefore, the control unit 40 achieves predetermined functions through cooperation between hardware and software. Note that the control unit 40 may be a hardware circuit designed specifically to achieve predetermined functions. That is, the control unit 40 can be realized not only by a general-purpose processor such as a CPU or an MPU, but also by various processors or circuits such as a GPU (Graphics Processing Unit), an FPGA (Field Programmable Gate Array), a DSP (Digital Signal Processor), or an ASIC (Application Specific Integrated Circuit).

[0042] The control unit 40 includes, as functional components, a pre-processing unit 41, a porosity definition unit 42, and a calculation unit 43. These are realized by the cooperation of the above hardware and software. These may also be interpreted as corresponding circuits.

[0043] The preprocessing unit 41 processes the shape data of the solid portion into a form suitable for analysis before the definition of porosity, which will be described later. Specifically, the preprocessing unit 41 may perform resizing processing to process the shape data of the solid portion into an arbitrary size by adding necessary margins to the shape data of the solid portion. Furthermore, if the shape data of the solid portion is color data, the preprocessing unit 41 may convert the color data into grayscale data. Furthermore, the preprocessing unit 41 may convert pixel values equal to or greater than a predetermined threshold value into the maximum pixel value (white). Furthermore, the preprocessing unit 41 may perform various filtering processes, such as noise filtering and diffusion filtering, on the shape data of the solid portion. Note that the preprocessing unit 41 is not an essential component and may be omitted as necessary.

[0044] The porosity definition unit 42 is a unit that sets a computational grid for the shape data of the solid part processed by the preprocessing unit 41 and defines the porosity corresponding to the computational grid. In this embodiment, the shape data of the solid part is a set of coordinate position data of the solid boundary. The porosity definition unit 42 derives a signed distance field from the set of coordinate position data and defines the porosity based on the signed distance field. The definition of the porosity will be described in detail later.

[0045] The calculation unit 43 uses the porosity defined by the porosity definition unit 42 to perform a fluid simulation according to the following modified Navier-Stokes equations (basic equations of flow with embedded solid boundaries):

[0046]

number

[0047] The modified Navier-Stokes equation (1) is the conventional Navier-Stokes equation with the porosity ε introduced. Here, the interfacial force f w is expressed as a vector related to the viscous drag at the interface. Also, the viscous stress T is expressed as a tensor related to the stress caused by viscosity. In the above equation, ∇p represents the pressure gradient vector, ∇v represents the gradient tensor matrix of the flow velocity vector, and ∇V tdenotes the transpose matrix of ∇V, and ∇·V denotes the divergence of the velocity vector. Here, we consider an incompressible flow and the density is constant. However, equation (1) and this procedure can also be applied to compressible flows where the density ρ is distributed.

[0048] Under slip conditions, the interfacial force f w will be 0.

[0049] Under no-slip conditions, the interfacial force f w is formulated as follows:

[0050] If we assume that the porosity ε does not change over time (the solid region is stationary), the above formula (1) can be mathematically transformed into the following formula (2).

[0051]

number

[0052] In the vicinity of the solid part and inside (1>ε→0), if the mean interface curvature κ=∇·n is finite, then |∇ε| / ε will be a finite value (constant 2 if the plane κ=0), and the term on the right side will become a penalty function, giving v·n→0 (i.e., wall condition) in the time evolution calculation. At this time, the law of conservation of mass is approximately satisfied. On the other hand, the interfacial force f in the above equation (2) w simulates the Darcy resistance in porous media, but when the flow is parallel to the wall (v·n→0), the shear stress with the solid wall is proportional to the velocity near the wall. The interfacial force acting on the assumed boundary thickness △ is expressed as follows:

[0053]

number

[0054] Under no-slip conditions, the wall shear stress τ w When this occurs, it is further defined as in the following equation (4).

[0055]

number

[0056] This can be evaluated by dimensional analysis as shown in equation (5) below.

[0057]

number

[0058] Here, following the conservative reinitialization of the level set method, the porosity ε can be given using a sigmoid function as shown in the following equation (6): Note that the signed distance function ξ is a parameter representing the signed distance field mentioned above.

[0059]

number

[0060] In this case, the interfacial force f is expressed as a polynomial approximation as shown in the following equation (7): w is required.

[0061]

number

[0062] The constant α can be given as α=32 by one-dimensional numerical analysis.

[0063] The calculation unit 43 calculates the interfacial force f formulated as above under the no-slip condition. w In the case of a slip condition, the calculation unit 43 calculates the interfacial force f wis set to 0 and calculations are performed based on the above formula (1). In addition, in the above formula (7), when the porosity ε changes over time (when the solid region moves), accurate calculations are possible by replacing the flow velocity vector V with the relative velocity Vr between the velocity of the fluid region and the velocity of the solid region (=flow velocity vector V - velocity Vw of the solid region).

[0064] FIG. 2 shows a flowchart of the fluid simulation method according to this embodiment.

[0065] When the fluid simulation method of this embodiment is started, shape data of a solid part is prepared (step S1), and the shape data is preprocessed into a form suitable for analysis (step S2). This preprocessing may be performed by a preprocessing unit 41. Next, an arbitrary computational grid is set for the shape data of the solid part, and porosity corresponding to the shape of the solid part is defined for the computational grid (step S3). This porosity definition may be performed by a porosity definition unit 42. Then, a calculation is performed based on the modified Navier-Stokes equation of the above formula (1) including the defined porosity (step S4). This calculation may be performed by a calculation unit 43. Then, the results of the fluid simulation are output (step S5). The results of the fluid simulation may be output to the output unit 20.

[0066] 3 and 4, the preprocessing (step S2) and the definition of porosity (step S3) will be described using a specific example. Here, the case where two-dimensional color image data is prepared as shape data of the solid part will be described as an example.

[0067] Fig. 3 is a flowchart showing an example of pre-processing in the fluid simulation method of this embodiment. Fig. 4 is a flowchart showing an example of porosity definition in the fluid simulation method of this embodiment.

[0068] 3, in the preprocessing (step S2), the color image data is binarized (step S2-11), and the binarized image data is resized to a size and shape suitable for analysis (for example, a square of a desired size) (step S2-12). Any method, such as adding margins, can be used for resizing.

[0069] Figure 5 shows the binary image obtained by the above preprocessing. The binary image is a plan view with a circular solid part placed in the center, and pixel values are shown in shades of gray. The black circular area with a pixel value of 0 represents the solid area (solid part), and the white area with a pixel value of 255 represents the fluid area. Note that if binary image data such as that shown in Figure 5 is given from the beginning, the preprocessing (step S2) can be omitted.

[0070] 4, the porosity definition unit 42 derives a signed distance field from the shape data of the solid part (step S3-11), and calculates a level set function based on the signed distance field (step S3-12). The calculated level set function is defined as the porosity for an appropriate computational grid (step S3-13). Note that the method for defining the porosity is not limited to this, and multiple methods can be adopted, as will be described in the second and third embodiments below.

[0071] Figure 6 is a contour plot showing the signed distance field of Figure 5 with a high pixel count. Figure 7 is a contour plot showing the signed distance field obtained by downscaling the number of pixels in Figure 6 and adding a computational grid corresponding to the downscaled number of pixels. Figure 8 is a contour plot showing the porosity distribution based on the signed distance field of Figure 7. Note that the dashed circles in Figures 6 to 8 indicate the part corresponding to the interface in Figure 5 (signed distance field 0 in Figures 6 and 7, porosity 0.5 in Figure 8). In addition, in Figure 8, the solid circle indicates the part with a porosity of 0.8.

[0072] Specifically, the porosity definition unit 42 derives a signed distance field for the binary image data shown in FIG. 5, with the interface between the solid and fluid regions set to 0 (see FIG. 6). A common method for this is the fast marching method. Here, the signed distance field (signed distance function ξ) derived using the fast marching method is applied to the image data by performing the operation of Equation (6) above to derive the porosity distribution stored in a two-dimensional array. This two-dimensional array may be treated as a regular, orthogonal computational grid, and the porosity corresponding to the shape of the solid part may be defined for that computational grid. However, in practice, an appropriate grid size is desirable when applying it to a fluid simulation. Therefore, downscaling is preferable when applying it to a fluid simulation at a resolution smaller than the image size. Here, to improve shape reproducibility, a signed distance field with a distance unit of 1 pixel is calculated, as shown in FIG. 7, and downscaling is performed by applying a scale factor when resizing to an arbitrary grid size. Then, the level set function was calculated by manipulating the above equation (6) and defined as the porosity (see Fig. 8). As a result, as shown by the shading in Fig. 8, the porosity is defined in the range from 0 to 1 for the computational grid.

[0073] The following describes the results of verification of the analysis time of the fluid simulation of this embodiment. In the verification, the behavior of fluid around a figure in which the Koch curve was converted into a simple closed curve was analyzed for each complexity level.

[0074] Figure 9 shows a diagram of a Koch curve at complexity level 0 that has been converted into a simple closed curve. Figure 10 shows a diagram of a Koch curve at complexity level 1 that has been converted into a simple closed curve. Figure 11 shows a diagram of a Koch curve at complexity level 2 that has been converted into a simple closed curve. Figure 12 shows a diagram of a Koch curve at complexity level 3 that has been converted into a simple closed curve.

[0075] The Koch curve is a type of fractal shape that has the property of gradually evolving from an initial shape into line segments. As shown in Figures 9 to 12, we considered an equilateral triangle as the initial shape (complexity level 0), and by performing the same operations as the Koch curve on each side, we formed a simple closed curve that gradually evolves according to the complexity level. In Figures 9 to 12, the inside of the closed curve is painted black for the purpose of fluid simulation, and is treated as a solid region.

[0076] As a preparation for the calculation, Koch curves with complexity levels 0 to 6 were converted into simple closed curves and generated as black and white image data, and the respective level set functions were calculated as described above. Note that illustrations of complexity levels 4 to 6 are omitted.

[0077] For the verification, a level set function (i.e., porosity) was defined for a 1215 x 842 pixel (computational grid point) level set with a boundary thickness Δ / dx = 2.5 for a 3645 x 2526 pixel image, equivalent in size to a 4K image. The porosity distribution was generated within several tens of seconds. The fluid simulation involved a uniform velocity inflow (right) with an angle of attack of 0°, a uniform pressure outflow (left), and periodic conditions (top and bottom) around the periphery of the computational domain, with a no-slip condition imposed on the boundary region. Discretization was performed using a second-order central difference method in space and a first-order explicit Euler method in time, with the collocated grid MAC (Marker and Cell) method used for time evolution calculations. For a flow field equivalent to a Reynolds number Re of 20 (uniform flow velocity, approximate radius), a converged solution was obtained within 2000 steps under the conditions of computational stability (CFL < 1, grid Pe < 2).

[0078] The analysis time under the above calculation conditions and the complexity of each level of the Koch curve are summarized in Table 1 below. As a result, the analysis time was approximately 420 seconds, regardless of the complexity level. Unlike this embodiment, when the computational grid is formed to fit the interface of the solid part, the computational grid also becomes more complex as the complexity level increases, resulting in a longer analysis time. However, in the fluid simulation of this embodiment, it was confirmed that the analysis time was roughly constant, regardless of the complexity level.

[0079] [Table 1]

[0080] As an example of the results of the fluid simulation, Fig. 13 shows the fluid behavior (visualized image of the streamlines) around a figure in which a Koch curve of complexity level 1 is converted into a simple closed curve.

[0081] 13 closely matches the actual fluid behavior, confirming that the fluid simulation of this embodiment has high analytical accuracy. Therefore, the fluid simulation of this embodiment does not require the creation of a complex computational grid that matches the interface, and achieves an analysis time that is not dependent on the shape of the solid part while maintaining high analytical accuracy.

[0082] The fluid simulation of this embodiment provides the following advantages.

[0083] Since the modified Navier-Stokes equation shown in equation (1) above includes porosity, the fluid and solid regions can be represented by the porosity. This eliminates the need to create a complex computational grid that matches the interface, allowing for easy fluid simulation. Although a computational grid is still created in the fluid simulation of this embodiment, it is no longer necessary to create a computational grid that matches the shape of the solid part. This allows for the use of a uniform computational grid, making the analysis time independent of the complexity of the shape of the solid part.

[0084] Moreover, under the no-slip condition, the interfacial force is defined as in equation (7), so that the modified Navier-Stokes equations that take the interfacial force into account can be specifically formulated.

[0085] Furthermore, by using a signed distance field, it is possible to easily define the porosity corresponding to the shape of the solid part.

[0086] (Second embodiment) In the second embodiment, the definitions of pretreatment and porosity differ from those in the first embodiment. Other than these, the second embodiment is substantially the same as the first embodiment. Therefore, the explanation of the parts shown in the first embodiment may be omitted. Furthermore, in this embodiment, black-and-white data is used as the shape information data of the solid part.

[0087] Fig. 14 is a flowchart showing an example of pre-processing of the fluid simulation method of this embodiment. Fig. 15 is a flowchart showing an example of the definition of porosity in the fluid simulation method of this embodiment. Note that the overall configuration of the fluid simulation device 1 of this embodiment is the same as in Fig. 1, and the overall flowchart of the fluid simulation method of this embodiment is the same as in Fig. 2, so they are not shown here.

[0088] 14, in pre-processing (step S2), the prepared black-and-white image data is resized to a size and shape suitable for analysis (for example, a square of a desired size) (step S2-21), and various filtering processes are performed (step S2-22). In this embodiment, as various filtering processes, a diffusion filtering process (for example, a Gaussian filter) is performed on the black-and-white image data to smooth the transition of pixel values at the outline of solid parts.

[0089] Figures 16 and 17 show the black-and-white image data obtained by the above preprocessing. The black-and-white image data is a front view with a human-shaped solid part placed in the center. Figure 16 shows the binarized image resized to a square, and Figure 17 shows the image of Figure 16 after applying diffusion filter processing.

[0090] In the case of black-and-white image data, if only the outline of a solid portion is given, the interior of the outline may be filled in during preprocessing. Furthermore, if an illustration image is given, the illustration image may be subjected to noise filtering during preprocessing to remove dot noise. In this way, the preprocessing may be changed as appropriate depending on the given black-and-white image data.

[0091] Referring to FIG. 15, in the definition of porosity (step S3), pixel values of the black and white image data that has been subjected to the above-mentioned filtering process are normalized within a range of 0 to 1 (step S3-21), and the normalized pixel values are defined as porosity (step S3-22). Specifically, porosity is normalized so that it takes a value of 1 in a white region (fluid region) with a pixel value of 255, and a value of 0 in a black region (solid region) with a pixel value of 0. In this embodiment, the image pixels of the black and white image data of the original image are used as a regularly spaced computational grid, and porosity is defined for this computational grid. Therefore, the original image is reproduced faithfully.

[0092] The results of verification of the analysis accuracy when the fluid simulation of this embodiment is implemented will be described.

[0093] For the verification, a 901 x 901 pixel (computational grid point) original image was subjected to a diffusion filter process equivalent to 2 pixels. For the fluid simulation, a uniform velocity inflow (right) with an angle of attack of 20°, a uniform pressure outflow (left), and periodic conditions (top and bottom) were imposed on the periphery of the computational domain, with a no-slip condition imposed on the boundary domain. For discretization, a second-order accuracy central difference method was applied to space, and a first-order accuracy explicit Euler method was applied to time, with time evolution calculations performed using the collocated grid MAC method. For a flow field equivalent to a Reynolds number Re of 80 (uniform flow velocity, solid part height), a converged solution was obtained in 2000 steps under conditions that satisfied calculation stability (CFL<1, grid Pe<2).

[0094] As an example of the results of the fluid simulation, a visualized image of the flow lines from the above calculation results is shown in FIG.

[0095] The streamlines shown in FIG. 18 closely match the actual fluid behavior, confirming that the fluid simulation of this embodiment has high analytical accuracy.

[0096] The fluid simulation according to this embodiment provides the following advantageous effects.

[0097] By using pixel values for black and white data as shape data of a solid part, it is possible to easily define the porosity corresponding to the shape of the solid part.

[0098] In addition, the contours of solid parts can be handled smoothly by using a diffusion filter.

[0099] (Third embodiment) In the third embodiment, the definitions of pretreatment and porosity differ from those in the first and second embodiments. Other than these, the third embodiment is substantially the same as the first and second embodiments. Therefore, the explanation of the parts shown in the first and second embodiments may be omitted. Furthermore, in this embodiment, color data is used as shape information data for solid parts.

[0100] Fig. 19 is a flowchart showing an example of pre-processing of the fluid simulation method of this embodiment. Fig. 20 is a flowchart showing an example of the definition of porosity in the fluid simulation method of this embodiment. Note that the configuration of the fluid simulation device 1 of this embodiment is the same as in Fig. 1, and the flowchart of the fluid simulation method of this embodiment is the same as in Fig. 2, so they are not shown here.

[0101] Referring to FIG. 19, in preprocessing (step S2), the prepared color image data is converted into a grayscale image (step S2-31) and resized to a size and shape (e.g., a square of the desired size) suitable for analysis (step S2-32). Next, regions with a predetermined pixel value or greater are corrected to white (step S2-33). Specifically, white correction is performed to set the pixel values of fluid regions (fluid regions excluding solid portions and the vicinity of the solid portions) with a predetermined pixel value or greater to the maximum value of 255. Then, various filter processes are performed on the image data that has undergone white correction (step S2-34). In this embodiment, noise filtering and diffusion filtering are performed as various filter processes. For example, the noise filtering is a median filter (Median blur processing). The diffusion filtering is the same as in the second embodiment.

[0102] Figures 21 to 23 show image data obtained by the above preprocessing. The image data is a plan view with a scissors-shaped solid part placed in the center. Figure 21 shows grayscaled image data, Figure 22 shows image data obtained by adding white space to Figure 21 and resizing it, and Figure 23 shows Figure 22 after white correction processing. As a result, the area other than the scissors has been corrected to white.

[0103] 20, the definition of porosity (step S3) involves normalizing the pixel values of the filtered image data within the range of 0 to 1 (step S3-31), and defining the normalized pixel values as porosity (step S3-32). In this embodiment, an arbitrary computational grid is set, and the porosity is defined for the computational grid.

[0104] The results of verification of the analysis accuracy when the fluid simulation of this embodiment is implemented will be described.

[0105] For the verification, a 1491 x 1491 pixel color image was subjected to the above preprocessing to generate a 400 x 400 grayscale image, followed by diffusion and noise filtering, and pixel values were normalized to define porosity. For the fluid simulation, a 0° angle of attack uniform velocity inflow (right), uniform pressure outflow (left), and periodic conditions (top and bottom) were imposed on the periphery of the computational domain, with a no-slip condition imposed on the boundary domain. Discretization was performed using second-order accuracy central differences in space and first-order accuracy explicit Euler methods in time, with the collocated grid MAC method used for time evolution calculations. For a flow field equivalent to a Reynolds number Re of 80 (uniform flow velocity, solid length), a converged solution was obtained within 1000 steps under conditions that satisfied computational stability (CFL<1, grid Pe<2).

[0106] As an example of the results of the fluid simulation, a visualized image of the flow lines from the above calculation results is shown in FIG.

[0107] The streamlines shown in FIG. 24 closely match the actual fluid behavior, confirming that the fluid simulation of this embodiment has high analytical accuracy.

[0108] According to this embodiment, the following advantageous effects are achieved.

[0109] By performing grayscale processing on color data, it is possible to easily define the porosity corresponding to the shape of the solid part, in the same way as with the black and white data in the second embodiment.

[0110] Furthermore, grayscale images may contain noise-like areas with erroneous pixel values, such as areas that are actually fluid but have pixel values that are not fluid, but these areas can be corrected and treated as fluid areas, thereby improving calculation accuracy.

[0111] Furthermore, noise filtering can remove dot noise, allowing for more accurate handling of shape data. This is effective when, for example, a fluid region is erroneously displayed as dots within a solid region, or when a solid region is erroneously displayed as dots within a fluid region.

[0112] (Fourth embodiment) The fourth embodiment differs from the first to third embodiments in that it is possible to calculate the time-varying porosity ε. Other than this, the fourth embodiment is substantially the same as the first to third embodiments. Therefore, the explanation of the parts shown in the first to third embodiments may be omitted.

[0113] In this embodiment, a method will be described that enables calculation when the porosity ε changes with time, that is, when the solid region moves.

[0114] The equation (1) relating to the law of conservation of momentum is transformed as follows:

[0115]

number

[0116] Here, for simplicity, we will first consider one dimension, and dimensional expansion will be discussed later. If we use the FTCS (Forward-Time Centered-Space) scheme, which applies second-order accuracy central differences to space and first-order accuracy explicit Euler method to time, we can discretize the above equation (8) as follows:

[0117]

number

[0118] (A) to (F) in the above equation (9) represent the advection term, diffusion term, viscosity term (diffusion term), additional term due to porosity distribution, interfacial force term, and pressure term, respectively, and all are discretized by second-order central difference. (A) to (E) represent the state at time k, but only (F) represents the state at time k+1. Therefore, the velocity field calculated from (A) to (E) is called the virtual velocity field v + Toki, v + The pressure field is calculated using

[0119] Assuming time-varying porosity ε and an incompressible fluid, the pressure Poisson equation can be expressed as follows:

[0120]

number

[0121] The above formula (10) is transformed as follows using second-order accuracy central difference:

[0122]

number

[0123] The simultaneous linear equations derived in this way are calculated using an iterative solution method, etc., to obtain a tentative velocity field v + Pressure field p k+1 is calculated. The resulting pressure field p k+1By substituting into (F) in equation (9), the velocity field corrected to take pressure into account is obtained as follows:

[0124]

number

[0125] Here, the SOR (Successive Over-Relaxation) method was used as the iterative solution for pressure. The SOR method can change the convergence by using the relaxation coefficient w, and it is known that setting an appropriate value improves the convergence. The equation used in the SOR method is expressed as follows using the above equation (11). Note that when implementing, for example, the relaxation coefficient w may be set to 1.7 and the number of iterations may be set to 100.

[0126]

number

[0127] In this way, a formulation was created that allows calculations to be made when the porosity ε changes over time.

[0128] Below, we explain the verification results of fluid simulation for two examples in which the porosity ε changes over time.

[0129] As a first example, we show the case where the solid region is rotating. Figure 25 shows a schematic diagram of the calculation domain. Figure 26 shows the distribution of porosity ε.

[0130] In the first example, we applied the no-slip condition to a two-dimensional circular cylinder with a rotating velocity and evaluated the accuracy of the numerical simulation. The diameter of the cylinder, D [m], the angular velocity of rotation, ω [rad / s], and the inflow velocity, U in [m / s], and the dimensionless rotational speed α (= Dω / 2U in ) and Reynolds number Re(=DU in / ν) was defined, where ν is the dynamic viscosity coefficient. The X and Y velocity components of the solid region were defined as (-ωr·sinωt, ωr·cosωt), where r and t represent the distance [m] and time [s] from the center O of the cylinder, respectively.

[0131] The computational domain was set to -25D ≤ X, Y ≤ 25D, and a cylinder with diameter D = 20Δx [m] and boundary thickness Δ / dx = 1.5 was placed at the center (X, Y) = (0, 0). The boundary conditions were a velocity Dirichlet condition (v = (1, 0) [m / s]) and a pressure Neumann condition (∂p / ∂n = 0) on the inflow boundary (left boundary), and a velocity Neumann condition (∂v / ∂n = 0) and a pressure Dirichlet condition (p = 0 [Pa]) on the outflow boundary (right boundary). Periodic boundary conditions were used on the upper and lower boundaries. In this example, a fluid simulation was performed for a uniform flow with a Reynolds number of Re = 20, with five non-dimensional rotational velocities of α = 0, 0.25, 0.5, 1, and 2. The global step count was 50,000, the SOR iteration count was 100, and the computation time was 10 seconds.

[0132] Figure 27 shows a visualized image of the calculation results (step = 50000). The left image shows the flow velocity, and the right image shows the pressure.

[0133] As a result of the fluid simulation, it was observed that as the non-dimensional rotational speed α increases, the velocity difference between the top and bottom of the cylinder increases, resulting in a pressure difference, the Magnus effect. This is in good agreement with actual behavior, confirming the high accuracy of the fluid simulation.

[0134] As a second example, we consider the case where a solid region undergoes periodic motion.

[0135] In the second example, the no-slip condition was applied to a two-dimensional circular cylinder undergoing periodic motion, and the accuracy of the numerical simulation was evaluated. The central position of the cylinder was defined as (-A sinωt,0) and the velocity as (-ωA cosωt,0), where A is the amplitude [m]. Other quantities were the same as in the first example. The maximum velocity was defined as Umax [m / s], and the Couligan-Carpenter number KC (= 2πUmax / Dω) and the Reynolds number Re (= DUmax / ν) were defined.

[0136] The computational domain is -20D≦X,Y≦20D, and a cylinder with diameter D=100Δx[m] and boundary thickness Δ / dx=1.5 is placed at the center (X,Y)=(0,0) at t=0[s]. The boundary conditions are a Neumann condition (∂v / ∂n=0) for velocity in all directions and a Dirichlet condition (p=0[Pa]) for pressure. In this example, a fluid simulation was performed with a Reynolds number Re=100 and a Couligan-Carpenter number KC=5, with a global step of 12000[steps], SOR iterations of 100[steps], and a computation time of 2.0[s].

[0137] Figure 28 shows visualized images of the calculation results (step = 9000 (180°), 9250 (210°), 10250 (330°)). The left image shows the flow velocity, and the right image shows the pressure.

[0138] As a result of the fluid simulation, it was observed that low-pressure regions were generated above and below the cylinder, and that these regions moved while generating vortices. This was in good agreement with actual behavior, confirming the high accuracy of the fluid simulation.

[0139] Furthermore, the calculation can be speeded up by parallelization as follows.

[0140] The SOR method allows for parallel calculations by dividing the computational grid into two regions, odd points and even points. Here, when expanding the dimensions, the pressure solver was parallelized to shorten the execution time. Below is the equation expanded to two dimensions, equation (11).

[0141]

number

[0142] According to equation (14), in the iterative calculation of pressure, as shown in Figure 29, when updating the value of (i,j) (white), only the four adjacent stencils (i-1,j), (i,j-1), (i,j+1), and (i+1,j) (black) are used. In other words, the update of the same color is not affected. Therefore, parallelization is possible by dividing the white and black regions as shown in Figure 29 and solving them. In addition, in the case of three dimensions, the number of stencils used in the calculation increases to six, but this can be easily implemented by dividing the grid into two regions using the same procedure as for two dimensions.

[0143] We will also explain how to extend the area division using two colors to four colors using the quadtree method.

[0144] Referring to Figure 30, when dividing into two regions, the two-dimensional space must be divided into an odd number of rows and columns, whereas when dividing into four regions, it must be divided into an even number of rows and columns. Therefore, as shown in Figure 30, the array is shifted based on the even / odd dimension. By storing values in the array for each region based on the prepared array index, a continuous array with a shorter vector length than when dividing into two regions is obtained. These arrays are executed for the entire region in the order [f00,f01] → [f10,f11] (however, the values within [] can be parallelized), and by repeating this process, a solution is obtained for the SOR iteration. In this algorithm, the frequency of data updates per point is the same as in the conventional coordinate system order or the two-color division method described above, so the convergence of the SOR iteration remains essentially unchanged.

[0145] The advantage of this method is that there are no self-references within each color (and within []), allowing for full parallelization. Up to half the number of calculation points can be parallelized, as with the two-color division. By combining parallelization within each color (1 / 4 of the calculation points) with the [] (2-parallelization), multiple parallelization methods (e.g., multiple GPU calculations, vector calculations, and core parallelization) can be implemented in the program. In this case, by defining an independent array for each color (f00,...), the number of arrays becomes 1 / 4 of the calculation points, and access becomes continuous, which is expected to improve memory speed for GPU and vector calculations. Furthermore, by expanding the number of colors (e.g., 4x4, 2x2x2), it is possible to support multi-point stencils and 3D grids in high-order spatial differentiation. It is also possible to implement an increase in the number of core parallelizations.

[0146] In this way, we have established a method for parallelizing computational algorithms to increase their speed.

[0147] The results of verifying the execution speed of the above method are explained below.

[0148] In a square region (X, Y) with a side length of L, Poisson equation ∇ 2 The following conditions are set for φ=c: the left end (X=0) is the Neumann condition (∂φ / ∂x=0), the right end (X=L) is the φ=0 (Dirichlet condition), the upper and lower boundaries Y=0,L (z=0,L) are periodic conditions, and the generator term on the right side is positive for the central point source and 0 elsewhere.

[0149] The variable φ is arranged in an equally spaced orthogonal grid (number of grids m x m) in each direction along the Cartesian coordinate system, and a numerical solution is obtained by solving the discretized algebraic equation at the grid points. In the program, the spacing between neighboring grid points dx = dy = L / m is normalized to 1, the central four points are generated with a value of 1, and the others are normalized to 0.

[0150] The execution speed was compared using a single core and a GPU. The specifications of the computer used are shown in Table 2 below.

[0151] [Table 2]

[0152] Table 3 below shows the results of a comparison of the calculation times for non-parallel (normal), parallel (two-color), and parallel (quadtree) methods at 100 iterations.

[0153] [Table 3]

[0154] Compared to the general coordinate-ordered array (normal), the improved method (quadtree method) was accelerated even in single-core calculations, confirming that the parallelized program improves not only the calculation speed but also the data transfer speed. When the calculation size is large (512x512 or more), GPU acceleration is further activated, resulting in an acceleration of up to approximately 90 times compared to the conventional scalar calculation method (more than twice as fast as the two-color parallel method).

[0155] While specific embodiments of the present invention and their modifications have been described above, the present invention is not limited to the above-described embodiments and can be implemented with various modifications within the scope of the present invention. For example, an appropriate combination of the contents of the individual embodiments may be considered as one embodiment of the present invention.

[0156] Furthermore, each of the above embodiments may be realized by a program that causes a computer to execute each step constituting the above method. By executing these programs, it is possible to obtain the effects of the above method. In other words, it can be said that the above method is being used. [Explanation of symbols]

[0157] 1. Fluid simulation device 10 Input section 20 Output section 30 Storage section 40 Control Unit 41 Pretreatment section 42 Porosity definition part 43 Arithmetic section

Claims

1. Prepare the shape data of the solid part, setting an arbitrary computational grid for the shape data; defining a porosity corresponding to the shape of the solid portion as an index of ease of fluid flow for the computational grid; Perform a fluid simulation based on the following modified Navier-Stokes equations, including the porosity: v: flow velocity vector p: pressure ρ: density ε: porosity f ｗ : Interface force T: Viscous stress tensor matrix μ: first viscosity coefficient λ: second viscosity coefficient I: identity matrix 12. A method for fluid simulation, comprising:

2. The interfacial force is expressed as follows in no-slip conditions for the porosity, which is time independent: α: constant Δ: boundary thickness μ: viscosity coefficient The method for fluid simulation according to claim 1 .

3. the shape data is a group of coordinate position data of a solid boundary, The coordinate position data group is converted into a signed distance function to the computational grid, and the porosity is defined by the following formula: ξ: Signed distance function Δ: boundary thickness A method for fluid simulation according to claim 1 or 2.

4. the shape data is black and white data, normalizing pixel values of the monochrome data to a range of 0 to 1; 3. The method for fluid simulation according to claim 1, wherein the normalized pixel value is defined as the porosity.

5. the shape data is color data, converting the color data into grayscale data; normalizing pixel values of the grayscale data to a range of 0 to 1; 3. The method for fluid simulation according to claim 1, wherein the normalized pixel value is defined as the porosity.

6. The method for fluid simulation according to claim 5 , wherein areas in the grayscale data where the pixel values are greater than a predetermined value are corrected to white.

7. 3. A method for fluid simulation according to claim 1 or 2, wherein the shape data is diffusion filtered.

8. 3. A method for fluid simulation according to claim 1 or 2, wherein the shape data is noise filtered.

9. For the time-dependent porosity, the momentum conservation expression of the modified Navier-Stokes equations is rearranged as follows:

2. The method for fluid simulation of claim 1, wherein the interfacial force is expressed in the following equation at no-slip conditions for the time-dependent porosity: α: constant Δ: boundary thickness μ: viscosity coefficient Vr: relative velocity between the flow velocity vector V and the velocity Vw of the solid part

10. 10. The method for fluid simulation according to claim 1 or 9, wherein a quadtree method is used to calculate the pressure field by dividing the computational grid area into four and performing parallel calculations.

11. a porosity definition unit that sets an arbitrary computational grid for shape data of a solid portion and defines a porosity corresponding to the shape of the solid portion for the computational grid as an index of ease of fluid flow; a calculation unit that performs a fluid simulation based on the following Navier-Stokes equations including the porosity: v: flow velocity vector p: pressure ρ: density ε: porosity f ｗ : Interface force T: Viscous stress tensor matrix μ: first viscosity coefficient λ: second viscosity coefficient I: identity matrix 1. An apparatus for fluid simulation comprising:

12. Set an arbitrary calculation grid for the shape data of the solid part, defining a porosity corresponding to the shape of the solid portion as an index of ease of fluid flow for the computational grid; Perform a fluid simulation based on the following modified Navier-Stokes equations, including the porosity: v: flow velocity vector p: pressure ρ: density ε: porosity f ｗ : Interface force T: Viscous stress tensor matrix μ: first viscosity coefficient λ: second viscosity coefficient I: identity matrix A program for fluid simulation that causes a computer to execute a fluid simulation method including the steps of:

Images