WO2014069421A1 - シミュレーション装置、シミュレーション方法およびプログラム - Google Patents

シミュレーション装置、シミュレーション方法およびプログラム Download PDF

Info

Publication number
WO2014069421A1
WO2014069421A1 PCT/JP2013/079175 JP2013079175W WO2014069421A1 WO 2014069421 A1 WO2014069421 A1 WO 2014069421A1 JP 2013079175 W JP2013079175 W JP 2013079175W WO 2014069421 A1 WO2014069421 A1 WO 2014069421A1
Authority
WO
WIPO (PCT)
Prior art keywords
equation
term
cell
velocity
information
Prior art date
Application number
PCT/JP2013/079175
Other languages
English (en)
French (fr)
Inventor
一貴 柳原
恒洋 齊藤
Original Assignee
旭硝子株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 旭硝子株式会社 filed Critical 旭硝子株式会社
Priority to CN201380056818.9A priority Critical patent/CN104769593B/zh
Priority to EP13852061.4A priority patent/EP2916247A4/en
Priority to JP2014544502A priority patent/JP6164222B2/ja
Publication of WO2014069421A1 publication Critical patent/WO2014069421A1/ja
Priority to US14/699,858 priority patent/US20150234784A1/en

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present invention relates to a simulation apparatus, a simulation method, and a program.
  • Non-Patent Document 1 the advection term is called a convection term.
  • Non-Patent Document 2 and 3 the advection term is called a convection term.
  • the inspection volume V (t) is taken in the continuum.
  • the inspection volume V (t) also moves and deforms as the continuum moves and deforms.
  • the inspection volume V (t) at the position of the origin O ′ (t) moves to the position of O ′ (t + ⁇ t)
  • the inspection volume V The state of moving and deforming to the state of (t + ⁇ t) and surface S (t + ⁇ t) is shown. This movement and deformation is caused by the movement of the surface S (t) surrounding V (t) at the velocity of the velocity vector u s (t).
  • a partial differential equation representing a continuum is obtained by the limit operation of V ⁇ 0 and ⁇ t ⁇ 0.
  • a Lagrangian coordinate system is a coordinate system having an origin that moves with a continuum, such as an origin O ′ (t) in FIG. 19.
  • the Euler coordinate system is a coordinate system set separately from the continuum as shown by the origin O in FIG. In the Euler coordinate system, the origin O does not move with the continuum.
  • u, v, and w are the x, y, and z components of the velocity vector u at the position of the coordinates (x, y, z) in the continuum, and the velocity vector u is expressed by Equation (1-4). It is expressed in
  • equations (1-2) and (1-3) become equations (1-6) and (1-7).
  • the substance differential D / Dt shown in the formula (1-6) is expressed as Euler as a term representing the movement and deformation of the continuum over time for all physical quantities expressed as a continuum as shown below. Appears in partial differential equations expressed in a coordinate system. As shown in the equations (1-3) and (1-7), the material differential includes the advection term, and therefore the partial differential equation in which the material differentiation appears includes the advection term.
  • Equation (1-8) Fluid dynamics (Navier-Stokes equations)
  • the Navier-Stokes equation which is the governing equation of fluid mechanics, is expressed as equation (1-8).
  • t time
  • ⁇ fluid density
  • ⁇ fluid viscosity coefficient
  • p pressure.
  • the substance differentiation appears on the left side. That is, it includes an advection term.
  • Equation (1-10) Energy conservation equations Energy conservation laws can be divided into thermal energy conservation and kinetic energy conservation. The conservation of kinetic energy is the same as in equation (1-9), so here the equation for thermal energy conservation is shown in equation (1-10).
  • U internal energy of continuum
  • ⁇ density of continuum
  • ⁇ heat source (source / sink of heat energy )
  • Dij (i, j 1, 2, 3): deformation rate tensor of the continuum.
  • the substance differentiation appears on the left side. That is, it includes an advection term.
  • Advection-diffusion equation The advection-diffusion phenomenon of a substance C in the continuum is expressed by the following equation (1-11).
  • C concentration of substance C
  • ⁇ c diffusion coefficient of substance C
  • qc source / sink of substance C
  • ⁇ density of continuum
  • the substance differentiation appears on the left side. That is, it includes an advection term.
  • the advection term is a term that appears not only in the Navier-Stokes equation but also in an equation including material differentiation (hereinafter referred to as a transport equation).
  • a transport equation an equation used in the case of the finite volume method described in Non-Patent Document 1 will be described as an example of a method for solving the advection term by conserving type.
  • the transport equation for the variable ⁇ is expressed as in the equation (2-1) or the equation (2-2).
  • the second term on the left side of equation (2-2) is the advection term.
  • the Navier-Stokes equation can be said to be a transport equation in which the variable ⁇ is the speed u or the momentum ⁇ u.
  • Equation (2-5) is a conservative advection term when the density is constant.
  • Expression (2-7) is added to Expression (2-7).
  • Expression (2-9) is obtained.
  • Expression (2-10) is obtained.
  • the second term on the left side of Equation (2-10) is a conservative advection term when the density is not constant (the density changes).
  • Equation (2-10) By replacing ⁇ ⁇ ⁇ ′ and ⁇ F ( ⁇ ) ⁇ F ′ ( ⁇ ′) in Equation (2-10), Equation (2-10) can be expressed in the same form as Equation (2-5). Therefore, the following description is based on the conservative transport equation of the formula (2-5).
  • the equation (2-15) is a discretized equation obtained by volume integration of the conservative transport equation for the variable ⁇ . That is, the physical quantity ⁇ is an expression representing a state of advection as the continuum moves and deforms. If ⁇ is the same value in all regions of the continuum at a certain point in time, even if the continuum is moving or deforming, the result of the simulation by the equation (2-15) is ⁇ in all regions of the continuum. Need to be the same value without changing (Condition A). For example, when the physical quantity ⁇ is the salinity concentration of water and the salinity concentration ⁇ is the same in all regions, the salinity concentration ⁇ does not change no matter how the water flows over time. For this condition A, consider the one-dimensional flow problem shown in FIG.
  • equation (2-22) is obtained from equation (2-22).
  • the equation (2-23) is phi a show to be a constant value without fluctuation even if progression of time. This satisfies the condition A described in the explanation of the expression (2-15).
  • the condition A is not satisfied even if the equation (2-15) is used.
  • the one-dimensional flow problem is considered as an example, but the same applies to the two-dimensional and three-dimensional transport equations. That is, the continuous equation (2-14) is not completely satisfied, and when it is the equation (2-25), even if ⁇ is constant in the entire calculation region, the equation (2-26) is obtained. phi a can will be increased or decreased with time progresses, condition A is not satisfied.
  • Non-Patent Documents 2 and 3 pose the problem that even if the integral equation analytically satisfies the conservation law, the conservation law is not satisfied when these are discretized. However, Non-Patent Documents 2 and 3 present solutions to the problem that the conservation law is not satisfied by modifying the discretization formula so as to satisfy the conservation law. That is, Non-Patent Documents 2 and 3 do not present any solution for the problem of error in numerical calculation.
  • the present invention has been made in view of such circumstances, and an object thereof is to provide a simulation apparatus, a simulation method, and a program capable of obtaining a good simulation result.
  • the present invention has been made to solve the above-described problems, and one aspect of the present invention is a simulation apparatus that calculates a solution of an equation including an advection term, and acquires a velocity field that acquires a velocity field. And a solution calculation unit that calculates a solution in the discretization equation by applying the velocity field acquired by the velocity field acquisition unit to a discretization equation obtained by adding at least two operations to the equation And the two operations include transforming the advection term into a conservative type and including a term for correcting an error in the conservative type of the advection term.
  • the term for correcting the error is obtained by performing at least the same integration as the integration for the equation on the continuous equation. It is a term.
  • the solution calculation unit calculates a velocity field at a next time step using at least a continuous equation, and acquires the velocity field.
  • the velocity field acquired by the unit is the velocity field calculated by the solution calculation unit.
  • the other aspect of this invention is the above-mentioned simulation apparatus, Comprising:
  • the said velocity field acquisition part acquires the velocity field memorize
  • a simulation method for calculating a solution of an equation including an advection term includes at least a first process for acquiring a velocity field and an equation including an advection term.
  • a second step of calculating a solution in the discretization equation by applying the velocity field obtained in the first step to the discretization equation obtained by adding two operations, and the two operations Is characterized in that the advection term is transformed into a conservative type and a term for correcting the conservative type error of the advection term is included.
  • a computer is added to the discretization equation obtained by adding at least two operations to an equation including a velocity field acquisition unit for acquiring a velocity field and an advection term.
  • a simulation apparatus for estimating a generation source of a diffusing material that estimates fluid generation source information based on information from a plurality of observation devices.
  • Observation information obtaining means for obtaining the position information and measured concentration information
  • virtual emission point setting means for setting a plurality of virtual emission points in a region where the concentration of the diffused substance is to be estimated, position information of the observer
  • the above-described simulation device that estimates the concentration of the diffusing substance at the virtual release point based on the measured concentration information and the information on the virtual release point.
  • An observation information obtaining step for obtaining the position information and the measured concentration information, a virtual release point setting step for setting a plurality of virtual release points in a region where the concentration of the diffused substance is to be estimated, position information of the observer, And a step of estimating the concentration of the diffusing substance at the virtual release point by the above-described simulation method based on the measured concentration information and the information on the virtual release point.
  • FIG. 1 is a schematic block diagram showing a configuration of a fluid simulation apparatus 10 according to a first embodiment of the present invention. It is a conceptual diagram explaining the cell in the same embodiment. It is a figure which shows the example of the cell information which the cell structure memory
  • Equation (3-1) the speed obtained by simulation and measurement includes an error, and therefore does not completely satisfy Equation (3-1), and the right side does not become zero. Therefore, here, a conservative transport equation is derived using an equation (3-2) in which the right side of the equation (3-1) is replaced from zero to D. That is, the conservative transport equation is derived on the assumption that the continuous equation is not completely satisfied and the right side of the equation (3-1) is not zero.
  • this equation (3-5) when the density ⁇ is constant, this equation (3-5) is used as a conservative transport equation. That is, in the present invention, the solution of Expression (3-5) is calculated.
  • This (3-5) is the same as the conventional conservative transport equation (2-5) except that it has the third term ( ⁇ ⁇ D) on the left side.
  • the third term on the left-hand side as shown in the equation (3-2), the right-hand side of the continuous equation is not zero, and a conservative transport equation is assumed that includes the error D. ing.
  • the pseudo source / sink that appears in the equation (2-24) can be corrected by ⁇ ⁇ D. That is, it is possible to prevent a pseudo source / sink from appearing.
  • Equation (3-5) As an example of simulation using Equation (3-5), a case where Equation (3-5) is discretized by the finite volume method and simulated will be described. In the following, the case of discretization by the finite volume method is described. However, any other discretization method can be used as long as it is a method of discretizing an equation having a conservative advection term such as a difference method. Also good.
  • equation (3-5) is integrated with the volume V of the control volume to obtain the equation (3-6).
  • equation (3-6) is the same as that in formula (2-11).
  • equation (3-2) is also integrated with respect to the volume V of the control volume, equation (3-7) is obtained.
  • Expression (3-9) is substituted into Expression (3-8), Expression (3-10) is obtained.
  • a good simulation result can be obtained by solving the discretized transport equation (3-10).
  • this equation (3-14) when the density ⁇ is not constant, this equation (3-14) is used as a conservative transport equation. That is, the solution of equation (3-14) is calculated.
  • This (3-14) is the same as the conventional conservative transport equation (2-10) except that it has the third term ( ⁇ ⁇ D) on the left side.
  • the third term on the left side as shown in the equation (3-11), the right side of the continuous equation does not become zero, and a conservative transport equation is assumed that includes the error D. ing.
  • the pseudo source / sink that appears in the equation (2-24) can be corrected by - ⁇ ⁇ D. That is, even when the density ⁇ is not constant, the pseudo source / sink can be prevented from appearing.
  • the equation (3-14) is discretized by the finite volume method and the simulation is performed in the same manner as when the density ⁇ is constant.
  • the expression (3-14) is obtained by integrating the expression (3-14) with the volume V of the control volume.
  • Each symbol in formula (3-15) is the same as that in formula (2-11).
  • the equation (3-16) is obtained.
  • equation (3-19) is obtained.
  • the expression (3-19) becomes the expression (3-21).
  • ⁇ ab is the value of ⁇ at the boundary surface E ab of the control volumes a and b.
  • ⁇ ab is calculated by interpolation using ⁇ a and ⁇ b . That is, ⁇ ab can be expressed by Expression (3-24).
  • x a and x b are coordinate vectors of control points C a and C b that are representative points of the control volumes a and b, respectively.
  • interpolation functions for calculating ⁇ ab are conventionally known. Which interpolation function (also referred to as a scheme) is used depends on the purpose of the simulation, the type of continuum to be used, a numerical calculation method, and the like, but any of them may be used. Here are some typical schemes.
  • the sigma with a tilde in the equations (3-30) and (3-31) is a continuum flow rate S ab ⁇ (n ab ⁇ u) when ⁇ b flows from the control volume b adjacent to the control volume a.
  • S ab ⁇ (n ab ⁇ u) is a continuum flow rate S ab ⁇ (n ab ⁇ u) when ⁇ b flows from the control volume b adjacent to the control volume a.
  • This is an equation indicating that ⁇ a flows out at the same flow rate as ab ). That is, the balance of the flow rate of inflow / outflow with respect to the control volume a is a formula. This is an example showing the effectiveness of the present invention.
  • control volume is divided into four control volumes a1, a2, a3, and a4 in the flow direction.
  • the velocities at the respective boundary surfaces are assumed to be 1.0, 1.1, 0.9, 0.8, and 1.0 sequentially from the upstream boundary surface of the control volume a1.
  • this velocity distribution is substituted into the continuous equation in each control volume, ⁇ 0.1, 0.2, 0 in order from the control volume a1 as shown in FIG. .1 and -0.2, both of which are not 0.
  • u i is the rate at the interface between the control volume a i-1 of the upstream side of the control volume a i, phi i, the physical quantity of the control volume a i (temperature), phi i ', the front
  • the physical quantity ⁇ i ⁇ 1 of the control volume at the time step is the physical quantity of the upstream control volume a i ⁇ 1 .
  • FIG. 3 shows an example of a conventional calculation result when the equation (2-20 ′) is used
  • FIG. 4 shows an example of the calculation result of the present invention when the equation (3-30 ′) is used.
  • the calculation result should be uniformly 100 ° C.
  • the distribution does not become uniform in the past, but in the calculation result example of the present invention, Uniformly 100 ° C.
  • ⁇ ab is represented by Formula (3-32). That is, an interpolation function using grad (gradient vector) of ⁇ a and ⁇ b is used.
  • Equation (3-28) the sigma term in equations (3-22) and (3-23) is the sum of the advection terms for all control volumes b adjacent to the control volume a.
  • FIG. 5 is a schematic block diagram showing the configuration of the fluid simulation apparatus 10 according to the first embodiment of the present invention.
  • the fluid simulation apparatus 10 according to the present embodiment is a simulation apparatus that simulates an incompressible fluid by a finite volume method and calculates a speed and a pressure in a simulation region. That is, in this embodiment, a simulation apparatus that calculates a solution of the Navier-Stokes equation, which is a transport equation related to momentum, will be described.
  • the fluid simulation apparatus 10 includes a cell configuration setting unit 101, a cell configuration storage unit 102, an initial condition setting unit 103, a velocity field storage unit 104, a velocity field acquisition unit 105, a velocity / pressure calculation unit 106, a pressure
  • the field storage unit 107 is included.
  • the cell configuration setting unit 101 sets information on a cell obtained by dividing an area (space) to be simulated, and stores the information in the cell configuration storage unit 102.
  • the cell configuration storage unit 102 stores information related to the cells set by the cell configuration setting unit 101.
  • the information regarding a cell contains the cell information which is the information of each cell, and the information between cells which is the information between adjacent cells.
  • the initial condition setting unit 103 sets the velocity vector of each cell as an initial condition and stores it in the velocity field storage unit 104.
  • the speed field storage unit 104 stores the speed vector of each cell set by the initial condition setting unit 103.
  • the velocity field storage unit 104 also stores the velocity vector of each cell calculated by the velocity / pressure calculation unit 106 for each time step.
  • the velocity field acquisition unit 105 reads the velocity vector of each cell at the time step n ⁇ 1 from the velocity field storage unit 104 when the velocity / pressure calculation unit 106 calculates a solution regarding the time step n.
  • the velocity / pressure calculation unit 106 uses the velocity vector of each cell at time step n-1 read by the velocity field acquisition unit 105 from the velocity field storage unit 104 and the cell information stored in the cell configuration storage unit 102. The velocity vector and pressure of each cell at the target time step n are calculated. The velocity / pressure calculation unit 106 stores the calculated velocity vector of each cell in the velocity field storage unit 104 and stores the pressure in the pressure field storage unit 107. The velocity / pressure calculation unit 106 discretely calculates the Navier-Stokes equation for the incompressible fluid and the continuous equation for the incompressible fluid (also referred to as a mass conservation equation) when calculating the velocity vector and the pressure. Use the discretized equation.
  • the pressure field storage unit 107 stores the pressure of each cell calculated by the speed / pressure calculation unit 106 for each time step.
  • FIG. 6 is a conceptual diagram illustrating a cell in the present embodiment.
  • each of the regions denoted by reference symbols R 1 to R 15 is the 1st to 15th cells.
  • each cell is a polyhedron.
  • the cells are arranged so as to fill the region to be simulated without any gaps.
  • the boundary surface E 12 is a boundary surface between the cell R 1 and R 2.
  • Unit normal vector n 12 is the unit normal vector of the boundary surface E 12.
  • the distance alpha 12 is when the one of the distance from the representative point C 1 to C 2, is the distance from the intersection of the boundary surface E 12 of the line connecting the points C 1 and C 2 to the representative point C 1.
  • FIG. 7 is a diagram illustrating an example of cell information stored in the cell configuration storage unit 102.
  • the cell information includes cell IDs, representative point coordinates, and volumes of all cells.
  • the cell information of the cell with the cell ID “0001” is the cell ID “0001”, the representative point coordinates “15, 225, 214”, and the volume “24.3”.
  • the cell information of the cell with the cell ID “0002” is the cell ID “0002”, the representative point coordinates “15, 225, 214”, the volume “24.3”, and the cell ID “2453”.
  • Cell information is a cell ID “2453”, representative point coordinates “215, 25, 44”, and volume “31.2”.
  • the three numerical values of the representative point coordinates are the x coordinate, the y coordinate, and the z coordinate.
  • X n indicates the representative point coordinates of the cell with the cell ID n
  • V n indicates the volume of the cell with the cell ID n.
  • FIG. 8 is a diagram illustrating an example of inter-cell information stored in the cell configuration storage unit 102.
  • the cell-to-cell information includes the cell ID of the reference cell, the cell ID of the cell adjacent to the cell, the area of the boundary surface between these two cells, the unit normal vector of the boundary surface, these It includes the distance ⁇ from the intersection of the straight line passing through the representative point of the two cells and the boundary surface to the representative point of the reference cell.
  • the unit normal vector is a unit vector having a length of “1”. Further, the distance is a distance when the distance between the representative points of the two cells was evaluated as "1", corresponding to alpha 12 in FIG.
  • the cell-to-cell information between the cell with the cell ID “0001” and the cell with the cell ID “0005” is cell ID “0001”, cell ID “0005”, and area “4.2458”.
  • cell-to-cell information regarding the cell between the cell with cell ID “0001” and the cell with cell ID “0102” is cell ID “0001”, cell ID “0102”, area “3.133352”, and unit normal.
  • the inter-cell information regarding the cell between the cell with the cell ID “2453” and the cell with the cell ID “2452” is the vector “0.584276, ⁇ 0.617123, 0.527049”, the distance “0.554878”, The cell ID is “2453”, the cell ID is “2452”, the area is “3.3215”, the unit normal vector is “0.335786, 0.814266, 0.473517”, and the distance is “0.515486”.
  • FIG. 9 is a diagram illustrating an example of speed field information stored in the speed field storage unit 104.
  • the velocity field storage unit 104 stores velocity field information at an initial time step set as an initial condition, and velocity field information at each time step calculated by the velocity / pressure calculation unit 106.
  • the example shown in FIG. 9 is velocity field information in one time step among them.
  • the velocity field information in one time step includes the cell ID of each cell and the velocity vector of the cell.
  • the example of the velocity field information shown in FIG. 9 includes a cell ID “0001”, a velocity vector “1.221354, 3.255546, 0.2513356”, a cell ID “0002”, and a velocity vector “1” of the cell. .240134, 3.15224, 0.662158 ", the cell ID" 2453 ", and the velocity vector" 3.22211, 3.00125, 1.658855 "of the cell.
  • FIG. 10 is a diagram illustrating an example of pressure field information stored in the pressure field storage unit 107.
  • the pressure field storage unit 107 stores the pressure field information at each time step calculated by the speed / pressure calculation unit 106.
  • the example shown in FIG. 10 is pressure field information at one of them.
  • the pressure field information in one time step includes the cell ID of each cell and the pressure of the cell.
  • the example of the pressure field information shown in FIG. 10 includes a cell ID “0001”, a pressure “1.2568” of the cell, a cell ID “0002”, a pressure “1.22325” of the cell, and a cell ID “2453”. And the pressure of the cell “1.56846”.
  • FIG. 11 is a flowchart for explaining the operation of the fluid simulation apparatus 10 in the present embodiment.
  • the cell configuration setting unit 101 causes the cell configuration storage unit 102 to store information (cell information and inter-cell information) related to cells acquired from the outside.
  • the initial condition setting unit 103 stores the velocity field information of the initial time step, which is the initial condition acquired from the outside, in the velocity field storage unit 104 (S1).
  • the speed field acquisition unit 105 sets the target time step to an initial value “1” (S2).
  • the speed field acquisition unit 105 acquires speed field information one step before the target time step from the speed field storage unit 104 (S3).
  • the target time step is “1”
  • the speed field acquisition unit 105 sets the initial time step one step before, and acquires the speed field information of the initial time step from the speed field storage unit 104.
  • the velocity / pressure calculation unit 106 calculates each element of the matrix A and the vector b using the velocity field information acquired by the velocity field acquisition unit 105 in step S3.
  • the vector ⁇ is a vector having the elements of the velocity vectors and pressures of all cells in the target time step as elements. That is, the vector ⁇ has the number of cells ⁇ 4 elements.
  • the simultaneous linear equations are composed of an equation obtained by discretizing the Navier-Stokes equations in each axial direction for each cell, and an equation obtained by discretizing an equation of continuous incompressible fluid for each cell. That is, the simultaneous linear equations are composed of the same number of expressions as the number of elements of the vector ⁇ , so that a solution can be calculated. Details of the method of calculating the matrix A and the vector b in step S4 will be described later.
  • the speed / pressure calculation unit 106 stores the speed vector among the calculated elements of the vector ⁇ in the speed field storage unit 104 as the speed vector of the target time step. Similarly, the speed / pressure calculation unit 106 stores the pressure in the calculated element of the vector ⁇ in the pressure field storage unit 107 as the pressure of the target time step (S5).
  • the speed / pressure calculation unit 106 determines whether or not an end condition is satisfied (S6).
  • Examples of the termination condition include that the target time step matches a preset value.
  • the process proceeds to step S7, where the speed / pressure calculation unit 106 increases the target time step by “1” (S7), and returns to step S3. . Thereby, processing for the next time step is performed. If it is determined in step S6 that the termination condition is satisfied (S6-Yes), the process is terminated.
  • the Navier-Stokes equation of the incompressible fluid in the i-th axis direction is expressed as in Equation (4-1).
  • the continuous equation of the incompressible fluid is represented by the equation (4-2).
  • an orthogonal coordinate system (Cartesian coordinate system) is used as the coordinate system.
  • Einstein's sum rules are used in the following formulas including formulas (4-1) and (4-2).
  • Equation (4-1) t is time.
  • p is the pressure.
  • ⁇ is the density.
  • ⁇ is a kinematic viscosity coefficient.
  • Equation (4-1) the second term on the left side is the advection term.
  • equation (4-1) is an equation for conservation of momentum in the i-th axis direction and because it is an incompressible fluid, the density ⁇ of the fluid is constant, so advection in equation (4-1) In the term, the advection is the velocity u i in the i-th axis direction.
  • the sum of the second term (viscosity term) and the third term (external force) on the right side of Equation (4-1) will be denoted as -F (u), This will be described using Expression (4-3) in which the term is moved to the left side.
  • the Navier-Stokes equation is volume-integrated with a control volume, and a discrete expression is solved.
  • an equation (4-4) obtained by multiplying both sides of the equation (4-2), which is a continuous equation, by the velocity u i is added to the Navier-Stokes equation to obtain an equation (4
  • Use a modified advection term called a conservative type as in -5) is a transformation on the premise that the velocity vector satisfies the continuous equation (4-2).
  • Equation (4-6) it is assumed that even if the velocity vector u is substituted into the left side of the continuous equation, it does not become “0”, and there is an error D as shown in Equation (4-6).
  • the size of the error D may not be negligible when the calculation cannot be performed until the speed is sufficiently converged when the speed is obtained by iterative calculation due to problems such as calculation capacity. Such an error is likely to occur when the cell shape is long in a specific axial direction.
  • Equation (4-9) volume integrated with the control volume
  • equation (4-10) is obtained.
  • equation (4-11) is obtained by volume integration with the control volume.
  • V is the volume of the control volume
  • S is the surface of the control volume.
  • n and u in bold font in the second term on the left side of Equation (4-10) and the first term on the left side of Equation (4-11) are the unit normal vector and velocity vector of the surface S of the control volume, respectively. is there.
  • the unit normal vector n is a vector that is perpendicular to the surface S of the control volume, faces outward, and has a length of “1”.
  • equations (4-10) and (4-11) are expressed by, for example, the control volume V 1 and the boundary surface E 1k (k ⁇ N 1 ) of the cell R 1 in FIG. 6, and the set N 1 is ⁇ 2, 3 4, 11, 12, 13, 14) are discretized with respect to the area S 1k and time ⁇ t, and converted into approximate equations by algebraic equations, equations (4-12) and (4-13) are obtained.
  • the boundary surface E 1k is a boundary surface between the cell R 1 and the cell R k .
  • the set N 1 is a set of cell IDs of cells adjacent to the cell R 1 .
  • ⁇ t is the width of the time step.
  • the area S 1k is an area of the boundary surface between the cell R 1 and the cell R k .
  • D 1 is the error in the continuous equation in cell R 1 .
  • u 1i is the i-th axial component of the velocity in cell R 1 .
  • u 1i (t ⁇ t) is the i-th axial component of the velocity in cell R 1 at the time step before the target time step.
  • u 1ki is the i-th axial component of the velocity at the boundary surface E 1k .
  • n 1k is a unit normal vector facing the cell R k in the boundary surface E 1k .
  • u 1k is a velocity vector at the boundary surface E 1k .
  • p k is the pressure in cell R k .
  • this equation (4-14) is one of the linear equations constituting the simultaneous linear equations.
  • the second term includes multiplication of unknown speeds and is not a linear equation. Therefore, approximation is performed in which the velocity vector of the previous time step is used as the velocity vector u 1k in the second term. For the discretization term of F (u), approximation is performed by using the velocity vector of the previous time step in order to simplify the explanation.
  • a primary upwind scheme is used.
  • the speed in the windward cell is used as the speed of the boundary surface.
  • the unit normal vector n 1k is an outward vector when viewed from the cell R1. Therefore, when the inner product of the unit normal vector n 1k and the velocity vector v 1k of the boundary surface is negative, the cell R k is windward at the boundary surface E 1k , and when positive, the cell R 1 is windward Above.
  • This primary upwind scheme is represented by the formula (4-15). Note that when the inner product is negative, it flows into the cell R 1 , and when it is positive, it is also called out from the cell R 1 .
  • the velocity vector u 1k of the boundary surface E 1k is set as a velocity vector u 1k ′ obtained by interpolation from the velocity vector of each cell in the previous time step.
  • the velocity vectors u 1 ′ and u 1 ′ in the previous time step are used. Is used.
  • Expression (4-14) becomes Expression (4-17).
  • the set N 1 'among the cells adjacent to the cell R 1 is a set of the cell number of the cell that flows into the cell R 1 and element.
  • the equation (4-17) Transforming equation (4-18), and the speed of the i-axis component u 1i cell R 1, and i-axis component u ki of the velocity of each cell R k adjacent cell R 1 Is a linear equation with the unknown pressure p.
  • an expression of a continuous equation (4-2) is also control volume V 1 of the cell R 1, a boundary surface E 1k (k ⁇ N 1, the set N 1 is ⁇ 2,3,4,11,12 , discretizing the area S 1k of 13,14 a ⁇ ), is converted into the approximate expression by algebraic equations, equation (4-19) is obtained.
  • equation (4-19) becomes the equation (4-20), and each component of the velocity of the cell R 1 and each component of the velocity of the adjacent cell R k are obtained. It becomes a linear equation with an unknown number.
  • the advection term is transformed into a conservative form by adding a continuous equation.
  • the continuous equation has an error D, and the error D is discretized including a term for correcting the error D. Therefore, even if the speed of each time step does not completely satisfy the continuous equation, the influence of the error can be suppressed, so that an excellent simulation result can be obtained.
  • the speed / pressure calculation unit 106 calculates the solution for each time step
  • an iterative method such as the conjugate gradient method
  • an initial value of the solution is required, but the solution of the previous time step is used as the initial value. May be used.
  • the velocity field acquisition unit 105 acquires the pressure of the previous time step from the pressure field storage unit 107 and inputs the pressure to the velocity / pressure calculation unit 106.
  • FIG. 12 is a schematic block diagram showing a configuration of a fluid simulation apparatus 10a according to the second embodiment of the present invention.
  • the fluid simulation apparatus 10a according to the present embodiment is a simulation apparatus that simulates a compressible fluid by a finite volume method and calculates speed, pressure, and density in a simulation region.
  • the fluid simulation apparatus 10a includes a cell configuration setting unit 101, a cell configuration storage unit 102, an initial condition setting unit 103a, a density / velocity field storage unit 104a, a density / velocity field acquisition unit 105a, and a solution calculation unit 106a.
  • the pressure field storage unit 107 is included.
  • the same reference numerals (101, 102, 107) are assigned to portions corresponding to the respective portions in FIG.
  • Initial condition setting unit 103a sets the velocity vector and density of each cell as initial conditions, and stores them in the density / velocity field storage unit 104a.
  • the density / velocity field storage unit 104a stores the velocity vector and density of each cell set by the initial condition setting unit 103a.
  • the density / velocity field storage unit 104a also stores the velocity vector and density of each cell calculated by the solution calculation unit 106a for each time step.
  • the solution calculation unit 106a calculates a solution related to the time step n
  • the density / speed field acquisition unit 105a reads the velocity vector and density of each cell at the time step n-1 from the density / speed field storage unit 104a.
  • the solution calculation unit 106a includes the velocity vector and density of each cell at time step n-1 read from the density / velocity field storage unit 104a by the density / velocity field acquisition unit 105a, cell information stored in the cell configuration storage unit 102, and Is used to calculate the velocity vector, pressure, and density of each cell at the target time step n.
  • the solution calculation unit 106 a stores the calculated velocity vector and density of each cell in the density / velocity field storage unit 104 and stores the pressure in the pressure field storage unit 107.
  • the solution calculation unit 106a discretizes a discretized equation obtained by discretizing the Navier-Stokes equation of the compressible fluid, a discretized equation obtained by discretizing the continuous equation of the compressible fluid, and a state equation representing the relationship between density and pressure.
  • the simultaneous linear equations including the discretized equations are solved, and the velocity vector, pressure, and density of each cell are calculated.
  • the solution calculation unit 106a uses a known simultaneous linear equation solving method such as a conjugate gradient method.
  • R is a gas constant
  • T temperature (thermodynamic temperature).
  • the temperature T is constant.
  • the temperature T is not constant and may be known.
  • a discretization formula using the temperature T as a variable may be included in the simultaneous linear equations described above.
  • the discretization equation obtained by discretizing the Navier-Stokes equation of the compressible fluid will be described.
  • the Navier-Stokes equation of the compressible fluid in the i-th axis direction is expressed as the following equation (4-21).
  • the continuous equation of the incompressible fluid is represented by the equation (4-22).
  • an orthogonal coordinate system (Cartesian coordinate system) is used as the coordinate system.
  • ⁇ is a viscosity coefficient.
  • equation (4-22) which is a continuous equation.
  • equation (4-22) the right side of equation (4-22) and the error D, is multiplied by u i both sides, equation (4-23) is obtained. Then, by adding this equation (4-23) to the Navier-Stokes equation (4-21), the equation (4-24) is obtained.
  • This equation (4-25) is the same as the equation used in the conventional finite volume method except that there is a third term on the left side.
  • the error in the conservative advection term is corrected in the third term.
  • equation (4-25) is volume integrated with the control volume
  • equation (4-26) is obtained.
  • equation (4-27) is obtained by volume integration with the control volume.
  • V is the volume of the control volume
  • S is the surface of the control volume.
  • n and u in bold font in the second term on the left side of the equation (4-26) are the unit normal vector and the velocity vector of the surface S of the control volume, respectively.
  • the unit normal vector n is a vector that is perpendicular to the surface S of the control volume, faces outward, and has a length of “1”.
  • equations (4-26) and (4-27) are expressed by, for example, the control volume V 1 and the boundary surface E 1k (k ⁇ N 1 ) of the cell R 1 in FIG. 6, and the set N 1 is ⁇ 2, 3 4, 11, 12, 13, 14) are discretized with respect to the area S 1k and time ⁇ t, and converted to approximate equations using algebraic equations, equations (4-28) and (4-29) are obtained. .
  • equation (4-29) Substituting equation (4-29) into the third term of equation (4-28) yields equation (4-30). Similar to the first embodiment, a linear upwind scheme or the like is applied to this equation (4-30) to form a linear equation with the velocity vector and pressure as unknowns.
  • This formula (4-30) is a formula related to the cell R1, but a similar formula can be established for each of the other cells.
  • the solution calculation unit 106a sets each unknown coefficient in the equation (4-30) for all cells as an element of the matrix A and the right side as an element of the vector b.
  • the advection term is transformed into a conservative type by adding a continuous equation.
  • the continuous equation has an error D
  • the error D is discretized including a term for correcting the error D. Therefore, even if the speed of each time step does not completely satisfy the continuous equation, the influence of the error can be suppressed, so that an excellent simulation result can be obtained.
  • the Navier-Stokes equation is given as an equation including the advection term, and the fluid simulation apparatuses 10 and 10a that solve this are described.
  • a continuum thermal energy storage equation 10b that solves this by a finite volume method will be described as an equation including the advection term as an equation including the advection term.
  • the continuum simulation apparatus 10b of this embodiment uses the value input from the outside for the velocity vector of each time step in the area
  • FIG. 13 is a schematic block diagram showing the configuration of the continuum simulation apparatus 1 according to the third embodiment of the present invention.
  • the continuum simulation apparatus 10b in this embodiment includes a cell configuration setting unit 101, a cell configuration storage unit 102, a velocity field setting unit 103b, a velocity field storage unit 104b, a velocity field acquisition unit 105b, a solution calculation unit 106b, and a solution storage unit 107b. It is comprised including.
  • portions corresponding to the respective portions in FIG. 5 are denoted by the same reference numerals (101, 102), and description thereof is omitted.
  • the velocity field setting unit 103b stores the velocity vector of each cell input from the outside in the velocity field storage unit 104b.
  • velocity vectors at all time steps to be simulated are input as velocity vectors.
  • the speed field storage unit 104b stores a speed vector of each cell.
  • the solution calculation unit 106b calculates a solution related to the time step n
  • the velocity field acquisition unit 105b reads the velocity vector of each cell at the time step n-1 from the velocity field storage unit 104b.
  • the solution calculation unit 106b includes a velocity vector of each cell at time step n-1 read by the velocity field acquisition unit 105b from the velocity field storage unit 104b, cell information stored in the cell configuration storage unit 102, and a solution storage unit 107b. Using the stored solution at time step n ⁇ 1, the solution of the thermal energy conservation equation at target time step n is calculated. The solution calculation unit 106b stores the calculated solution in the solution storage unit 107b as a solution for the target time step n.
  • the solution calculating unit 106b uses a discretized equation obtained by discretizing the energy conservation equation shown in Equation (4-31) when calculating the solution. As shown in the equation (4-31), the energy conservation equation has an advection term in the second term on the left side. Since the density ⁇ is not constant, the solution calculation unit 106b uses a discretized equation obtained by discretizing the equation (4-31) in the same manner as the discretization of the Navier-Stokes equation in the second embodiment.
  • the advection diffusion equation describing the advection diffusion phenomenon of the substance in the continuum also has an advection term in which the concentration of the substance is advected as shown in the equation (4-32). For this reason, you may make it solve the discretization equation which discretized this advection diffusion equation similarly to the formula of energy conservation.
  • the stress field equation describing the stress field in the continuum also has an advection term in which the deformation speed of the continuum is advected as shown in Expression (4-33). Therefore, a discretized equation obtained by discretizing the stress field equation in the same manner as the Navier-Stokes equation may be solved.
  • a material derivative including an advection term is used for an equation representing a phenomenon involving movement and deformation of a continuum. You may make it solve the discretization equation which discretized this equation like each above-mentioned embodiment.
  • a fourth embodiment of the present invention will be described with reference to the drawings.
  • a fluid simulation apparatus 10c that solves an equation including a plurality of advection terms by a finite volume method will be described.
  • a fluid simulation apparatus 10c that simultaneously solves the Navier-Stokes equation and the thermal energy conservation equation will be described as an example of a simulation apparatus that simultaneously solves a plurality of equations having advection terms.
  • FIG. 14 is a schematic block diagram showing a configuration of a fluid simulation apparatus 10c according to the fourth embodiment of the present invention.
  • the fluid simulation apparatus 10c in this embodiment includes a cell configuration setting unit 101, a cell configuration storage unit 102, an initial condition setting unit 103, a velocity field storage unit 104, a velocity field acquisition unit 105, a solution calculation unit 106c, and a solution storage unit 107c. Consists of including. In the figure, portions corresponding to the respective portions in FIG. 5 are denoted by the same reference numerals (101 to 105), and description thereof is omitted.
  • the solution calculation unit 106c includes the velocity vector of each cell at time step n-1 read by the velocity field acquisition unit 105 from the velocity field storage unit 104, cell information stored in the cell configuration storage unit 102, and the solution storage unit 107c. Using the stored solution at time step n-1, the solution of the Navier-Stokes equation, the continuity equation, and the thermal energy conservation equation at the target time step n is calculated. The solution calculation unit 106c stores the calculated solution in the solution storage unit 107c as the solution of the target time step n. The solution calculation unit 106c stores the velocity in the calculated solution in the velocity field storage unit 104.
  • the solution calculation unit 106c generates a discretized equation obtained by discretizing the Navier-Stokes equation and the continuity equation in the first embodiment and a discretized equation obtained by discretizing the thermal energy conservation equation in the third embodiment. Use.
  • the Navier-Stokes equation and the thermal energy conservation equation are solved simultaneously, but other combinations may be used. .
  • FIG. 15 is a contour diagram (contour map) showing a velocity field to be simulated. That is, FIG. 15 is a contour diagram of the magnitude of the velocity vector of each cell that is stored in the velocity field storage unit 104b by the velocity field setting unit 103b. The velocity field was calculated on the assumption that there was no change over time.
  • the horizontal axis is the x-axis
  • the vertical axis is the y-axis
  • the magnitude of the velocity field is 4 m in the x-axis direction and 1 m in the y-axis direction.
  • FIG. 16 is a contour diagram showing the magnitude of velocity divergence in the velocity field shown in FIG. Also in FIG. 16, the horizontal axis is the x-axis and the vertical axis is the y-axis.
  • the divergence of velocity coincides with the left side of the continuous equation shown in Equation (2-3). That is, if the velocity field in FIG. 15 completely satisfies the continuity equation, the contour diagram in FIG. 16 should be uniformly “0”. However, in the velocity field shown in FIG. 15, since there is a region of “err” having a velocity of 1.2 m / s, “0” is not obtained at the left and right ends of the region.
  • FIG. 17 is a contour diagram showing an example in which the temperature distribution of the velocity field shown in FIG. 15 is calculated using a conventional method.
  • the initial value of the temperature distribution is uniformly 100 ° C.
  • the temperature distribution shown in FIG. 17 is a temperature distribution in a steady state.
  • the horizontal axis is the x-axis
  • the vertical axis is the y-axis
  • the contour lines are lines obtained by dividing 83 ° C. to 107 ° C. into 13 stages.
  • the initial value of the temperature distribution is uniformly 100 ° C., but exceeds 100 ° C. downstream from the left and right ends of the shape of “err”. Or, an area below is created. This is because the divergence of the velocity is not “0” at the left and right ends of the shape “err”, and a pseudo heat sink / source is generated at these locations.
  • FIG. 18 is a contour diagram showing an example when the temperature distribution of the velocity field shown in FIG. 15 is calculated by the continuum simulation apparatus 10b of the third embodiment.
  • the initial value of the temperature distribution is uniformly 100 ° C., which is the temperature distribution in the steady state.
  • the horizontal axis is the x-axis
  • the vertical axis is the y-axis
  • the contour lines are lines obtained by dividing 83 ° C. to 107 ° C. into 13 stages.
  • the temperature distribution calculated by the continuum simulation apparatus 10b is uniformly 100 ° C. because no pseudo heat sink / source is generated. If the fluid at 100 ° C. is flowing uniformly, the temperature is uniformly 100 ° C. even if time passes, and it can be seen that a good simulation result is obtained by the continuum simulation apparatus 10b.
  • the program for realizing the function of each part in FIG. 5, FIG. 12, FIG. 13, or FIG. 14 or a part of the function is recorded on a computer-readable recording medium, and the program recorded on this recording medium is recorded.
  • the “computer system” includes an OS and hardware such as peripheral devices.
  • the “computer-readable recording medium” means a storage device such as a flexible disk, a magneto-optical disk, a portable medium such as a ROM and a CD-ROM, and a hard disk incorporated in a computer system. Furthermore, the “computer-readable recording medium” dynamically holds a program for a short time like a communication line when transmitting a program via a network such as the Internet or a communication line such as a telephone line. In this case, a volatile memory in a computer system serving as a server or a client in that case, and a program that holds a program for a certain period of time are also included.
  • the program may be a program for realizing a part of the functions described above, and may be a program capable of realizing the functions described above in combination with a program already recorded in a computer system.
  • the simulation apparatus of the present invention can be used for a generation source estimation apparatus as described in, for example, JP 2012-132824 A.
  • This source estimation device obtains the position information and measured concentration information of the observation devices installed at multiple locations, and estimates the release point from a plurality of preset virtual release points based on this information is doing.
  • the virtual emission point is the generation source, and what concentration is measured by each observer when the unit amount is released (the influence function in the document) ) And the estimation result is used.
  • the continuum simulation apparatus that solves the advection diffusion equation shown in the third embodiment can be used for the estimation of the measurement value of the concentration, regardless of the setting method of the virtual discharge point.
  • a fluid simulation apparatus that simultaneously solves the Navier-Stokes equation and the advection diffusion equation can be used.
  • the generation source may be gas or liquid.

Abstract

　優れたシミュレーション結果を得ることができるシミュレーション装置を提供すること。　移流項を含む方程式の解を算出するシミュレーション装置であって、速度場を取得する速度場取得部と、方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、取得した速度場を適用して、離散化方程式における解を算出する解算出部とを具備し、上述の２つの操作は、移流項を保存型に変形することと、移流項の保存型の誤差を修正する項を含めることであることを特徴とする。

Description

シミュレーション装置、シミュレーション方法およびプログラム
　本発明は、シミュレーション装置、シミュレーション方法およびプログラムに関する。
　従来、流体などの連続体の数値シミュレーションに用いられる方法として、有限体積法がある。この有限体積法では、シミュレーションする領域をセルあるいはコントロールボリュームと呼ばれる領域に分割しておく。そして、各セルについて、支配方程式を体積積分し、さらに離散化した式を用いて、解を算出する。通常、流体をシミュレーションするときは、支配方程式であるナビエ・ストークス方程式を体積積分した後、離散化した式を用いる。この体積積分をする際に、移流項と呼ばれる項については、保存型と呼ばれる形式に変換しておく。これにより、移流項の体積積分の結果は表面積分となる（例えば、非特許文献１のＰ．２７参照）。なお、この非特許文献１においては、移流項のことを対流項と呼んでいる。また、非圧縮性流れに対する完全保存型については、離散化の差分スキームについて研究がされている（非特許文献２、３）。
　上述の移流項について説明する。図１９に示すように連続体中に検査体積Ｖ（ｔ）を取る。連続体が移動、変形するときは、検査体積Ｖ（ｔ）も連続体の移動、変形に伴って移動、変形する。図１９は、時刻ｔ→ｔ＋Δｔと、Δｔだけ時間進行するとき、原点Ｏ’（ｔ）の位置にあった検査体積Ｖ（ｔ）が、Ｏ’（ｔ＋Δｔ）の位置へ移動し、検査体積Ｖ（ｔ＋Δｔ）、表面Ｓ（ｔ＋Δｔ）の状態へ、移動、変形する様子を示す。この移動、変形は、Ｖ（ｔ）を囲む表面Ｓ（ｔ）が速度ベクトルｕｓ（ｔ）の速度で移動することにより、生じている。
　Ｖ→０、Δｔ→０の極限操作により、連続体を表す偏微分方程式が得られる。この偏微分方程式の記述方法には、ラグランジュ座標系による記述方法と、オイラー座標系による記述方法がある。ラグランジュ座標系は、図１９の原点Ｏ’（ｔ）のように、連続体とともに移動する原点を持つ座標系である。オイラー座標系は、図１９の原点Ｏのように、連続体とは別に設定された座標系である。オイラー座標系では、原点Ｏは、連続体とともに移動しない。
　連続体として表される物理量をφとし、φの偏微分方程式が、ラグランジュ座標系で以下の式（１－１）のように表されるものとする。式（１－１）をオイラー座標系で表すと、式（１－２）となる。
Figure JPOXMLDOC01-appb-M000001
　ここで、ｕ、ｖ、ｗは、連続体内の座標（ｘ、ｙ、ｚ）の位置の速度ベクトルｕのｘ、ｙ、ｚ成分であり、速度ベクトルｕは、式（１－４）のように表される。
Figure JPOXMLDOC01-appb-M000002
　式（１－３）の左辺にあるＤ／Ｄｔは、全微分あるいは物質微分と呼ばれる。式（１－３）の右辺第２項から第４項（式（１－５））が、移流項と呼ばれる。
Figure JPOXMLDOC01-appb-M000003
　添え字とアインシュタインの総和規則を用いると、式（１－２）、（１－３）は、式（１－６）、（１－７）となる。ここで添え字ｉ、ｊは、ｉ＝１、２、３、ｊ＝１、２、３のインデックスを示し、１，２，３→ｘ、ｙ、ｚに対応する。同じ添え字が２度表れたときは、総和規則に従う。
Figure JPOXMLDOC01-appb-M000004
　式（１－６）に示されている物質微分Ｄ／Ｄｔは、次に示すように、連続体として表現される全ての物理量の、時間進行による連続体の移動、変形を表す項として、オイラー座標系で表された偏微分方程式の中に現れる。式（１－３）、（１－７）で示したように、物質微分は、移流項を含んでいるので、物質微分が現れる偏微分方程式は、移流項を含んでいる。
　物質微分を含む偏微分方程式の例をいくつか以下に示す。
１）流体力学（ナビエ・ストークス方程式）
　流体力学の支配方程式であるナビエ・ストークス方程式は、式（１－８）の様に表される。式（１－８）において、ｔ：時間、ｘｉ（ｉ＝１、２、３）：カーテシアン系における座標、ρ：流体の密度、μ：流体の粘性係数、ｕｉ（ｉ＝１、２、３）：流体の速度成分、ｐ：圧力である。また、添え字ｉ（ｉ=１、２、３）、ｊ（ｊ＝１、２、３）はカーテシアン座標系における各方向成分を示す。式（１－８）では、左辺に物質微分が現れている。すなわち、移流項を含んでいる。
Figure JPOXMLDOC01-appb-M000005
２）応力場の方程式
　オイラー座標系における応力場の方程式は、式（１－９）である。式（１－９）において、ρ：連続体の密度、ｕｉ（ｉ＝１、２、３）：連続体の変形速度、σｉｊ（ｉ、ｊ＝１、２、３）：連続体の内部応力、ｆｉ＝連続体に作用する外力（例えば、重力）である。式（１－９）でも、左辺に物質微分が現れている。すなわち、移流項を含んでいる。
Figure JPOXMLDOC01-appb-M000006
３）エネルギー保存の方程式
　エネルギー保存則は、熱エネルギー保存と運動エネルギー保存とに分けられる。運動エネルギーの保存は、式（１－９）と同じであるので、ここでは熱エネルギー保存の方程式を式（１－１０）に示す。式（１－１０）において、Ｕ：連続体の内部エネルギー、ｑｊ（ｊ＝１、２、３）：熱流束ベクトル成分、ρ：連続体の密度、γ：熱源（熱エネルギーのソース／シンク）、σｉｊ（ｉ、ｊ＝１、２、３）：連続体の内部応力テンソル、Ｄｉｊ（ｉ、ｊ＝１、２、３）：連続体の変形速度テンソルである。式（１－１０）でも、左辺に物質微分が現れている。すなわち、移流項を含んでいる。
Figure JPOXMLDOC01-appb-M000007
４）移流拡散方程式
　ある物質Ｃの連続体内での移流拡散現象は、以下の式（１－１１）である。式（１－１１）において、Ｃ：物質Ｃの濃度、μｃ：物質Ｃの拡散係数、ｑｃ：物質Ｃのソース・シンク、ρ：連続体の密度、ｕｊ（ｊ＝１、２、３）：連続体の変形速度成分である。式（１－１１）でも、左辺に物質微分が現れている。すなわち、移流項を含んでいる。
Figure JPOXMLDOC01-appb-M000008
Ｈ．Ｋ．Ｖｅｒｓｔｅｅｇ、Ｗ．Ｍａｌａｌａｓｅｋｅｒａ（原著）、松下洋介、齋藤泰洋、青木秀之、三浦隆利（共訳）、"数値流体力学［第２版］"、（Ａｎ　Ｉｎｔｒｏｄｕｃｔｉｏｎ　ｔｏ　Ｃｏｍｐｕｔａｔｉｏｎａｌ　Ｆｌｕｉｄ　Ｄｙｎａｍｉｃｓ　Ｓｅｃｏｎｄ　Ｅｄｉｔｉｏｎ）、第２版、森北出版株式会社、２０１１年５月３０日． 森西洋平、〔竜門賞受賞記念解説〕非圧縮性流体解析のための完全保存形高次精度差分スキーム、ながれ １９（２０００）、Ｐ．１６４－１７０． Ｍｏｒｉｎｉｓｈｉ，Ｙ．，Ｌｕｎｄ，Ｔ．Ｓ．，Ｖａｓｉｌｙｅｖ，Ｏ．Ｖ．，Ｍｏｉｎ，Ｐ．，"Ｆｕｌｌｙ　ｃｏｎｓｅｒｖａｔｉｖｅ　ｈｉｇｈｅｒ　ｏｒｄｅｒ　ｆｉｎｉｔｅ　ｄｉｆｆｅｒｅｎｃｅ　ｓｃｈｅｍｅｓ　ｆｏｒ　ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　ｆｌｏｗ，Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．，Ｖｏｌ．１４３（１９９８），９０－１２４．
　上述したように移流項は、ナビエ・ストークス方程式だけでなく、物質微分を含む方程式（以下、輸送方程式という）に現れる項である。本発明が解決しようとする課題を説明する前に、移流項を保存型にして解く方法の一例として、非特許文献１に記載の有限体積法の場合に用いられている式を説明する。変数φに対する輸送方程式は、式（２－１）あるいは式（２－２）のように表される。上述したように、この式（２－２）の左辺第２項が、移流項である。なお、ナビエ・ストークス方程式は、変数φを速度ｕあるいは運動量ρｕとした輸送方程式と言うことが出来る。
Figure JPOXMLDOC01-appb-M000009
　密度ρが一定であると仮定すると、一般に連続体は、連続の式と呼ばれる式（２－３）を満たしている。
Figure JPOXMLDOC01-appb-M000010
　そこで、この式（２－３）の添え字ｉをｊに変えて、両辺にφを乗じた式を、式（２－２）に加えると、式（２－４）が得られる。さらに、式（２－４）の左辺第２項と第３項との和を変形すると式（２－５）が得られる。式（２－５）の左辺第２項が、密度が一定のときの保存型の移流項である。
Figure JPOXMLDOC01-appb-M000011
　一方、連続体の密度ρが一定でなく、ρ＝ρ（ｔ，ｘ，ｙ，ｚ）である場合は、まず、式（２－１）、式（２－２）の両辺にρを乗じて、式（２－６）、式（２－７）を得る。また、密度ρが一定でない場合の連続の式は、式（２－８）である。
Figure JPOXMLDOC01-appb-M000012
　式（２－４）のときと同様に、式（２－８）の添え字ｉをｊに変えて、両辺にφを乗じた式を、式（２－７）に加えると、式（２－９）が得られる。さらに、式（２－９）の左辺を変形すると式（２－１０）が得られる。式（２－１０）の左辺第２項が、密度が一定でない（密度が変化する）ときの保存型の移流項である。
Figure JPOXMLDOC01-appb-M000013
　式（２－１０）において、ρφ→φ’、ρＦ（φ）→Ｆ’（φ’）と置き換えると、式（２－１０）は、式（２－５）と同じ形で表すことができるので、以降では、式（２－５）の保存型の輸送方程式を基に説明する。
　式（２－５）を、コントロールボリュームの体積で積分すると、式（２－１１）が得られる。また、連続の式（２－３）についても、同様にコントロールボリュームの体積で積分すると、式（２－１２）が得られる。
Figure JPOXMLDOC01-appb-M000014
　積分方程式である式（２－１１）、式（２－１２）の各々を、ｍ個の境界面Ｅａｂ（ｂは１～ｍ）を持つコントロームボリュームａについて離散化すると、以下の式（２－１３）、式（２－１４）が得られる。なお、式（２－１３）の左辺第３項である[Ｆ（φ）の離散化項］は、Ｆ（φ）を体積積分して離散化した項であるが、詳細な説明は省略する。
Figure JPOXMLDOC01-appb-M000015
　本発明が解決しようとする課題を分かり易くするために、式（２－１３）において、Ｆ（φ）＝０とした式（２－１５）を用いて説明する。
Figure JPOXMLDOC01-appb-M000016
　上述したように、式（２－１５）は、変数φに対する保存型の輸送方程式を、体積積分し、離散化した式である。すなわち、物理量φが、連続体の移動、変形に伴って、移流する様子を表した式である。ある時点で、連続体の全ての領域でφが同じ値であれば、連続体が移動、変形していても、式（２－１５）でシミュレーションした結果は、連続体の全ての領域でφが変化せず、同じ値となる必要がある（条件Ａ）。例えば、物理量φが水の塩分濃度であり、全ての領域で塩分濃度φが同じであるときには、時間が経過して、どのように水が流れても、塩分濃度φは変化しない。この条件Ａについて、図２０に示す１次元流れ問題を例に考える。
　図２０に示す１次元流れでは、流体はｘ軸方向の正の向きに流れており、この流れにより物理量φが移流している。この１次元流れの流路を、図２０に示すように、ｘ軸に沿ってコントロールボリュームａ、ｂ、ｃに分割する。なお、コントロールボリュームｃは、コントロールボリュームａの上流側に隣接しており、コントロールボリュームｂは、コントロールボリュームａの下流側に隣接している。流路の断面積は一定であり、１とする。すなわち、式（２－１６）が成り立つ。
　Ｓａｂ＝Ｓａｃ＝１　・・・（２－１６）
また、コントロールボリューム間の境界面は、ｘ軸に垂直であるとすると、式（２－１７）が成り立つ。各コントロールボリュームの幅は同じで、Δｘであるとすると、式（２－１８）が成り立つ。また、コントロールボリュームａの境界面における速度ベクトルｕａｂ、ｕａｃは、式（２－１９）で表されるとする。
Figure JPOXMLDOC01-appb-M000017
　これらを用いて、式（２－１５）を表すと、式（２－２０）となる。同様に、これらを用いて連続の式（２－１４）を表すと、式（２－２１）となる。
Figure JPOXMLDOC01-appb-M000018
 したがって、ｕａｂ＝ｕａｃ＝Ｕ０として、式（２－２０）の両辺をΔｘで割ると、式（２－２２）が得られる。
Figure JPOXMLDOC01-appb-M000019
　ここで、流路全体で、φが一定である、すなわち、φａｂ＝φａｃとすると、式（２－２２）から、式（２－２３）が得られる。この式（２－２３）は、φａは、時間が進行しても増減せずに一定の値となることを示す。これは、式（２－１５）の説明で述べた条件Ａを満たしている。
Figure JPOXMLDOC01-appb-M000020
　ところが、実際のシミュレーション計算の結果は、数値計算上の誤差を含んでいる。また、測定値は測定誤差を含んでいる。そのため、数値計算により得られた速度であっても、測定により得られた速度であっても、離散化された連続の式（２－１４）を完全に満たしていない。すなわち、式（２－２１）は満たされず、ｕａｂ－ｕａｃ≠０である。
　そこで、その誤差をδＵ０とし、ｕａｂ－ｕａｃ＝δＵ０とする。これを、式（２－２０）に代入すると、φａｂ＝φａｃ＝φ＝一定としても、式（２－２３）とは異なる式（２－２４）が得られる。
Figure JPOXMLDOC01-appb-M000021
　この式（２－２４）は、δＵ0が正の値であるときは、コントロールボリュームａには、｜δＵ0・φ／Δｘ｜の大きさの擬似的なシンク（sink）があり、時間が進行するとφａは減少していくことを示す。また、δＵ0が負の値であるときは、コントロールボリュームａには、｜δＵ0・φ／Δｘ｜の大きさの擬似的なソース（source）があり、時間が進行するとφａは増加していくことを示す。
　すなわち、速度が連続の式を満たしていないときは、式（２－１５）を用いても条件Ａは満たされない。
　ここでは、１次元流れ問題を例に考えたが、２次元、３次元の輸送方程式でも同様である。つまり、連続の式（２－１４）が完全には満足されず、式（２－２５）であるときは、計算領域全てでφが一定であっても、式（２－２６）となる。φａは、時間が進行すると増減してしまい、条件Ａは満たされない。
Figure JPOXMLDOC01-appb-M000022
　以上を纏めると、以下の１）、２）が言える。
１)数値シミュレーションによって算出された連続体の速度は、完全に連続の式を満たすことは無い。有効桁数を大きくし、十分に高精度な計算を行っても、誤差は発生する。また、計測により得られた連続体の速度も、計測誤差を必ず含んでいるので、完全に連続の式を満たすことは無い。
２）移流項を保存型にした輸送方程式では、速度が連続の式を満足していないと、その誤差に比例した大きさの擬似的なソース／シンク項が発生してしまう。
　さらに、１）で述べたように、擬似的なソース／シンク項の大きさは、誤差に比例していることから以下の３）、４）が言える。
３）誤差が小さいときは、擬似的なソース／シンク項の大きさも小さくなる。この場合、対象の方程式に含まれている拡散項や、メッシュ分割による数値拡散があるため、擬似的なソース／シンク項により発生した増減は拡散され、目立たなくなる。したがって、良好な解が得られているように見える。
４）解析領域が複雑な形状の場合は、メッシュ分割をすると、歪んだ形状のコントロールボリュームが生成されることが多い。コントロールボリュームの形状が歪んでいると、数値シミュレーションによって速度を算出する際に、誤差が大きくなりやすい。このため、擬似的なソース／シンク項により発生した増減も大きくなり、数値拡散では抑えきれずに、解の精度を悪化させて異常な解を生じさせる。さらに、著しい場合には、計算を発散させてしまうこともある。
　このように、移流項を保存型にした輸送方程式を用いたシミュレーションでは、速度が連続の式を満足していないと、異常な解が生じてしまうことがあるという問題がある。
非特許文献２、３では、積分方程式が解析的には保存則を満足していても、これらを離散化した時点で保存則が満足されない問題点について提起している。しかし、非特許文献２、３では、離散化式において保存則を満たすように変形することによって、この保存則が満足されない問題についての解決方法を提示している。すなわち、非特許文献２、３は、数値計算上の誤差の問題については何ら解決策を提示していない。このため、非特許文献２、３の方法によって離散化段階で保存則を満足することによって精度を向上しても、離散化方程式の入力データとなる速度場が誤差を持っている場合には、依然として異常な解が生じてしまう問題がある。
　本発明は、このような事情に鑑みてなされたもので、その目的は、良好なシミュレーション結果を得ることができるシミュレーション装置、シミュレーション方法およびプログラムを提供することにある。
　上述したように、従来、輸送方程式を保存型に変形する際には、連続の式（２－３）、（２－８）が完全に満たされていることを前提にして、連続の式（２－３）、（２－８）にφを乗じたものを、輸送方程式に加算している。しかし、上述したように、数値シミュレーションによって算出したものでも、測定により得たものでも、速度は誤差を含んでおり、連続の式を完全には満足していない。すなわち、この前提は成り立っていないために、保存型にした移流項に誤差が生じてしまう。そこで、本発明では、この保存型にした移流項で生じる誤差を修正する項を、方程式に含めておき、該方程式の解を算出する。
　従来の研究は、離散化上の問題に着目して行われてきた。しかし、工業的に数値シミュレーションを利用する場合には、通常、莫大な容量の計算で、計算対象の形状が複雑で、複雑な現象の再現が求められる。このため、本発明者らは、数値シミュレーションを利用する場合に、単に離散化上の精度向上の対策では不十分であり数値計算上の誤差の存在を前提にした対策が重要であることに着目して、この発明をした。
（１）この発明は上述した課題を解決するためになされたもので、本発明の一態様は、移流項を含む方程式の解を算出するシミュレーション装置であって、速度場を取得する速度場取得部と、前記方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部とを具備し、前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであることを特徴とする。
（２）また、本発明の他の態様は、上述のシミュレーション装置であって、前記移流項の保存型の誤差を修正する項が、前記速度場が連続の式を満たさないことにより生じる誤差を修正する項である。
（３）また、本発明の他の態様は、上述のシミュレーション装置であって、前記誤差を修正する項が、前記連続の式に対して、前記方程式に対する積分と同じ積分を少なくとも行うことで得られる項である。
（４）また、本発明の他の態様は、上述のシミュレーション装置であって、前記解算出部は、少なくとも連続の式を用いて、次の時間ステップにおける速度場を算出し、前記速度場取得部が取得する速度場は、前記解算出部が算出した速度場である。
（５）また、本発明の他の態様は、上述のシミュレーション装置であって、前記速度場取得部は、予め記憶された速度場を取得する。
（６）また、本発明の他の態様は、移流項を含む方程式の解を算出するシミュレーション方法であって、速度場を取得する第１の過程と、移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記第１の過程で取得した速度場を適用して、前記離散化方程式における解を算出する第２の過程とを有し、前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであることを特徴とする。
（７）また、本発明の他の態様は、コンピュータを、速度場を取得する速度場取得部、移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部として機能させるためのプログラムであって、前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであることを特徴とするプログラムである。
（８）また、本発明の他の態様は、複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション装置であって、前記観測器の位置情報および計測した濃度情報を入手する観測情報入手手段と、拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定手段と、前記観測器の位置情報、計測した濃度情報、および前記仮想放出地点の情報に基づいて前記仮想放出地点での拡散物質の濃度の推定をする上述のシミュレーション装置と、を有することを特徴とする。
（９）また、本発明の他の態様は、複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション方法であって、前記観測器の位置情報および計測した濃度情報を入手する観測情報入手ステップと、拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定ステップと、前記観測器の位置情報、計測した濃度情報および前記仮想放出地点の情報に基づいて上述のシミュレーション方法によって前記仮想放出地点での拡散物質の濃度の推定をするステップと、を有することを特徴とする。
　この発明によれば、良好なシミュレーション結果を得ることができる。
簡単なシミュレーション例におけるコントロールボリュームを示す図である。 簡単なシミュレーション例における連続の式の値を示す図である。 従来の方法による簡単なシミュレーション例を示す図である。 本発明の簡単なシミュレーション例を示す図である。 この発明の第１の実施形態による流体シミュレーション装置１０の構成を示す概略ブロック図である。 同実施形態におけるセルを説明する概念図である。 同実施形態におけるセル構成記憶部１０２が記憶するセル情報の例を示す図である。 同実施形態におけるセル構成記憶部１０２が記憶するセル間情報の例を示す図である。 同実施形態における速度場記憶部１０４が記憶する速度場情報の例を示す図である。 同実施形態における圧力場記憶部１０７が記憶する圧力場情報の例を示す図である。 同実施形態における流体シミュレーション装置１０の動作を説明するフローチャートである。 この発明の第２の実施形態による流体シミュレーション装置１０ａの構成を示す概略ブロック図である。 この発明の第３の実施形態による連続体シミュレーション装置１０ｂの構成を示す概略ブロック図である。 この発明の第４の実施形態による連続体シミュレーション装置１０ｃの構成を示す概略ブロック図である。 シミュレーション対象の速度場を示すコンター図である。 シミュレーション対象の速度場における速度の発散の大きさを示すコンター図である。 温度分布を、従来の方法を用いて算出した場合の例を示すコンター図である。 温度分布を、本発明を用いて算出した場合の例を示すコンター図である。 オイラー座標系とラグランジュ座標系とを説明する図である。 一次元流れにおける各パラメータを説明する図である。
　各実施形態を説明する前に、本発明の概略を説明する。まず、連続体の密度ρが一定の場合について説明する。上述したように、従来、輸送方程式（２－１）、（２－２）を、保存型の式（２－５）に変形する際には、式（２－３）でも示した連続の式（３－１）を用いている。
Figure JPOXMLDOC01-appb-M000023
　しかし、上述したように、シミュレーションや測定により得られた速度は、誤差を含んでいるために、式（３－１）を完全には満たさず、右辺はゼロにならない。そこで、ここでは、式（３－１）の右辺をゼロからＤに置き換えた式（３－２）を用いて、保存型の輸送方程式を導出する。すなわち、連続の式が完全には満たされず、式（３－１）の右辺がゼロにはならないことを前提として、保存型の輸送方程式を導出する。
Figure JPOXMLDOC01-appb-M000024
　まず、輸送方程式（２－２）の左辺第２項と第３項の間に、+φ・Ｄ－φ・Ｄを挿入すると、式（３－３）が得られる。次に、式（３－３）の左辺第３項のＤに、式（３－２）を代入すると、式（３－４）が得られる。この式（３－４）の左辺第２項と第３項とは、式（２－４）と同じであるので、同様に変形すると、式（３－５）が得られる。
Figure JPOXMLDOC01-appb-M000025
　本発明では、密度ρが一定のときは、この式（３－５）を保存型の輸送方程式として用いる。すなわち、本発明では、式（３－５）の解を算出する。なお、この（３－５）は、左辺第３項（－φ・Ｄ）を有することを除けば、従来の保存型の輸送方程式（２－５）と同一である。この左辺第３項を有することにより、式（３－２）に示したように、連続の式の右辺がゼロとならず、誤差Ｄを含むことを前提とした、保存型の輸送方程式となっている。この式（３－５）を用いたシミュレーションでは、式（２－２４）に現れたような擬似的なソース／シンクを、－φ・Ｄにより修正することができる。すなわち、擬似的なソース／シンクが現れないようにすることができる。
　式（３－５）を用いたシミュレーションの例として、式（３－５）を有限体積法で離散化して、シミュレーションする場合を説明する。なお、以下では、有限体積法で離散化する場合を説明しているが、差分法など、保存型の移流項を持つ方程式を離散化する方法であれば、その他の離散化の方法であってもよい。
　まず、式（２－１１）を導出したときと同様に、式（３－５）を、コントロールボリュームの体積Ｖに対して積分すると、式（３－６）が得られる。式（３－６）における各記号は、式（２－１１）と同様である。また、式（３－２）も、同様にコントロールボリュームの体積Ｖに対して積分すると、式（３－７）が得られる。
Figure JPOXMLDOC01-appb-M000026
　これらの式（３－６）、（３－７）を、式（２－１３）、（２－１４）を導出した際と同様に、ｍ個の境界面Ｅａｂ（ｂは１～ｍ）を持つコントロームボリュームａについて離散化すると、以下の式（３－８）、式（３－９）が得られる。
Figure JPOXMLDOC01-appb-M000027
さらに、式（３－８）に、式（３－９）を代入すると、式（３－１０）が得られる。この離散化した輸送方程式（３－１０）を解くことで、良好なシミュレーション結果を得ることが出来る。
Figure JPOXMLDOC01-appb-M000028
　次に、連続体の密度ρが一定でない場合を説明する。密度ρが一定でない場合は、密度ρが一定のときと同様に、密度ρが一定でないときの連続の式の右辺をＤとした式（３－１１）を用いて、保存型の輸送方程式を導出する。
Figure JPOXMLDOC01-appb-M000029
　密度ρが一定でないときの輸送方程式（２－７）の左辺第１項と第２項の間に、+φ・Ｄ－φ・Ｄを挿入すると、式（３－１２）が得られる。次に、式（３－１２）の左辺第３項のＤに、式（３－１１）を代入すると、式（３－１３）が得られる。この式（３－１３）の左辺第１項、第２項は、式（２－９）と同様であるので、同様に変形すると、式（３－１４）が得られる。
Figure JPOXMLDOC01-appb-M000030
　本発明では、密度ρが一定でないときは、この式（３－１４）を保存型の輸送方程式として用いる。すなわち、式（３－１４）の解を算出する。なお、この（３－１４）は、左辺第３項（－φ・Ｄ）を有することを除けば、従来の保存型の輸送方程式（２－１０）と同一である。この左辺第３項を有することにより、式（３－１１）に示したように、連続の式の右辺がゼロとならず、誤差Ｄを含むことを前提とした、保存型の輸送方程式となっている。この式（３－１４）を用いたシミュレーションでは、式（２－２４）に現れたような擬似的なソース／シンクを、－φ・Ｄにより修正することができる。すなわち、密度ρが一定でないときでも、擬似的なソース／シンクが現れないようにすることができる。
　次に、式（３－１４）を用いたシミュレーションの例として、密度ρが一定のときと同様に、式（３－１４）を有限体積法で離散化して、シミュレーションする場合を説明する。式（２－１１）を導出したときと同様に、式（３－１４）を、コントロールボリュームの体積Ｖに対して積分すると、式（３－１５）が得られる。式（３－１５）における各記号は、式（２－１１）と同様である。また、式（３－１１）も、同様にコントロールボリュームの体積Ｖに対して積分すると、式（３－１６）が得られる。
Figure JPOXMLDOC01-appb-M000031
　これらの式（３－１５）、（３－１６）を、式（３－８）、（３－９）を導出した際と同様に、ｍ個の境界面Ｅａｂ（ｂは１～ｍ）を持つコントロームボリュームａについて離散化すると、以下の式（３－１７）、式（３－１８）が得られる。
Figure JPOXMLDOC01-appb-M000032
さらに、式（３－１７）に、式（３－１８）を代入すると、式（３－１９）が得られる。
Figure JPOXMLDOC01-appb-M000033
　ここで、式（３－２０）から、式（３－１９）は、式（３－２１）となる。この離散化した輸送方程式（３－２１）を解くことで、密度ρが一定でない場合でも、良好なシミュレーション結果を得ることが出来る。
Figure JPOXMLDOC01-appb-M000034
　ところで、式（３－１０）、（３－２１）をそれぞれ変形すると、以下の式（３－２２）、（３－２３）が得られる。
Figure JPOXMLDOC01-appb-M000035
　式（３－２２）、（３－２３）の各々の左辺第２項のΣの中にある、（φａｂ－φａ）に注目する。このφａｂは、コントロールボリュームａ、ｂの境界面Ｅａｂにおけるφの値である。φａｂは、φａとφｂとを用いた補間により算出する。すなわち、φａｂは、式（３－２４）で表すことができる。ここで、ｘａ、ｘｂは、コントロールボリュームａ、ｂ各々の代表点であるコントロールポイントＣa、Ｃｂの座標ベクトルである。
Figure JPOXMLDOC01-appb-M000036
　このφａｂを算出するための補間関数は、様々なものが従来から知られている。いずれの補間関数（スキームともいう）を用いるかは、シミュレーションの目的、対象とする連続体の種類、数値計算方法などにより決められるが、いずれを用いるようにしてもよい。ここでは、代表的なスキームをいくつか示す。
（１）一次風上スキーム
　流体を対象とするときに良く使用されるスキームである。φａｂを導出するための補間関数として、以下の式（３－２５）、（３－２６）を用いる。単位法線ベクトルｎａｂは、コントロールボリュームａに対して外向きのベクトルであるので、式（３－２５）は、コントロールボリュームａに対して、流入の場合の補間関数であり、式（３－２６）は、流出の場合の補間関数である。
Figure JPOXMLDOC01-appb-M000037
　これらを、式（３－２２）、（３－２３）の（φａｂ－φａ）に代入すると、以下の式（３－２７）、（３－２８）となる。
Figure JPOXMLDOC01-appb-M000038
　従って、式（３－２２）、（３－２３）のシグマのうち、流出する境界面の項目は０になる。そこで、流入する境界面についての総和を式（３－２９）のように記載すると、式（３－２２）、（３－２３）は、式（３－３０）、（３－３１）となる。
Figure JPOXMLDOC01-appb-M000039
　これら、式（３－３０）、（３－３１）におけるチルダ付きのシグマは、コントロールボリュームａと隣接するコントロールボリュームｂからφｂが流入する際の連続体の流量Ｓａｂ・（ｎａｂ・ｕａｂ）と同じ流量でφａが流出することを示す式となっている。すなわち、コントロールボリュームａに対する流入/流出の流量の収支が合う式となっている。本発明の有効性を示す一例となっている。
　一次風上スキームを用いた場合の簡単なシミュレーション例を示す。１次元流れで、密度は一定であるとする。図１に示すように、流れ方向に４つのコントロールボリュームａ１、ａ２、ａ３、ａ４に分割する。各境界面における速度は、コントロールボリュームａ１の上流側の境界面から順に、１．０、１．１、０．９、０．８、１．０であったとする。各コントロールボリュームの体積を１、境界面を１とし、この速度分布を各コントロールボリュームにおいて連続の式に代入すると、図２のようにコントロールボリュームａ１から順に、－０．１、０．２、０．１、－０．２となり、いずれも０とならない。
　式（３－３０）は、簡単のためにＦ（φ）＝０、時間微分を１次差分で近似し、Δｔ＝１とすると、式（３－３０’）となる。
　φｉ＋ｕｉ（φｉ－１－φｉ）＝
　（１－ｕｉ）φａ＋ｕｉφｉ－１＝φｉ’　・・・（３－３０’）
　ここで、ｕｉは、コントロールボリュームａｉの上流側のコントロールボリュームａｉ－１との境界面での速度、φｉは、当該コントロールボリュームａｉの物理量（温度）、φｉ’は、前の時間ステップにおける当該コントロールボリュームの物理量、φｉ―１は、上流側のコントロールボリュームａｉ－１の物理量である。
　一方、従来は、式（２－２０）に同様の条件を当てはめて、式（２－２０’）となる。
　φｉ＋ｕｉ・φｉ－１－ｕｉ＋１・φｉ＝
　（１－ｕｉ＋１）φｉ＋ｕｉ・φｉ－１＝φｉ’　・・・（２－２０’）
　ここで、初期条件として全てのコントロールボリュームで温度が１００℃、すなわちφｉ’＝１００℃（ｉ＝１～４）として、式（２－２０’）、（３－３０’）を用いて、次の時間ステップにおける各コントロールボリュームの温度を算出する。式（２－２０’）を用いたとき、すなわち従来の算出結果例を図３に、式（３－３０’）を用いたとき、すなわち本発明の算出結果例を図４に示す。このように、一様に１００℃であったので、本来は算出結果も一様に１００℃となるべきであるが、従来は、一様な分布にならないが、本発明の算出結果例では、一様に１００℃となる。
（２）二次風上スキーム、三次風上スキーム
　二次風上スキーム、三次風上スキームとも、流体計算で良く使用されるスキームである。これらのスキームの詳細は、ここでは省略するが、φａｂは、式（３－３２）となる。すなわち、φａ、φｂのｇｒａｄ（勾配ベクトル）を使用した補間関数を用いる。
Figure JPOXMLDOC01-appb-M000040
（１）の風上スキームと異なり、式（３－２８）のように、流出する場合に移流項が単純にゼロにはならない。したがって、式（３－２２）、（３－２３）におけるシグマの項は、コントロールボリュームａに隣接する全てのコントロールボリュームｂに関する移流項の総和となる。
（３）その他のスキーム
　一次、二次、三次風上スキーム以外にも、指数関数ｅｘなど、非線形の関数を利用したスキームがある。また、固体や流れの遅い流体などでは、中心差分と同じ考え方から、式（３－３３）が用いられることがある。
Figure JPOXMLDOC01-appb-M000041
　以上、いくつかのφａｂを導出するためのスキームの例を示したが、どのようなスキームを用いても、連続体の移動、変形速度が連続の式を満たさず誤差を含んだ状態であっても、良好なシミュレーション結果を得ることができる。
［第１の実施形態］
　以下、図面を参照して、本発明の第１の実施形態について説明する。図５は、この発明の第１の実施形態による流体シミュレーション装置１０の構成を示す概略ブロック図である。本実施形態における流体シミュレーション装置１０は、有限体積法により、非圧縮性流体をシミュレーションして、シミュレーションする領域における速度と圧力とを算出するシミュレーション装置である。すなわち、本実施形態では、運動量に関する輸送方程式であるナビエ・ストークス方程式の解を算出するシミュレーション装置を説明する。
　図５に示すように流体シミュレーション装置１０は、セル構成設定部１０１、セル構成記憶部１０２、初期条件設定部１０３、速度場記憶部１０４、速度場取得部１０５、速度・圧力算出部１０６、圧力場記憶部１０７を含んで構成される。セル構成設定部１０１は、シミュレーションする領域（空間）を分割したセルに関する情報を設定し、セル構成記憶部１０２に記憶させる。セル構成記憶部１０２は、セル構成設定部１０１が設定したセルに関する情報を記憶する。なお、セルに関する情報は、各セルの情報であるセル情報と、隣接するセル間の情報であるセル間情報とを含む。
　初期条件設定部１０３は、初期条件として、各セルの速度ベクトルを設定し、速度場記憶部１０４に記憶させる。速度場記憶部１０４は、初期条件設定部１０３が設定した各セルの速度ベクトルを記憶する。また、速度場記憶部１０４は、各時間ステップについて、速度・圧力算出部１０６が算出した各セルの速度ベクトルも記憶する。速度場取得部１０５は、速度・圧力算出部１０６が時間ステップｎに関する解を算出するときは、時間ステップｎ－１における各セルの速度ベクトルを、速度場記憶部１０４から読み出す。
　速度・圧力算出部１０６は、速度場取得部１０５が速度場記憶部１０４から読み出した時間ステップｎ－１における各セルの速度ベクトルと、セル構成記憶部１０２が記憶するセル情報とを用いて、対象時間ステップｎにおける各セルの速度ベクトルと圧力とを算出する。速度・圧力算出部１０６は、算出した各セルの速度ベクトルを、速度場記憶部１０４に記憶させ、圧力を圧力場記憶部１０７に記憶させる。なお、速度・圧力算出部１０６は、速度ベクトルと圧力とを算出する際に、非圧縮性流体のナビエ・ストークス方程式と、非圧縮性流体の連続の式（質量保存の式ともいう）を離散化した離散化方程式を用いる。本実施形態における非圧縮性流体のナビエ・ストークス方程式の離散化方程式は、移流項を後述する保存型に変形することと、移流項の保存型の誤差を修正する項を含める操作を非圧縮性流体のナビエ・ストークス方程式に対して行っている。これらの詳細については後述する。圧力場記憶部１０７は、各時間ステップについて、速度・圧力算出部１０６が算出した各セルの圧力を記憶する。
　図６は、本実施形態におけるセルを説明する概念図である。図６では、符号Ｒ１～Ｒ１５を付した領域の各々が１～１５番目のセルである。図６に示すように、各セルは、多面体である。また、シミュレーションする領域を隙間なく埋めるように、セルは配置されている。同図において、代表点Ｃｎ（ｎ＝１～４）は、それぞれｎ番目のセルＲｎの代表点である。また、境界面Ｅ１２は、セルＲ１とＲ２との境界面である。単位法線ベクトルｎ１２は、境界面Ｅ１２の単位法線ベクトルである。距離α１２は、代表点Ｃ１からＣ２までの距離を１としたときの、点Ｃ１とＣ２を結ぶ線の境界面Ｅ１２との交点から代表点Ｃ１までの距離である。
　図７は、セル構成記憶部１０２が記憶するセル情報の例を示す図である。図７に示すようにセル情報は、全てのセルのセルＩＤ、代表点座標、体積を含む。図７に示す例では、セルＩＤが「０００１」のセルのセル情報は、セルＩＤ「０００１」、代表点座標「１５，２２５，２１４」、体積「２４．３」である。同様に、セルＩＤが「０００２」のセルのセル情報は、セルＩＤ「０００２」、代表点座標「１５，２２５，２１４」、体積「２４．３」であり、セルＩＤが「２４５３」のセルのセル情報は、セルＩＤ「２４５３」、代表点座標「２１５、２５、４４」、体積「３１．２」である。なお、本実施形態において、代表点座標の３つの数値は、ｘ座標と、ｙ座標と、ｚ座標である。以降、Ｘｎは、セルＩＤがｎのセルの代表点座標を示し、Ｖｎは、セルＩＤがｎのセルの体積を示す。
　図８は、セル構成記憶部１０２が記憶するセル間情報の例を示す図である。図８に示すようにセル間情報は、基準とするセルのセルＩＤ、該セルに隣接するセルのセルＩＤ、これら２つのセル間の境界面の面積、該境界面の単位法線ベクトル、これら２つのセルの代表点を通る直線と境界面との交点から基準とするセルの代表点までの距離αを含む。なお、単位法線ベクトルは、長さが「１」の単位ベクトルである。また、該距離は、２つのセルの代表点間の距離を「１」としたときの距離であり、図６のα１２に対応する。
　図８に示す例では、セルＩＤ「０００１」のセルとセルＩＤ「０００５」のセルとのセル間に関するセル間情報は、セルＩＤ「０００１」、セルＩＤ「０００５」、面積「４．２１４５８」、単位法線ベクトル「０．７００２４３，０．５５２４６６，０．８９１９４１」、距離「０．４２２５６１」である。同様に、セルＩＤ「０００１」のセルとセルＩＤ「０１０２」のセルとのセル間に関するセル間情報は、セルＩＤ「０００１」、セルＩＤ「０１０２」、面積「３．１３３５２」、単位法線ベクトル「０．５８４２７６，－０．６１７１２３，０．５２７０４９」、距離「０．５５４８７８」であり、セルＩＤ「２４５３」のセルとセルＩＤ「２４５２」のセルとのセル間に関するセル間情報は、セルＩＤ「２４５３」、セルＩＤ「２４５２」、面積「３．３１２１５」、単位法線ベクトル「０．３３５７８６，０．８１４２６６，０．４７３５１７」、距離「０．５１５４８６」である。
　図９は、速度場記憶部１０４が記憶する速度場情報の例を示す図である。速度場記憶部１０４は、初期条件として設定された初期時間ステップにおける速度場情報と、速度・圧力算出部１０６が算出した各時間ステップにおける速度場情報とを記憶する。図９に示す例は、それらのうちの１つの時間ステップにおける速度場情報である。１つの時間ステップにおける速度場情報は、各セルのセルＩＤと該セルの速度ベクトルとを含む。図９に示す速度場情報の例は、セルＩＤ「０００１」と該セルの速度ベクトル「１．２２１３５４，３．２５５４６，０．２５１３５６」と、セルＩＤ「０００２」と該セルの速度ベクトル「１．２４０１３４，３．１５２２４，０．６６２１５８」と、セルＩＤ「２４５３」と該セルの速度ベクトル「３．２２１５１，３．００１２５，１．６５８５５」とを含む。
　図１０は、圧力場記憶部１０７が記憶する圧力場情報の例を示す図である。圧力場記憶部１０７は、速度・圧力算出部１０６が算出した各時間ステップにおける圧力場情報を記憶する。図１０に示す例は、それらのうちの１つの時間ステップにおける圧力場情報である。１つの時間ステップにおける圧力場情報は、各セルのセルＩＤと該セルの圧力とを含む。図１０に示す圧力場情報の例は、セルＩＤ「０００１」と該セルの圧力「１．２８５６８」と、セルＩＤ「０００２」と該セルの圧力「１．２２３２５」と、セルＩＤ「２４５３」と該セルの圧力「１．５６８４６」とを含む。
　図１１は、本実施形態における流体シミュレーション装置１０の動作を説明するフローチャートである。まず、セル構成設定部１０１が、外部から取得したセルに関する情報（セル情報およびセル間情報）をセル構成記憶部１０２に記憶させる。また、初期条件設定部１０３が、外部から取得した初期条件である初期時間ステップの速度場情報を、速度場記憶部１０４に記憶させる（Ｓ１）。次に、速度場取得部１０５は、対象時間ステップを初期値「１」とする（Ｓ２）。次に、速度場取得部１０５が、速度場記憶部１０４から、対象時間ステップの１ステップ前の速度場情報を取得する（Ｓ３）。なお、速度場取得部１０５は、対象時間ステップが「１」のときは、初期時間ステップを１ステップ前とし、初期時間ステップの速度場情報を、速度場記憶部１０４から取得する。
　次に、速度・圧力算出部１０６は、非圧縮性流体のナビエ・ストークス方程式を離散化した式と、非圧縮性流体の連続の式を離散化した式とからなる連立一次方程式を、Ａ・φ＝ｂとしたときの行列Ａとベクトルｂとを算出する（Ｓ４）。速度・圧力算出部１０６は、ステップＳ３にて速度場取得部１０５が取得した速度場情報を用いて、行列Ａおよびベクトルｂの各要素を算出する。
　なお、ベクトルφは、対象時間ステップにおける全てのセルの速度ベクトルの各要素と圧力とを要素とするベクトルである。すなわち、ベクトルφは、セルの数×４個の要素を持つ。一方、前記連立一次方程式は、各セルに関する、各軸方向のナビエ・ストークス方程式を離散化した式と、各セルに関する非圧縮性流体の連続の式を離散化した式とからなる。すなわち、前記連立一次方程式は、ベクトルφの要素数と同じ数の式からなるため、解の算出が可能である。なお、ステップＳ４の行列Ａと、ベクトルｂとの算出方法についての詳細は、後述する。
　次に、速度・圧力算出部１０６は、前記Ａ・φ＝ｂの解であるベクトルφを算出する。前記Ａ・φ＝ｂは、ベクトルφの各要素を未知数とする連立一次方程式であるので、速度・圧力算出部１０６は、ベクトルφの算出には、例えば、共役勾配法などの公知の連立一次方程式の解法を用いる。速度・圧力算出部１０６は、算出したベクトルφの要素のうちの速度ベクトルを、対象時間ステップの速度ベクトルとして速度場記憶部１０４に記憶させる。同様に、速度・圧力算出部１０６は、算出したベクトルφの要素のうちの圧力を、対象時間ステップの圧力として圧力場記憶部１０７に記憶させる（Ｓ５）。
　次に、速度・圧力算出部１０６は、終了条件が満たされているか否かを判定する（Ｓ６）。この終了条件としては、例えば、対象時間ステップが、予め設定された値と一致することなどが挙げられる。終了条件が満たされていないと判定したときは（Ｓ６－Ｎｏ）、ステップＳ７に遷移して、速度・圧力算出部１０６は、対象時間ステップを「１」増加させ（Ｓ７）、ステップＳ３に戻る。これにより、次の時間ステップについての処理を行う。また、ステップＳ６にて、終了条件が満たされていると判定したときは（Ｓ６－Ｙｅｓ）、処理を終了する。
　次に、速度・圧力算出部１０６による行列Ａ、ベクトルｂの算出処理、すなわち図１１のステップＳ４の詳細を説明する。前述のように、Ａ・φ＝ｂは、非圧縮性流体のナビエ・ストークス方程式を離散化した式と、非圧縮性流体の連続の式を離散化した式との連立一次方程式である。一般に、第ｉ軸方向の非圧縮性流体のナビエ・ストークス方程式は、式（４－１）のように表される。また、一般に、非圧縮性流体の連続の式は、式（４－２）のように表される。なお、本実施形態では、座標系として直交座標系（カーテシアン座標系）を用いる。また、式（４－１）（４－２）を含む以降の各式では、アインシュタインの総和則を用いる。
Figure JPOXMLDOC01-appb-M000042
　式（４－１）において、ｔは時間である。ｕｉ（ｉ＝１，２，３）は、速度の各軸方向の成分である。ｘｉ（ｉ＝１，２，３）は、座標の各軸成分である。pは圧力である。ρは密度である。νは、動粘性係数である。Ｆｉ（ｉ＝１，２，３）は、外力の各軸方向の成分である。式（４－１）では、左辺の第２項が移流項である。なお、式（４－１）は第ｉ軸方向の運動量保存の式であることと、非圧縮性流体であるため流体の密度ρは一定であることとから、式（４－１）における移流項では、移流されるのは第ｉ軸方向の速度ｕｉとなっている。以降では、説明の簡易のために、式（４－１）の右辺の第２項（粘性項）と第３項（外力）との和を－Ｆ（ｕ）と表記し、右辺の全ての項を左辺に移動させた式（４－３）を用いて説明する。
Figure JPOXMLDOC01-appb-M000043
　従来の有限体積法では、ナビエ・ストークス方程式をコントロールボリュームで体積積分し、さらに離散化した式を解いている。有限体積法では、この体積積分の際に、連続の式である式（４－２）の両辺に速度ｕｉを乗じた式（４－４）をナビエ・ストークス方程式に加えて、式（４－５）のように移流項を保存型と呼ばれる形に変形したものを用いる。この移流項の保存型への変形は、速度ベクトルが連続の式（４－２）を満たしていること前提とした変形である。
Figure JPOXMLDOC01-appb-M000044
　一方、本実施形態では、速度ベクトルｕを、連続の式の左辺に代入しても「０」にはならず、式（４－６）に示すように誤差Ｄが存在すると仮定する。例えば、計算容量などの問題から、反復計算により速度を求める際に充分に収束するまで計算できないときなどに、誤差Ｄの大きさが無視できないものとなることがある。このような誤差は、セルの形状が、特定の軸方向に長いときなどに発生しやすい。
　次に、式（４－６）の両辺に速度ｕｉを乗じた式（４－７）をナビエ・ストークス方程式（４－１）に加えると、式（４－８）が得られる。この式（４－８）の移流項を保存型に変形し、さらに式（４－８）の右辺を左辺に移動させると、式（４－９）が得られる。この式（４－９）は、左辺の第３項があることを除けば、従来の有限体積法で用いている式（４－５）と同じであり、第３項に、保存型の移流項における誤差を修正する項（すなわち、－ｕｉＤ）がある。以降、該誤差を修正する項を含んだままで、体積積分および離散化を行う。
Figure JPOXMLDOC01-appb-M000045
　この式（４－９）を、コントロールボリュームで体積積分すると、式（４－１０）が得られる。また、式（４－６）についても、コントロールボリュームで体積積分すると、式（４－１１）が得られる。なお、式（４－１０）、（４－１１）において、Ｖはコントロールボリュームの体積であり、Ｓはコントロールボリュームの表面である。また、式（４－１０）の左辺第２項および式（４－１１）の左辺第１項における太字体のｎとｕは、それぞれコントロールボリュームの表面Ｓの単位法線ベクトルと、速度ベクトルである。単位法線ベクトルｎは、コントロールボリュームの表面Ｓに対して垂直で、外向き、かつ、長さが「１」のベクトルである。
Figure JPOXMLDOC01-appb-M000046
　これらの式（４－１０）、（４－１１）を、例えば、図６のセルＲ１のコントロールボリュームＶ１、境界面Ｅ１ｋ（ｋ∈Ｎ１であり、集合Ｎ１は{２、３、４、１１、１２、１３、１４}である）の面積Ｓ１ｋ、時間Δｔについて離散化し、代数方程式による近似式に変換すると、式（４－１２）、式（４－１３）が得られる。なお、境界面Ｅ１ｋは、セルＲ１とセルＲｋとの間の境界面である。集合Ｎ１は、セルＲ１に隣接するセルのセルＩＤの集合である。また、Δｔは、時間ステップの幅である。
Figure JPOXMLDOC01-appb-M000047
 面積Ｓ１ｋは、セルＲ１とセルＲｋとの間の境界面の面積である。また、Ｄ１は、セルＲ１における連続の式の誤差である。ｕ１ｉは、セルＲ１における速度のｉ番目の軸方向の成分である。ｕ１ｉ（ｔ－Δｔ）は、対象時間ステップの前の時間ステップでのセルＲ１における速度のｉ番目の軸方向の成分である。ｕ１ｋｉは、境界面Ｅ１ｋにおける速度のｉ番目の軸方向の成分である。ｎ１ｋは、境界面Ｅ１ｋにおけるセルＲｋに向いた単位法線ベクトルである。u１ｋは、境界面Ｅ１ｋにおける速度ベクトルである。ｐｋは、セルＲｋにおける圧力である。
　次に、式（４－１３）を式（４－１２）の第３項に代入し、Σの項（第２項と第３項）を纏めると、式（４－１４）が得られる。
Figure JPOXMLDOC01-appb-M000048
　本実施形態では、この式（４－１４）を、前記連立一次方程式を構成する一次方程式の一つとしている。ただし、式（４－１４）は、第２項に未知数である速度同士の乗算を含んでおり、一次方程式となっていない。そこで、第２項における速度ベクトルｕ１ｋとして、前の時間ステップの速度ベクトルを用いるという近似を行う。また、Ｆ（ｕ）の離散化項については、説明を簡単にするため、前の時間ステップの速度ベクトルを用いるという近似を行う。
　また、境界面Ｅ１ｋにおける速度ｕ１ｋｉを、セルＲ１における速度ｕ１ｉと、セルＲｋにおける速度ｕｋｉとから補間する。補間方法としては、様々な方法があり、いずれの方法を用いてもよいが、本実施形態では一次風上スキームを用いる。一次風上スキームでは、境界面の速度として、風上側のセルにおける速度を用いる。単位法線ベクトルｎ１ｋは、セルＲ１から見て外向きのベクトルである。このため、単位法線ベクトルｎ１ｋと境界面の速度ベクトルｖ１ｋとの内積が負のときは、境界面Ｅ１ｋにおいてはセルＲｋが風上であり、正のときはセルＲ１が風上である。この一次風上スキームを式に表すと、式（４－１５）となる。なお、前記内積が負のときをセルＲ１へ流入、正のときをセルＲ１から流出ともいう。
Figure JPOXMLDOC01-appb-M000049
　この式（４－１５）を式（４－１４）に代入すると、セルＲ１から流出のときはΣの中は「０」になる。さらに、一次方程式とするために、境界面Ｅ１ｋの速度ベクトルｕ１ｋを、前の時間ステップにおける各セルの速度ベクトルから補間して求めた速度ベクトルｕ１ｋ’とする。本実施形態では、補間方法としては、式（４－１６）に示すように、セル構成記憶部１０２が記憶する距離α１ｋを用いて、前の時間ステップにおける速度ベクトルｕ１’とｕ１’とを比例配分する方法を用いる。
Figure JPOXMLDOC01-appb-M000050
　これらにより、式（４－１４）は、式（４－１７）のようになる。ただし、集合Ｎ１’は、セルＲ１に隣接するセルのうち、セルＲ１へ流入しているセルのセル番号を要素とする集合である。この式（４－１７）を変形すると式（４－１８）となり、セルＲ１の速度のｉ軸成分ｕ１ｉと、隣接する各セルＲｋの速度のｉ軸成分ｕｋｉと、セルＲ１の圧力ｐとを未知数とする一次方程式となる。なお、３次元であれば、ｉ＝１，２，３の各々について、式（４－１８）が立てられる。
Figure JPOXMLDOC01-appb-M000051
　また、連続の式である式（４－２）も、セルＲ１のコントロールボリュームＶ１、境界面Ｅ１ｋ（ｋ∈Ｎ１であり、集合Ｎ１は{２、３、４、１１、１２、１３、１４}である）の面積Ｓ１ｋについて離散化し、代数方程式による近似式に変換すると、式（４－１９）が得られる。前記式（４－１６）と同様の補間を用いると、式（４－１９）は式（４－２０）となり、セルＲ１の速度の各成分と、隣接セルＲｋの速度の各成分を未知数とする一次方程式となる。
Figure JPOXMLDOC01-appb-M000052
　これらの式（４－１８）、（４－２０）は、セルＲ１に関する式であるが、その他のセル各々についても同様の式を立てることができる。速度・圧力算出部１０６は、全てのセルについての式（４－１８）、（４－２０）における各未知数の係数を、行列Ａの要素とし、右辺をベクトルｂの要素とする。
　上述のように、非圧縮性流体のナビエ・ストークス方程式を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、本実施形態では、連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、各時間ステップの速度が、連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
　なお、本実施形態では、ナビエ・ストークス方程式と、連続の式とを連立させて解く場合を説明したが、これらに加えて他の方程式を連立させて解くようにしてもよい。その場合、他の方程式も移流項を含むときは、ナビエ・ストークス方程式と同様に移流項の誤差Ｄを修正する項を含めて離散化するようにしてもよい。これにより、該方程式についても誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
　また、速度・圧力算出部１０６が各時間ステップの解を算出する際に、共役勾配法などの反復法を用いる場合、解の初期値が必要になるが、初期値として前の時間ステップの解を用いるようにしてもよい。その場合、速度場取得部１０５は、圧力場記憶部１０７から前の時間ステップの圧力を取得し、速度・圧力算出部１０６に入力する。
［第２の実施形態］
　以下、図面を参照して、本発明の第２の実施形態について説明する。第２の実施形態では、圧縮性流体をシミュレーションする場合を説明する。図１２は、この発明の第２の実施形態による流体シミュレーション装置１０ａの構成を示す概略ブロック図である。本実施形態における流体シミュレーション装置１０ａは、有限体積法により、圧縮性流体をシミュレーションして、シミュレーションする領域における速度と圧力と密度とを算出するシミュレーション装置である。
　図１２に示すように流体シミュレーション装置１０ａは、セル構成設定部１０１、セル構成記憶部１０２、初期条件設定部１０３ａ、密度・速度場記憶部１０４ａ、密度・速度場取得部１０５ａ、解算出部１０６ａ、圧力場記憶部１０７を含んで構成される。同図において、図５の各部に対応する部分には、同一の符号（１０１、１０２、１０７）を付し、説明を省略する。
　初期条件設定部１０３ａ、初期条件として、各セルの速度ベクトルと密度とを設定し、密度・速度場記憶部１０４ａに記憶させる。密度・速度場記憶部１０４ａは、初期条件設定部１０３ａが設定した各セルの速度ベクトルと密度とを記憶する。また、密度・速度場記憶部１０４ａは、各時間ステップについて、解算出部１０６ａが算出した各セルの速度ベクトルと密度も記憶する。密度・速度場取得部１０５ａは、解算出部１０６ａが時間ステップｎに関する解を算出するときは、時間ステップｎ－１における各セルの速度ベクトルと密度を、密度・速度場記憶部１０４ａから読み出す。
　解算出部１０６ａは、密度・速度場取得部１０５ａが密度・速度場記憶部１０４ａから読み出した時間ステップｎ－１における各セルの速度ベクトルおよび密度と、セル構成記憶部１０２が記憶するセル情報とを用いて、対象時間ステップｎにおける各セルの速度ベクトルと圧力と密度とを算出する。解算出部１０６ａは、算出した各セルの速度ベクトルと密度とを、密度・速度場記憶部１０４に記憶させ、圧力を圧力場記憶部１０７に記憶させる。
　解算出部１０６ａは、圧縮性流体のナビエ・ストークス方程式を離散化した離散化方程式と、圧縮性流体の連続の式を離散化した離散化方程式と、密度と圧力の関係を表す状態方程式を離散化した離散化方程式とを含む連立一次方程式を解き、各セルの速度ベクトル、圧力、密度を算出する。解算出部１０６ａは、連立一次方程式を解く際に、例えば、共役勾配法などの公知の連立一次方程式の解法を用いる。なお、状態方程式としては、例えば、ｐ・１／ρ＝ＲＴを用いる。ここで、Ｒは気体定数であり、Ｔは温度（熱力学温度）である。本実施形態では、温度Ｔは一定である。なお、温度Ｔは、一定でなく、既知であってもよい。また、温度Ｔが、既知でないときは、温度Ｔを変数とする離散化方定式を、上述の連立一次方程式に含めてもよい。
　上述の圧縮性流体のナビエ・ストークス方程式を離散化した離散化方程式を説明する。一般に、第ｉ軸方向の圧縮性流体のナビエ・ストークス方程式は、式（４－２１）のように表される。また、一般に、非圧縮性流体の連続の式は、式（４－２２）のように表される。なお、本実施形態でも、座標系として直交座標系（カーテシアン座標系）を用いる。式（４－２１）においてμは、粘性係数である。
Figure JPOXMLDOC01-appb-M000053
　本実施形態でも、連続の式である式（４－２２）に速度ベクトルｕを代入しても、誤差Ｄが存在すると仮定する。第１の実施形態と同様に、式（４－２２）の右辺を誤差Ｄとし、両辺にｕｉを乗じると、式（４－２３）が得られる。そして、この式（４－２３）を、ナビエ・ストークス方程式（４－２１）に加えると、式（４－２４）が得られる。
Figure JPOXMLDOC01-appb-M000054
この式（４－２４）を整理すると、保存型である式（４－２５）となる。この式（４－２５）は、左辺の第３項があることを除けば、従来の有限体積法で用いている式と同様であり、第３項に、保存型の移流項における誤差を修正する項（すなわち、－ｕｉＤ）がある。
Figure JPOXMLDOC01-appb-M000055
　この式（４－２５）を、コントロールボリュームで体積積分すると、式（４－２６）が得られる。また、式（４－２２）の右辺を誤差Ｄとした式についても、コントロールボリュームで体積積分すると、式（４－２７）が得られる。第１の実施形態と同様に、Ｖはコントロールボリュームの体積であり、Ｓはコントロールボリュームの表面である。また、式（４－２６）の左辺第２項における太字体のｎとｕは、それぞれコントロールボリュームの表面Ｓの単位法線ベクトルと、速度ベクトルである。単位法線ベクトルｎは、コントロールボリュームの表面Ｓに対して垂直で、外向き、かつ、長さが「１」のベクトルである。
Figure JPOXMLDOC01-appb-M000056
　これらの式（４－２６）、（４－２７）を、例えば、図６のセルＲ１のコントロールボリュームＶ１、境界面Ｅ１ｋ（ｋ∈Ｎ１であり、集合Ｎ１は{２、３、４、１１、１２、１３、１４}である）の面積Ｓ１ｋ、時間Δｔについて離散化し、代数方程式による近似式に変換すると、式（４－２８）、式（４－２９）が得られる。
Figure JPOXMLDOC01-appb-M000057
　式（４－２８）の第３項に式（４－２９）を代入すると、式（４－３０）が得られる。この式（４－３０）に第１の実施形態と同様に、一次風上スキームなどを適用して、速度ベクトルと、圧力とを未知数とする一次方程式にする。
Figure JPOXMLDOC01-appb-M000058
　この式（４－３０）は、セルＲ１に関する式であるが、その他のセル各々についても同様の式を立てることができる。解算出部１０６ａは、全てのセルについての式（４－３０）における各未知数の係数を、行列Ａの要素とし、右辺をベクトルｂの要素とする。
　上述のように、圧縮性流体のナビエ・ストークス方程式を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、本実施形態では、連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、各時間ステップの速度が、連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
　なお、本実施形態では、ナビエ・ストークス方程式と、連続の式と、状態方程式とを連立させて解く場合を説明したが、これらに加えて他の方程式を連立させて解くようにしてもよい。その場合、他の方程式も移流項を含むときは、ナビエ・ストークス方程式と同様に移流項の誤差Ｄを修正する項を含めて離散化するようにしてもよい。これにより、該方程式についても誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
［第３の実施形態］
　以下、図面を参照して、本発明の第３の実施形態について説明する。第１および第２の実施形態では、移流項を含む方程式として、ナビエ・ストークス方程式を挙げ、これを解く流体シミュレーション装置１０、１０ａを説明した。本実施形態では、移流項を含む方程式として連続体の熱エネルギー保存の式を挙げ、これを有限体積法により解く連続体シミュレーション装置１０ｂを説明する。なお、本実施形態の連続体シミュレーション装置１０ｂは、シミュレーションする領域における各時間ステップの速度ベクトルは、外部から入力された値を用いる。この外部から入力された速度ベクトルはシミュレーションにより、予め算出されたものでもよいし、観測データ、実験データであってもよい。
　図１３は、この発明の第３の実施形態による連続体シミュレーション装置１の構成を示す概略ブロック図である。本実施形態における連続体シミュレーション装置１０ｂは、セル構成設定部１０１、セル構成記憶部１０２、速度場設定部１０３ｂ、速度場記憶部１０４ｂ、速度場取得部１０５ｂ、解算出部１０６ｂ、解記憶部１０７ｂを含んで構成される。同図において、図５の各部に対応する部分には、同一の符号（１０１、１０２）を付し、説明を省略する。
　速度場設定部１０３ｂは、外部から入力された各セルの速度ベクトルを、速度場記憶部１０４ｂに記憶させる。なお、本実施形態では、速度ベクトルとして、シミュレーション対象となる全ての時間ステップにおける速度ベクトルが入力される。速度場記憶部１０４ｂは、各セルの速度ベクトルを記憶する。速度場取得部１０５ｂは、解算出部１０６ｂが時間ステップｎに関する解を算出するときは、時間ステップｎ－１における各セルの速度ベクトルを、速度場記憶部１０４ｂから読み出す。
　解算出部１０６ｂは、速度場取得部１０５ｂが速度場記憶部１０４ｂから読み出した時間ステップｎ－１における各セルの速度ベクトルと、セル構成記憶部１０２が記憶するセル情報と、解記憶部１０７ｂが記憶する時間ステップｎ－１における解とを用いて、対象時間ステップｎにおける熱エネルギー保存の式の解を算出する。解算出部１０６ｂは、算出した解を、対象時間ステップｎの解として解記憶部１０７ｂに記憶させる。
　なお、解算出部１０６ｂは、解を算出する際に、式（４－３１）に示すエネルギー保存の式を離散化した離散化方程式を用いる。式（４－３１）に示すように、エネルギー保存の式には、左辺の第２項に移流項がある。密度ρは一定ではないため、解算出部１０６ｂは、第２の実施形態におけるナビエ・ストークス方程式の離散化と同様にして、式（４－３１）を離散化した離散化方程式を用いる。
Figure JPOXMLDOC01-appb-M000059
　上述のように、エネルギー保存の式を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、本実施形態では、連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、外部から入力された各時間ステップの速度が、離散化された連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
　なお、連続体中の物質の移流拡散現象を記述した移流拡散方程式も、式（４－３２）に示すように物質の濃度が移流される移流項を有する。このため、この移流拡散方程式を、エネルギー保存の式と同様に離散化した離散化方程式を解くようにしてもよい。
Figure JPOXMLDOC01-appb-M000060
　また、連続体における応力場を記述した応力場の方程式も、式（４－３３）に示すように連続体の変形速度が移流される移流項を有する。このため、この応力場の方程式を、ナビエ・ストークス方程式と同様に離散化した離散化方程式を解くようにしてもよい。
Figure JPOXMLDOC01-appb-M000061
　なお、連続体の移動・変形を伴う現象を表す方程式には、移流項を含む物質微分が用いられる。該方程式を、上述の各実施形態と同様に離散化した離散化方程式を解くようにしてもよい。
［第４の実施形態］
　以下、図面を参照して、本発明の第４の実施形態について説明する。第４の実施形態では、複数の移流項を含む方程式を有限体積法により解く流体シミュレーション装置１０ｃを説明する。なお、本実施形態では、移流項を持つ方程式を複数同時に解くシミュレーション装置の一例として、ナビエ・ストークス方程式と、熱エネルギー保存の式とを同時に解く流体シミュレーション装置１０ｃを説明する。
　図１４は、この発明の第４の実施形態による流体シミュレーション装置１０ｃの構成を示す概略ブロック図である。本実施形態における流体シミュレーション装置１０ｃは、セル構成設定部１０１、セル構成記憶部１０２、初期条件設定部１０３、速度場記憶部１０４、速度場取得部１０５、解算出部１０６ｃ、解記憶部１０７ｃを含んで構成される。同図において、図５の各部に対応する部分には、同一の符号（１０１～１０５）を付し、説明を省略する。
　解算出部１０６ｃは、速度場取得部１０５が速度場記憶部１０４から読み出した時間ステップｎ－１における各セルの速度ベクトルと、セル構成記憶部１０２が記憶するセル情報と、解記憶部１０７ｃが記憶する時間ステップｎ－１における解とを用いて、対象時間ステップｎにおけるナビエ・ストークス方程式、連続の式、熱エネルギー保存の式の解を算出する。解算出部１０６ｃは、算出した解を、対象時間ステップｎの解として解記憶部１０７ｃに記憶させる。なお、解算出部１０６ｃは、算出した解のうち速度については、速度場記憶部１０４に記憶させる。
　なお、解算出部１０６ｃは、第１の実施形態におけるナビエ・ストークス方程式および連続の式を離散化した離散化方程式と、第３の実施形態における熱エネルギー保存の式を離散化した離散化方程式を用いる。
　なお、本実施形態では、移流項を持つ方程式を複数同時に解くシミュレーション装置の一例として、ナビエ・ストークス方程式と、熱エネルギー保存の式とを同時に解く場合を示したがその他の組み合わせであってもよい。
　本実施形態においても、移流項を持つ方程式各々を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、第１から第３の実施形態と同様に連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、算出した各時間ステップの速度が、離散化された連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
　次に、本発明を用いたシミュレーション結果の例として、第３の実施形態の連続体シミュレーション装置１０ｂにより温度分布を算出した例を説明する。図１５は、シミュレーション対象の速度場を示すコンター図（等高線図）である。すなわち、図１５は、速度場設定部１０３ｂが、速度場記憶部１０４ｂに記憶させる各セルの速度ベクトルの大きさのコンター図である。なお、速度場は、時間変化がないものとして計算を行った。図１５において、横軸はｘ軸、縦軸はｙ軸であり、該速度場の大きさはｘ軸方向が４ｍ、ｙ軸方向が１ｍである。図１５のコンター図は、流体が一様に１ｍ／ｓで左から右へ（ｘ軸の正方向）と流れているが、速度場中に「ｅｒｒ」という形状において１．２ｍ／ｓで左から右へと流れている領域があることを示している。なお、ｙ軸方向の速度は０である。
　図１６は、図１５に示す速度場における速度の発散の大きさを示すコンター図である。図１６においても、横軸はｘ軸、縦軸はｙ軸である。速度の発散は、式（２－３）に示した連続の式の左辺と一致する。すなわち、図１５の速度場が連続の式を完全に満たしていれば、図１６のコンター図は一様に「０」となるはずである。しかし、図１５に示す速度場では、速度が１．２ｍ／ｓの「ｅｒｒ」という形状の領域があるため、該領域の左右端において「０」とはならない。
　図１７は、図１５に示す速度場の温度分布を、従来の方法を用いて算出した場合の例を示すコンター図である。なお、温度分布の初期値は、一様に１００℃としており、図１７に示す温度分布は、定常状態における温度分布である。図１７において、横軸はｘ軸、縦軸はｙ軸であり、等高線は８３℃から１０７℃を１３段階に分割した線である。図１７に示すように、従来の方法にて算出した場合、温度分布の初期値が一様に１００℃であるにも関わらず、「ｅｒｒ」という形状の左右端から下流に、１００℃を超える、あるいは、下回る領域が出来てしまう。これは、「ｅｒｒ」という形状の左右端において、速度の発散が「０」になっていないため、これらの箇所に擬似的な熱のシンク／ソースが発生しているからである。
　図１８は、図１５に示す速度場の温度分布を、第３の実施形態の連続体シミュレーション装置１０ｂにて算出した場合の例を示すコンター図である。図１８の例も、温度分布の初期値は、一様に１００℃としており、定常状態における温度分布である。図１８においても、横軸はｘ軸、縦軸はｙ軸であり、等高線は８３℃から１０７℃を１３段階に分割した線である。図１８に示すように、連続体シミュレーション装置１０ｂにより算出された温度分布は、擬似的な熱のシンク／ソースが発生していないため、一様に１００℃となる。一様に１００℃の流体が流れていれば、時間が経過しても一様に１００℃であるので、連続体シミュレーション装置１０ｂにより良好なシミュレーション結果が得られていることがわかる。
　また、図５、図１２、図１３、図１４における各部の機能もしくはその一部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより流体シミュレーション装置１０、１０ａ、連続体シミュレーション装置１０ｂを実現するようにしてもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。
　また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
　また、本発明のシミュレーション装置は、例えば、特開２０１２－１３２８２４号公報に記載されているような発生源推定装置に用いることができる。この発生源推定装置では、複数個所に設置された観測器の位置情報および計測した濃度情報を入手し、これらの情報に基づき、予め設定されている複数の仮想放出地点の中から放出地点を推定している。このとき、仮想放出地点の一つ一つについて、その仮想放出地点が発生源であり、単位量を放出したときに、各観測器でどのような濃度が測定されるか（該文献における影響関数）を推定し、その推定結果を用いている。この濃度の測定値の推定に、仮想放出地点の設定方法によらず、上述の第３の実施形態で示した移流拡散方程式を解く連続体シミュレーション装置を用いることができる。あるいは、第４の実施形態で示したようにナビエ・ストークス方程式と移流拡散方程式とを同時に解く流体シミュレーション装置を用いることができる。それにより、各観測器で測定される濃度を精度良く推定できるので、発生源も精度よく推定することができる。なお、発生源は、ガス、液体を問わない。
　以上、この発明の実施形態に対して図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、この発明の要旨を逸脱しない範囲の設計変更等も含まれる。
　本出願を詳細にまた特定の実施態様を参照して説明したが、本発明の精神と範囲を逸脱することなく様々な変更や修正を加えることができることは当業者にとって明らかである。
　本出願は、２０１２年１０月３１日出願の日本特許出願（特願２０１２-２４１２６０）に基づくものであり、その内容はここに参照として取り込まれる。
　１０、１０ａ、１０ｃ…流体シミュレーション装置
　１０ｂ…連続体シミュレーション装置
　１０１…セル構成設定部
　１０２…セル構成記憶部
　１０３、１０３ａ…初期条件設定部
　１０３ｂ…速度場設定部
　１０４、１０４ｂ…速度場記憶部
　１０４ａ…密度・速度場記憶部
　１０５、１０５ｂ…速度場取得部
　１０５ａ…密度・速度場取得部
　１０６…速度・圧力算出部
　１０６ａ、１０６ｂ、１０６ｃ…解算出部
　１０７…圧力場記憶部
　１０７ｂ、１０７ｃ…解記憶部

Claims (9)

  1. 　移流項を含む方程式の解を算出するシミュレーション装置であって、
    　速度場を取得する速度場取得部と、
    　前記方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部と
    　を具備し、
    　前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであること
    　を特徴とするシミュレーション装置。
  2. 　前記移流項の保存型の誤差を修正する項が、前記速度場が連続の式を満たさないことにより生じる誤差を修正する項である、請求項１に記載のシミュレーション装置。
  3. 　前記誤差を修正する項が、前記連続の式に対して、前記方程式に対する積分と同じ積分を少なくとも行うことで得られる項である、請求項２に記載のシミュレーション装置。
  4. 　前記解算出部は、少なくとも連続の式を用いて、次の時間ステップにおける速度場を算出し、
    　前記速度場取得部が取得する速度場は、前記解算出部が算出した速度場である、
    　請求項２に記載のシミュレーション装置。
  5. 　前記速度場取得部は、予め記憶された速度場を取得する、請求項２に記載のシミュレーション装置。
  6. 　移流項を含む方程式の解を算出するシミュレーション方法であって、
    　速度場を取得する第１の過程と、
    　移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記第１の過程で取得した速度場を適用して、前記離散化方程式における解を算出する第２の過程と
    　を有し、
    　前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであること
    　を特徴とするシミュレーション方法。
  7. 　コンピュータを、
    　速度場を取得する速度場取得部、
    　移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部
    　として機能させるためのプログラムであって、
    　前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであること
    　を特徴とするプログラム。
  8. 　複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション装置であって、
    　前記観測器の位置情報および計測した濃度情報を入手する観測情報入手手段と、
    　拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定手段と、
    　前記観測器の位置情報、計測した濃度情報、および前記仮想放出地点の情報に基づいて前記仮想放出地点での拡散物質の濃度の推定をする請求項１に記載の装置と、
    　を有するシミュレーション装置。
  9. 　複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション方法であって、
    　前記観測器の位置情報および計測した濃度情報を入手する観測情報入手ステップと、
    　拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定ステップと、
    　前記観測器の位置情報、計測した濃度情報および前記仮想放出地点の情報に基づいて請求項６に記載の方法によって前記仮想放出地点での拡散物質の濃度の推定をするステップと、
    　を有するシミュレーション方法。
PCT/JP2013/079175 2012-10-31 2013-10-28 シミュレーション装置、シミュレーション方法およびプログラム WO2014069421A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN201380056818.9A CN104769593B (zh) 2012-10-31 2013-10-28 仿真装置、仿真方法以及程序
EP13852061.4A EP2916247A4 (en) 2012-10-31 2013-10-28 SIMULATION DEVICE, SIMULATION METHOD, AND PROGRAM
JP2014544502A JP6164222B2 (ja) 2012-10-31 2013-10-28 シミュレーション装置、シミュレーション方法およびプログラム
US14/699,858 US20150234784A1 (en) 2012-10-31 2015-04-29 Simulation device, simulation method, and program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2012-241260 2012-10-31
JP2012241260 2012-10-31

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US14/699,858 Continuation US20150234784A1 (en) 2012-10-31 2015-04-29 Simulation device, simulation method, and program

Publications (1)

Publication Number Publication Date
WO2014069421A1 true WO2014069421A1 (ja) 2014-05-08

Family

ID=50627326

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2013/079175 WO2014069421A1 (ja) 2012-10-31 2013-10-28 シミュレーション装置、シミュレーション方法およびプログラム

Country Status (5)

Country Link
US (1) US20150234784A1 (ja)
EP (1) EP2916247A4 (ja)
JP (1) JP6164222B2 (ja)
CN (1) CN104769593B (ja)
WO (1) WO2014069421A1 (ja)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017162269A (ja) * 2016-03-10 2017-09-14 ソニー株式会社 情報処理装置、電子機器、情報処理方法、及びプログラム
CN108805949B (zh) * 2018-05-25 2022-04-12 杭州晟视科技有限公司 一种血流速度的修正方法、装置、终端和计算机存储介质
CN112560326B (zh) * 2019-09-26 2023-07-18 腾讯科技（深圳）有限公司 压力场的确定方法及装置
CN111930491B (zh) * 2020-09-29 2020-12-25 中国人民解放军国防科技大学 一种全局通信优化加速方法、装置和计算机设备
CN115169260A (zh) * 2022-07-06 2022-10-11 中国南方电网有限责任公司超高压输电公司检修试验中心 灭弧室直流燃弧耐受能力的判断方法、装置和计算机设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11232250A (ja) * 1998-02-10 1999-08-27 Mitsubishi Electric Corp ３次元流れの有限要素解析法、解析装置、成形品を製造する製造方法および解析法のプログラムを記録した媒体
JP2003255055A (ja) * 2002-03-01 2003-09-10 Hitachi Eng Co Ltd 大気汚染発生源推定の環境シミュレーション方法およびシステム
WO2010150758A1 (ja) * 2009-06-25 2010-12-29 旭硝子株式会社 物理量計算方法、数値解析方法、物理量計算プログラム、数値解析プログラム、物理量計算装置及び数値解析装置
JP2012132824A (ja) 2010-12-22 2012-07-12 Mitsubishi Heavy Ind Ltd 拡散物質の発生源推定装置および発生源推定方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7117138B2 (en) * 2003-03-14 2006-10-03 Seiko Epson Corporation Coupled quadrilateral grid level set scheme for piezoelectric ink-jet simulation
JP2006172136A (ja) * 2004-12-15 2006-06-29 Canon Inc 情報処理装置、情報処理方法
EP2481550B1 (en) * 2007-07-02 2014-10-08 MAGMA Giessereitechnologie GmbH Method for describing the statistical orientation distribution of particles in a simulation of a mould filling process and computer software product including a software code for performing said method.
US8428922B2 (en) * 2010-02-05 2013-04-23 Seiko Epson Corporation Finite difference level set projection method on multi-staged quadrilateral grids
KR101113301B1 (ko) * 2010-03-23 2012-02-24 한국과학기술원 부피비를 이용한 다종 유체 해석 방법 및 기록 매체
US20130013277A1 (en) * 2011-07-08 2013-01-10 Jiun-Der Yu Ghost Region Approaches for Solving Fluid Property Re-Distribution

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11232250A (ja) * 1998-02-10 1999-08-27 Mitsubishi Electric Corp ３次元流れの有限要素解析法、解析装置、成形品を製造する製造方法および解析法のプログラムを記録した媒体
JP2003255055A (ja) * 2002-03-01 2003-09-10 Hitachi Eng Co Ltd 大気汚染発生源推定の環境シミュレーション方法およびシステム
WO2010150758A1 (ja) * 2009-06-25 2010-12-29 旭硝子株式会社 物理量計算方法、数値解析方法、物理量計算プログラム、数値解析プログラム、物理量計算装置及び数値解析装置
JP2012132824A (ja) 2010-12-22 2012-07-12 Mitsubishi Heavy Ind Ltd 拡散物質の発生源推定装置および発生源推定方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
H. K. VERSTEEG; W. MALALASEKERA; MATSUSHITA YOSUKE; YASUHIRO SAITO; HIDEYUKI AOKI: "Numerical Fluid Dynamics, 2nd ed.", 30 May 2011, MORIKITA PUBLISHING CO., LTD., article "An Introduction to Computational Fluid Dynamics"
MORINISHI YOHEI: "Special Reviews by RYUMON-Prize Award", FULLY CONSERVATIVE HIGHER ORDER FINITE DIFFERENCE SCHEMES FOR INCOMPRESSIBLE FLUID ANALYSIS, vol. 19, 2000, pages 164 - 170
MORINISHI, Y.; LUND, T. S.; VASILYEV, O. V.; MOIN, P.: "Fully Conservative Higher Order Finite Difference Schemes for incompressible flow", J. COMPUT. PHYS., vol. 143, 1998, pages 90 - 124, XP055257887, DOI: doi:10.1006/jcph.1998.5962
See also references of EP2916247A4
YOHEI MORINISHI: "Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow", NAGARE, vol. 28, no. 3, June 2009 (2009-06-01), pages 197 - 200, XP055257887 *

Also Published As

Publication number Publication date
US20150234784A1 (en) 2015-08-20
EP2916247A4 (en) 2016-04-20
CN104769593A (zh) 2015-07-08
CN104769593B (zh) 2018-09-21
JP6164222B2 (ja) 2017-07-19
JPWO2014069421A1 (ja) 2016-09-08
EP2916247A1 (en) 2015-09-09

Similar Documents

Publication Publication Date Title
Ballarin et al. POD–Galerkin monolithic reduced order models for parametrized fluid‐structure interaction problems
JP6164222B2 (ja) シミュレーション装置、シミュレーション方法およびプログラム
Réthoré et al. A discrete force allocation algorithm for modelling wind turbines in computational fluid dynamics
JP2008165804A (ja) 流れのシミュレーション計算方法およびシステム
Xie et al. An accurate and robust HLLC‐type Riemann solver for the compressible Euler system at various Mach numbers
Abdo et al. Steady supercritical flow in a straight-wall open-channel contraction
Altuna et al. Application of a fast loosely coupled fluid/solid heat transfer method to the transient analysis of low-pressure-turbine disk cavities
JP2006258391A (ja) 浮力を伴う乱流の流体的及び熱的諸特性の推定方法及び推定プログラム
Gallerano et al. Central WENO scheme for the integral form of contravariant shallow‐water equations
Coppola‐Owen et al. A free surface finite element model for low Froude number mould filling problems on fixed meshes
Jarosch Icetools: A full Stokes finite element model for glaciers
Wang et al. Numerical simulation of one-dimensional two-phase flow using a pressure-based algorithm
Constantine et al. A hybrid collocation/Galerkin scheme for convective heat transfer problems with stochastic boundary conditions
Dauricio et al. A wall model for external laminar boundary layer flows applied to the Wall-Modeled LES framework
Moguen et al. Pressure–velocity coupling for unsteady low Mach number flow simulations: An improvement of the AUSM+-up scheme
Wu et al. Simulation of thermal flow problems via a hybrid immersed boundary-lattice Boltzmann method
Liberson et al. Numerical solution for the boussinesq type models with application to arterial flow
Yang et al. Numerical investigation on performance of three solution reconstructions at cell interface in DVM simulation of flows in all Knudsen number regimes
Karimian et al. A thermal periodic boundary condition for heating and cooling processes
Zhang et al. A conservative pressure based solver with collocated variables on unstructured grids for two-fluid flows with phase change
Spiering Coupling of TAU and TRACE for parallel accurate flow simulations
Yang et al. An immersed boundary method based on parallel adaptive Cartesian grids for high Reynolds number turbulent flow
Zdanski et al. A numerical method for simulation of incompressible three-dimensional Newtonian and non-Newtonian flows
Darbandi et al. Efficient multilevel restriction–prolongation expressions for hybrid finite volume element method
Brown et al. Efficient numerical differentiation of implicitly-defined curves for sparse systems

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 13852061

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2014544502

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 2013852061

Country of ref document: EP

NENP Non-entry into the national phase

Ref country code: DE