JP2014119882A - Analysis method, analysis device, and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To highly accurately perform numerical analysis in an advective diffusion question.SOLUTION: An analysis method includes: a reverse tracking step of, for each calculation lattice point in calculation lattices set in an analysis target field, retroactively tracking a locus of a virtual particle which has a variable field value and advectively flows to the calculation lattice point; and a field value acquisition step of, for each calculation lattice point, on the basis of a field value of a termination position of the reverse tracking and a result of the reverse tracking, acquiring a field value of the virtual particle when the virtual particle reaches the calculation lattice point from the termination position of the reverse tracking in a diffuse manner, as a field value at the calculation lattice point.

Description

  The present invention relates to an analysis method, an analysis apparatus, and a program.

  There is an advection diffusion problem as one of the main analysis objects in the numerical analysis field. The advection diffusion problem has been applied to a wide range including heat transfer and mixing problems involving fluids, and the sophistication and efficiency of numerical analysis methods in the advection diffusion problem are regarded as important from the viewpoint of industrial use. Yes. On the other hand, the accuracy required for the analysis of the advection-diffusion problem is generally very high, and it has been difficult to correctly predict the phenomenon except for a few ideal cases.

There is a lattice method as a numerical analysis method frequently used in the advection diffusion problem (for example, Patent Document 1). The lattice method is a general term for an analysis method that captures a physical phenomenon by Euler's way of thinking and discretely expresses the phenomenon by a calculation lattice provided in space.
Further, there is a particle method as a numerical analysis method compared with the lattice method in the advection diffusion problem (for example, Patent Document 2). The particle method is an analysis method that incorporates a Lagrangian way of thinking.

JP 2004-13442 A JP 7-334484 A

In the numerical analysis of the advection diffusion problem using the grid method, the accuracy of the analysis result decreases due to the occurrence of numerical errors due to the numerical diffusion. Numerical diffusion here is a pseudo-diffusion that is added with the discretization of the convection term in the advection-diffusion equation. Numerical diffusion is more dominant than actual molecular diffusion. The analysis results are not uncommon.
On the other hand, in the particle method, since a phenomenon is expressed using a particle group, numerical diffusion does not fundamentally occur. However, although the particle method is excellent in numerical analysis of a phenomenon involving a free surface represented by a gas-liquid two-phase flow, it does not have a means for quantitatively evaluating a physical quantity distribution in a space such as a gradient or a Laplacian. It has been pointed out that molecular diffusion cannot be accurately evaluated in the advection diffusion problem for the interparticle interaction model used as a solution.

  The present invention has been made in view of such circumstances, and an object of the present invention is to provide an analysis method, an analysis apparatus, and a program capable of performing a numerical analysis with higher accuracy in an advection diffusion problem.

  The present invention has been made to solve the above-described problems, and an analysis method according to an aspect of the present invention provides a variable field value for each calculation grid point in a calculation grid set as a field to be analyzed. A reverse tracking step of performing reverse tracking for tracking the trajectory of virtual particles that advect to the calculation grid point by tracing back the time, and a field value at the end position of the reverse tracking for each of the calculation grid points; Based on the result of the reverse tracking, the field value of the virtual particle is obtained as the field value at the calculation grid point when the calculation grid point is reached from the end position of the reverse tracking while diffusing. And a field value acquisition step.

  Further, the analysis method according to one aspect of the present invention is the above-described analysis method, wherein in the reverse tracking step, a trajectory of the virtual particle is obtained for each of the calculation lattice points, and in the field value acquisition step, A field value at each of the calculation lattice points is obtained from simultaneous equations obtained based on the locus of the virtual particles and the diffusion equation for each of the calculation lattice points.

  Further, the analysis method according to an aspect of the present invention is the above-described analysis method, wherein in the reverse tracking step, a trajectory of the virtual particle is obtained for each of the calculation lattice points, and the field value acquisition step includes For each of the calculation lattice points, integrating a diffusion equation along the trajectory of the virtual particle, and calculating a field value at each of the calculation lattice points, and a field value calculation step A convergence determination step for determining whether or not the calculated field value has converged, and the convergence determination step determines that the field value calculated in the field value calculation step has converged Until this is done, the process in the field value calculation step is repeated.

  The analysis method according to an aspect of the present invention is the above-described analysis method, wherein the field value acquisition step includes, for each of the calculation grid points, the virtual particles that reach the calculation grid point are the backtracked A field value setting step of setting a field value at the end position of the field to the calculation grid point, and for each of the calculation grid points, the virtual particles reach the calculation grid point from the end position of the reverse tracking And a field value updating step for updating the field value at the calculation grid point in accordance with the residence time required.

  Further, the analysis apparatus according to one aspect of the present invention provides a trajectory of virtual particles having a variable field value and advancing to the calculation grid point for each of the calculation grid points in the calculation grid set in the field to be analyzed. , While diffusing, based on the field value at the end position of the reverse tracking and the result of the reverse tracking for each of the calculation grid points A field value acquisition unit that acquires a field value in the virtual particle when reaching the calculation grid point from the end position of the reverse tracking, as a field value in the calculation grid point; To do.

  In addition, the program according to one aspect of the present invention provides a computer with a virtual particle having a variable field value and advancing to a calculation grid point for each calculation grid point in a calculation grid set as a field to be analyzed. A backtracking step for performing backtracking that traces the trajectory back in time, and for each of the calculation grid points, diffusion is performed based on the field value at the end position of the backtracking and the result of the backtracking. In order to execute the field value acquisition step of acquiring the field value of the virtual particle as the field value at the calculation grid point when reaching the calculation grid point from the end position of the reverse tracking It is a program.

  According to the present invention, numerical analysis can be performed with higher accuracy in a short time and with a smaller storage capacity than in the conventional method in the advection diffusion problem.

It is a schematic block diagram which shows the function structure of the analyzer in the 1st Embodiment of this invention. It is explanatory drawing which shows the example of the reverse tracking of the locus | trajectory of the virtual particle which the reverse tracking part in the embodiment performs. In the same embodiment, it is a flowchart which shows the procedure of the process which an analyzer performs. It is a schematic block diagram which shows the function structure of the analyzer in the 2nd Embodiment of this invention. In the same embodiment, it is a flowchart which shows the procedure of the process which an analyzer performs. It is explanatory drawing which shows the example of the scalar quantity calculation which the field value calculation part in the embodiment performs. It is a schematic block diagram which shows the function structure of the analyzer in the 3rd Embodiment of this invention. It is a flowchart which shows the procedure of the process which an analysis apparatus performs in the same embodiment. An example of a typical T-type mixer used for fluid mixing is shown. It is explanatory drawing which shows the example of the reverse tracking of the virtual particle in the analysis of the concentration field in the 1st analysis example of this invention. It is explanatory drawing which shows the result compared with the lattice method, changing the example of the mixing degree MI in the target plane in the 1st analysis example of this invention, changing the space | interval between calculation lattice points. It is explanatory drawing which shows the example of the analysis result of the density | concentration distribution in the maximum spatial resolution of the lattice method and the method in the 1st analysis example of this invention. It is explanatory drawing which shows the result of having compared the example of the calculation result of the mixing degree MI in the target plane A in the 2nd analysis example of this invention with the lattice method, changing the space | interval between calculation lattice points. It is explanatory drawing which shows the example of the analysis result of a density | concentration distribution with the method in the grid method and the 1st analysis example of this invention. It is explanatory drawing which shows the comparative example of the analysis result in the 1st analysis example of this invention, the analysis result in a 2nd analysis example, and the analysis result by a lattice method.

<First Embodiment>
Embodiments of the present invention will be described below with reference to the drawings.
In the description of the specification, the notation of an arrow indicating a vector or a matrix, and the superscript and subscript after the second stage are omitted.
FIG. 1 is a schematic block diagram showing a functional configuration of an analysis apparatus according to the first embodiment of the present invention. In the figure, the analysis apparatus 100 includes a precondition acquisition unit 110, a reverse tracking unit 120, a field value acquisition unit 130, and an analysis result output unit 140. The field value acquisition unit 130 includes a field value setting unit 131 and a field value update unit 132.

  The analysis device 100 performs numerical analysis to obtain a field value in the advection diffusion problem. The analysis device 100 is constructed by a computer, for example. Alternatively, the analysis apparatus 100 may be constructed by a device other than a computer, such as a dedicated circuit.

The precondition acquisition unit 110 acquires a precondition for performing numerical analysis. For example, the precondition acquisition unit 110 indicates, as preconditions for performing the numerical analysis, the size of the calculation grid set in the space to be analyzed, the interval between the calculation grid points, and the advection velocity field (the speed at which the field value advects). Get speed field data.
Various forms can be used as a form in which the precondition acquisition unit 110 acquires a precondition for performing numerical analysis. Alternatively, the precondition acquisition unit 110 may store in advance preconditions for performing numerical analysis, or when the analysis apparatus 100 performs numerical analysis, it receives user input of preconditions for performing numerical analysis. You may do it. In addition, the precondition acquisition unit 110 may obtain the advection velocity field by analysis, or may obtain data of the advection velocity field that has been analyzed in advance.

The reverse tracking unit 120 reversely tracks the trajectory of the virtual particles that advect (also referred to as convection or transport) to the calculation grid point for each calculation grid point in the calculation grid set in the analysis target field. In particular, the reverse tracking unit 120 reversely tracks the trajectory of the virtual particle to a region where the field value is known, such as an inflow boundary. Further, the reverse tracking unit 120 obtains the residence time required for the virtual particles to reach the calculation lattice point from the reverse tracking end position.
The process performed by the reverse tracking unit 120 corresponds to an example of a process in the reverse tracking step.
The virtual particle here is a virtual particle (fluid particle) having a variable field value and advancing to a calculation lattice point and having no mass (inertia). In addition, reverse tracking here refers to tracking a trajectory back in time. Accordingly, in the reverse tracking of the trajectory of the virtual particle, the trajectory is traced in the direction opposite to the direction in which the virtual particle flows.

FIG. 2 is an explanatory diagram illustrating an example of reverse tracking of the trajectory of the virtual particle performed by the reverse tracking unit 120. In the figure, a point P11 indicates a calculation lattice point, and a line L11 indicates a locus of virtual particles.
While the virtual particles have reached the calculation grid point (point P11) by advection, the reverse tracking unit 120 has a trajectory (line L11) in a direction opposite to the direction in which the virtual particles advect as indicated by an arrow V11. ).

By tracing back the trajectory of the virtual particle, the field value based on the convection term (also referred to as the advection term) at the calculation lattice point can be obtained. In addition, by tracing back the trajectory of the virtual particles, the diffusion when the virtual particles advect can be calculated. Thereby, the influence of the diffusion term can be reflected on the field value at the calculation lattice point.
Further, since virtual particles are virtual particles, energy advection and diffusion such as heat conduction can be analyzed using virtual particles as well as advection and diffusion by existing particles.

  The field value acquisition unit 130 reaches the calculation grid point from the end position of the reverse tracking while diffusing based on the field value at the end position of the reverse tracking and the result of the reverse tracking for each of the calculation grid points. At this time, the field value of the virtual particle is acquired as the field value at the calculation lattice point. The process performed by the field value acquisition unit 130 corresponds to an example of a process in the field value acquisition step.

The field value setting unit 131 sets, for each calculation lattice point, the field value that the virtual particle that reaches the calculation lattice point has at the end position of the reverse tracking as the field value at the calculation lattice point. More specifically, the reverse tracking unit 120 reversely tracks the trajectory of the virtual particles up to a region where the field value is known, such as an inflow boundary. Then, the field value acquisition unit 130 sets the field value known at the end position of the reverse tracking as the field value at the calculation grid point. Thereby, the field value setting unit 131 sets the field value based on the convection term at the calculation grid point.
The process performed by the field value setting unit 131 corresponds to an example of a process in the field value setting step.

The field value updating unit 132 sets, for each of the calculation grid points, the field value setting unit 131 set according to the residence time required for the virtual particles to reach the calculation grid point from the end position of the reverse tracking. Update the field value at the computational grid point. Thereby, the field value updating unit 132 reflects the influence of the diffusion term on the field value at the calculation lattice point.
The process performed by the field value update unit 132 corresponds to an example of a process in the field value update step.

  The analysis result output unit 140 outputs the field value acquired by the field value acquisition unit 130. Various modes can be used as a mode in which the analysis result output unit 140 outputs the field value. For example, the analysis result output unit 140 may have a display screen, and the obtained field value may be displayed visually, such as a graph display or a numerical display. Alternatively, the analysis result output unit 140 may output (send) data indicating the obtained field value to another device connected to the analysis device.

  Next, the operation of the analysis apparatus 100 will be described with reference to FIG. In the following, a case where the analysis apparatus 100 analyzes a scalar field will be described as an example, but the place where the analysis apparatus 100 performs an analysis is not limited to a scalar field. In addition, the analysis apparatus 100 can set a space of various dimensions as an analysis target. The same applies to other embodiments.

  Hereinafter, when a velocity field and a scalar field are formed in a space to be analyzed, and the scalar quantity constituting the scalar field diffuses while advancing on the velocity field, the problem of obtaining the scalar quantity distribution will be considered. More specifically, the analysis apparatus 100 sets a calculation grid in a space to be analyzed, and obtains a scalar amount at each calculation grid point. In the following description, the total number of calculation grid points is N, and the calculation grid points are identified by a number k (0 ≦ k ≦ N−1).

FIG. 3 is a flowchart showing a procedure of processing performed by the analysis apparatus 100.
In the process shown in the figure, the precondition acquisition unit 110 first acquires preconditions for performing numerical analysis, such as calculation grid setting conditions and advection velocity field data (step S101).
Next, the reverse tracking unit 120 performs reverse tracking of the trajectory of the virtual particle, and the field value setting unit 131 evaluates the convection term based on the result of the reverse tracking (step S102). Hereinafter, the reverse tracking performed by the reverse tracking unit 120 and the evaluation of the convection term performed by the field value setting unit 131 will be described in more detail.

First, a virtual particle is assumed to advect on the velocity field from a region where the scalar amount is known to the calculation lattice point k. The reverse tracking unit 120 calculates the trajectory x _k (t) of the virtual particle by solving the equation (1) in the time direction reverse direction.

However, t shows time and u shows the speed in a speed field. Note that the velocity field may be stationary or may change according to time. Furthermore, the scalar field may change according to the time as the speed field changes. That is, the analysis apparatus 100 can analyze both a stationary scalar field and a non-stationary scalar field. When analyzing an unsteady scalar field, the analysis apparatus 100 obtains the distribution of the scalar quantity at a certain time (for example, the time t ₀ when the virtual particle reaches the calculation lattice point).

  As a method of solving the equation (1), it is conceivable to solve the equation discretely. At this time, various schemes such as Euler explicit method, complete implicit method, Crank Nicholson method, Runge-Kutta method, Runge-Kutta-Ferberg method or Adams Bash force method can be used as the discretization scheme performed by the inverse tracking unit 120.

  Further, when the speed u is defined discretely, it is conceivable to use an interpolation value of the speed u from a nearby speed definition point. At that time, various schemes such as primary accuracy interpolation, bilinear interpolation, trilinear interpolation, secondary accuracy interpolation, or cubic accuracy interpolation can be used as the interpolation scheme performed by the inverse tracking unit 120.

  In addition, when the locus | trajectory obtained by calculation error turns into an unrealistic locus | trajectory, you may make it the reverse tracking part 120 take a means (for example, correction | amendment of a speed field) which prevents it. For example, when the virtual particle reaches the boundary (wall surface) where the velocity condition is non-slip due to the influence of numerical error, the reverse tracking unit 120 corrects the velocity used for reverse tracking to avoid stopping the virtual particle You may make it do.

The reverse tracking unit 120 performs the reverse tracking until the virtual particle reaches a region where the scalar amount is known (for example, the inflow boundary), so that the field value setting unit 131 is a scalar of the calculation lattice point k when diffusion is not considered. The amount can be acquired. An example of an area where the scalar quantity is known is an inflow boundary in which the scalar quantity is given in advance.
The reverse tracking unit 120 performs the reverse tracking described above for each calculation grid point. Then, the field value setting unit 131 obtains a scalar distribution when diffusion is not considered, based on the result of reverse tracking at each calculation grid point.

  Note that the reverse tracking unit 120 sets a region where the scalar amount is known, such as an inflow boundary, as the reverse tracking end position, and the virtual particle scalar amount at the reverse tracking end position and the virtual particle reaches the calculation lattice point from the reverse tracking end position. What is necessary is just to acquire time (dwelling time) required to do. Therefore, the reverse tracking unit 120 does not necessarily need to obtain a specific trajectory such as obtaining the coordinates of the trajectory of the virtual particle.

Next, the field value update unit 132 evaluates the diffusion term using the model (step S103). Hereinafter, the evaluation of the diffusion term performed by the field value update unit 132 will be described in more detail.
The field value update unit 132 uses the scalar distribution obtained by the field value setting unit 131 when the diffusion is not considered as a value in the initial state (ξ = 0), and the model expression shown in the equation (2) is ξ = [0 , 1] is integrated numerically.

Where φ represents a scalar field. Further, τ represents the residence time of virtual particles. That is, τ indicates the time required for the virtual particle to reach the calculation lattice point k from the area where the scalar quantity is known. D indicates a diffusion coefficient, and ▽ ² indicates a Laplacian operator. Ξ is a temporarily set variable.
The field value updating unit 132 obtains a scalar distribution in consideration of diffusion in the integration result ξ = 1. In the calculation of the scalar quantity using the equation (2), the initial distribution is a scalar distribution when diffusion is not taken into consideration, the diffusion in the locus of the virtual particles is approximated by the diffusion at the calculation lattice point k, and the retention of the virtual particles Simulates diffusion corresponding to time τ. In this respect, the field value updating unit 132 evaluates the diffusion term using a model.

  In addition, as the discretization scheme performed by the field value updating unit 132 when performing the integration with respect to the equation (2), various methods such as the Euler explicit method, the complete implicit method, the Crank Nicholson method, the Runge-Kutta method, the Runge-Kutta-Ferberg method, and the Adams Bash force method are used. Any discretization scheme can be used.

After step S103, the analysis result output unit 140 outputs a scalar field as the analysis result (step S104). More specifically, the analysis result output unit 140 outputs a scalar quantity at each calculation grid point.
Thereafter, the process of FIG. 3 is terminated.

As described above, the reverse tracking unit 120 performs the reverse tracking of the trajectory of the virtual particles that advect to the calculation grid point, and the field value acquisition unit 130, when diffusing, reaches the calculation grid point from the end position of the reverse tracking. The field value of the virtual particle is obtained as the field value at the calculation lattice point.
Thereby, in the analysis apparatus 100, since it is not necessary to interpolate the field value which a virtual particle has in a calculation lattice point, numerical diffusion does not arise and it can perform a numerical analysis with higher precision.

Further, the field value setting unit 131 sets, for each calculation grid point, the field value that the virtual particle that reaches the calculation grid point has at the end position of the reverse tracking in the calculation grid point. Then, the field value updating unit 132 calculates the field value at the calculation grid point for each calculation grid point according to the residence time required for the virtual particle to reach the calculation grid point from the end position of the reverse tracking. Update.
Thereby, the field value setting unit 131 can evaluate the convection term by a simple process of setting the field value at the end position of the reverse tracking as the field value at the calculation grid point. In addition, the field value update unit 132 calculates the field value when the virtual particle stays at the calculation lattice point reached by the virtual particle for the above-described residence time by a simple process, for example. Evaluation can be made.

<Second Embodiment>
FIG. 4 is a schematic block diagram showing a functional configuration of the analysis apparatus according to the second embodiment of the present invention. In the figure, the analysis apparatus 200 includes a precondition acquisition unit 110, a reverse tracking unit 220, a field value acquisition unit 230, and an analysis result output unit 140. The field value acquisition unit 230 includes a field value calculation unit 231 and a convergence determination unit 232. 4, parts having the same functions corresponding to the parts in FIG. 1 are denoted by the same reference numerals (110, 140), and description thereof is omitted.

  Similar to the analysis device 100, the analysis device 200 performs numerical analysis to obtain a field value in the advection diffusion problem. However, the specific process performed by the analysis apparatus 200 is different from the process performed by the analysis apparatus 100. The analysis device 200 is constructed by a computer, for example. Alternatively, the analysis device 200 may be constructed by a circuit other than a computer, such as a dedicated circuit.

For each of the calculation grid points in the calculation grid set in the analysis target field, the reverse tracking unit 220 traces the trajectory of the virtual particles having a variable field value and advancing to the calculation grid point back in time. Chase. In particular, the inverse tracking unit 220 specifically obtains the virtual particle trajectory, such as obtaining the coordinates of the virtual particle trajectory.
The process performed by the reverse tracking unit 220 corresponds to an example of a process in the reverse tracking step.

Similar to the field value acquisition unit 130, the field value acquisition unit 230 performs reverse tracking while spreading based on the field value at the end position of the reverse tracking and the result of reverse tracking for each of the calculation grid points. The field value of the virtual particle when the calculation lattice point is reached from the end position is acquired as the field value at the calculation lattice point. However, the specific process performed by the field value acquisition unit 230 is different from the process performed by the field value acquisition unit 130.
The process performed by the field value acquisition unit 230 corresponds to an example of a process in the field value acquisition step.

The field value calculation unit 231 integrates the diffusion equation along the trajectory of the virtual particle for each calculation lattice point, and calculates the field value at each calculation lattice point. The process performed by the field value calculation unit 231 corresponds to an example of a process in the field value calculation step.
The convergence determination unit 232 determines whether or not the field value calculated by the field value calculation unit 231 has converged. The field value calculation unit 231 repeats the calculation of the field value until the convergence determination unit 232 determines that the field value calculated by the field value calculation unit 231 has converged. The process performed by the convergence determination unit 232 corresponds to an example of a process in the convergence determination step.

Next, the operation of the analysis apparatus 200 will be described with reference to FIG.
In the first embodiment, diffusion is not taken into account when advancing virtual particles. For this reason, in the first embodiment, the scalar amount of virtual particles does not change during advection. On the other hand, in the second embodiment, the scalar amount of the virtual particles is changed in consideration of diffusion when the virtual particles advect. Accordingly, in the second embodiment, the scalar quantity of the calculation lattice point k is expressed as φ _k (t), and the scalar quantity of the virtual particle is expressed as φ (tilde) _k (t) to distinguish between them. .

FIG. 5 is a flowchart illustrating a procedure of processing performed by the analysis apparatus 200. Step S201 in FIG. 5 is the same as step S101 in FIG.
After step S201, the reverse tracking unit 220 reverse-tracks the trajectory of the virtual particle and calculates the trajectory (step S202). More specifically, the reverse tracking unit 220 obtains the trajectory x _k (t) of the virtual particle that advects to the calculation lattice point k on the velocity field by solving Equation (1) in the reverse direction in the time direction. Here, the reverse tracking unit 220 may perform the reverse tracking until it reaches the area where the scalar quantity is known, but is not limited to this, and it is only necessary to reversely track the virtual particles over a relatively long time (distance). Therefore, the analysis apparatus 200 can also analyze a field in a system (closed system) without an inflow boundary or a field in a dead water region.

  As a method of solving the equation (1), it is conceivable to solve the equation discretely. At this time, various schemes such as Euler explicit method, complete implicit method, Crank Nicholson method, Runge-Kutta method, Runge-Kutta-Ferberg method or Adams Bash force method can be used as the discretization scheme performed by the inverse tracking unit 220.

Further, when the speed u is defined discretely, it is conceivable to use an interpolation value of the speed u from a nearby speed definition point. At that time, various schemes such as primary accuracy interpolation, bilinear interpolation, trilinear interpolation, secondary accuracy interpolation, or cubic accuracy interpolation can be used as the interpolation scheme performed by the inverse tracking unit 220.
The reverse tracking unit 220 performs the reverse tracking described above for each calculation grid point.

  In addition, when the locus | trajectory obtained by a calculation error turns into an unrealistic locus | trajectory, you may make it the reverse tracking part 220 take a means (for example, correction | amendment of a speed field) which prevents it. For example, when the virtual particle reaches the boundary (wall surface) where the velocity condition is non-slip due to the influence of numerical error, the reverse tracking unit 220 corrects the velocity used for the reverse tracking to avoid stopping the virtual particle. You may make it do.

Next, the field value calculation unit 231 calculates a scalar field based on the trajectory and the diffusion equation (step S203). Hereinafter, the calculation of the scalar field performed by the field value calculation unit 231 will be described in more detail.
First, similarly to the field value setting unit 131 (FIG. 1), the field value calculation unit 231 also obtains a scalar distribution when diffusion is not considered as an initial value for scalar distribution calculation. However, the field value calculation unit 231 can use various scalar distributions as initial values for the scalar distribution calculation without being limited to the scalar distribution when diffusion is not considered. This is because, as will be described below, the field value calculation unit 231 repeats the calculation of the scalar field until the calculation result converges. On the other hand, the field value setting unit 231 uses the scalar distribution when diffusion is not considered as an initial value, whereby the number of iterations until the calculation converges can be reduced, and the processing time and processing load can be reduced.
Next, unlike the case of the first embodiment, the field value calculation unit 231 discretely solves the diffusion equation represented by Equation (3).

However, t shows time. D represents a diffusion coefficient, 、 ² represents a Laplacian operator, φ represents a scalar field, and s represents a generation term.
Specifically, the field value calculation unit 231 integrates the equation (3) from t = t _M to t = t ₀ and diffuses the scalar amount φ (tilde of the virtual particle that reaches the calculation lattice point k while diffusing. ) _k (t ₀ ) is obtained, and the scalar quantity φ _k (t ₀ ) of the calculation lattice point k is obtained by the equation (4).

However, t ₀ is a virtual particles, indicates the time to reach the calculation lattice point k is the inverse tracking start position. T _M indicates the time at which the virtual particle leaves the reverse tracking end position. Here, various positions included in the trajectory can be set as reverse tracking end positions. For example, the position where the virtual particle departs from the area of known scalar quantity such as the inflow boundary may be set as the reverse tracking end position. Alternatively, one point on the trajectory between the area where the scalar quantity is known and the calculation grid point k may be set as the reverse tracking end position.

As the discretization scheme related to time performed by the field value calculation unit 231 when integrating the expression (3), various methods such as the Euler explicit method, the complete implicit method, the Crank Nicholson method, the Runge-Kutta method, the Runge-Kutta-Ferberg method, and the Adams Bash force method are used. Schemes can be used.
At this time, if the scalar quantity at the reverse tracking end position x _k (t _M ) is given by the boundary condition or the like, this is substituted as the initial condition φ (tilde) _k (t _M ) of the scalar quantity of the virtual particle. . On the other hand, when the scalar quantity at the reverse tracking end position x _k (t _M ) is not given, interpolation is performed from a nearby calculation grid point. At this time, various schemes such as first-order accuracy interpolation, bilinear interpolation, trilinear interpolation, second-order accuracy interpolation, or third-order accuracy interpolation can be used as the interpolation scheme performed by the field value calculation unit 231.

  For Laplacian, the value on the locus is used as the scalar quantity at the center position, and the scalar quantity at the calculation grid point is used as the scalar distribution in the vicinity. As a Laplacian discretization scheme performed by the field value calculation unit 231, various schemes such as various Laplacian models, a secondary accuracy central difference method, and a fourth accuracy central difference method can be used. For example, when the secondary accuracy central difference method is applied, the equation (5) is obtained.

Here, e _x , e _y , and e _z represent unit vectors (base vectors) in the x direction, the y direction, and the z direction, respectively.
Further, δx, δy, and δz respectively indicate the width between the center position and the vicinity in the difference method. All of δx, δy, and δz are set to a value equal to or greater than the distance between the calculation grid points. This is to prevent overestimation of Laplacian.

  Φ (bar) (t, x) indicates the interpolation of the scalar quantity from the calculation grid point to the position x. At this time, various schemes such as first-order accuracy interpolation, bilinear interpolation, trilinear interpolation, second-order accuracy interpolation, or third-order accuracy interpolation can be used as the interpolation scheme performed by the field value calculation unit 231. However, from the viewpoint of analysis accuracy, it is preferable to use an interpolation scheme of spatial accuracy second order or higher.

FIG. 6 is an explanatory diagram illustrating an example of scalar amount calculation performed by the field value calculation unit 231. In the figure, a point P21 indicates a calculation lattice point at which the virtual particle reaches, and a line L21 indicates a locus of the virtual particle.
The field value calculation unit 231 uses the point P22 on the locus (line L21) as the center position in the difference method, and calculates the scalar amount of a nearby point such as the point P23 from the surrounding calculation grid points as shown in the figure. Calculate by insertion. Then, the field value calculation unit 231 substitutes the scalar quantity of the obtained nearby points into Expression (5).
In this way, the field value calculation unit 231 obtains a change in the scalar amount of the virtual particles due to diffusion based on the difference in scalar amount between the point on the trajectory and the surroundings. By integrating this along the trajectory, the scalar quantity of the calculation lattice point (point P21) at which the virtual particle has reached can be obtained by the displacement from the scalar quantity at the end position of the reverse tracking.
The field value calculation unit 231 obtains a scalar field by solving Equation (3) for each calculation grid point.

After step S203, the convergence determination unit 232 determines whether the scalar field obtained by the field value calculation unit 231 has converged (step S204). Specifically, the convergence determination unit 232 calculates the magnitude of the difference between the scalar field newly obtained by the field value calculation unit 231 and the scalar field obtained by the field value calculation unit 231 last time. When the magnitude is less than or equal to a predetermined threshold value, it is determined that it has converged.
The convergence determination unit 232 determines the magnitude of the change in the scalar field as the magnitude of the difference between the scalar field newly calculated by the field value calculation unit 231 and the scalar field previously calculated by the field value calculation unit 231. Various values indicating can be used. For example, the convergence determination unit 232 may use a value obtained by adding the absolute value of the change in the scalar amount for each calculation grid point for all the calculation grid points. Alternatively, the convergence determination unit 232 may use the root mean square of the change in scalar quantity for each calculation grid point.
When the field value calculation unit 231 first obtains a scalar field, the convergence determination unit 232 determines that the scalar field obtained by the field value calculation unit 231 has not converged.

If it is determined in step S204 that the scalar field obtained by the field value calculation unit 231 has not converged (step S204: NO), the process returns to step S203.
In this case, in step S203, the field value calculation unit 231 solves Equation (3) for each calculation grid point using the scalar field obtained in the previous process in step S203. Thereby, the field value calculation unit 231 updates the scalar field (scalar amount of each calculation grid point). Hereinafter, until the convergence determination unit 232 determines that the scalar field obtained by the field value calculation unit 231 has converged, that is, until the change in the scalar field becomes sufficiently small, the field value calculation unit 231 determines the obtained scalar field. The process of resolving Equation (3) for each calculation grid point using the field is repeated.

On the other hand, if it is determined in step S204 that the scalar field obtained by the field value calculation unit 231 has converged (step S204: YES), the process proceeds to step S205.
Step S205 is the same as step S104 in FIG. Thereafter, the process of FIG.

As described above, the reverse tracking unit 220 backtracks the trajectory of the virtual particle that advects to the calculation grid point, and the field value acquisition unit 230 reaches the calculation grid point from the end position of the reverse tracking while diffusing. The field value of the virtual particle is obtained as the field value at the calculation lattice point.
Thereby, in the analysis apparatus 200, since it is not necessary to interpolate the field value of the virtual particle into the calculation lattice point, numerical diffusion is less likely to occur, and numerical analysis can be performed with higher accuracy.

Further, the field value calculation unit 231 integrates the diffusion equation along the trajectory of the virtual particle for each calculation grid point, and calculates the field value at each calculation grid point.
The field value calculation unit 231 can more accurately evaluate the diffusion term in that the diffusion of the virtual particle trajectory is calculated according to the field value around the trajectory.
Further, the convergence determination unit 232 determines whether or not the field value calculated by the field value calculation unit 231 has converged, and the field value calculation unit 231 calculates the field value until it is determined that the field value has converged. repeat.
Thus, the field value calculation unit 231 can evaluate the diffusion term based on the field value that sufficiently reflects the influence of diffusion, and can more accurately evaluate the diffusion term in this respect.

As in the case of the analysis apparatus 100, the analysis apparatus 200 can analyze both a stationary scalar field and an unsteady scalar field. When analyzing an unsteady scalar field, the analysis apparatus 200 obtains the distribution of the scalar quantity at the time t ₀ as described with reference to the equation (4).

<Third Embodiment>
FIG. 7 is a schematic block diagram showing a functional configuration of the analysis apparatus according to the third embodiment of the present invention. In the figure, the analysis apparatus 300 includes a precondition acquisition unit 110, a reverse tracking unit 220, a field value acquisition unit 330, and an analysis result output unit 140. The field value acquisition unit 330 includes an equation acquisition unit 331 and an analysis execution unit 332. In FIG. 7, parts having the same functions corresponding to the parts in FIG. 4 are denoted by the same reference numerals (110, 220, 140), and description thereof is omitted.

  Similarly to the analysis devices 100 and 200, the analysis device 300 performs numerical analysis to obtain a field value in the advection diffusion problem. However, the specific processing performed by the analysis device 300 is different from the processing performed by the analysis devices 100 and 200. The analysis apparatus 300 is constructed by a computer, for example. Alternatively, the analysis apparatus 300 may be constructed by a circuit other than a computer, such as a dedicated circuit.

Similarly to the field value acquisition units 130 and 230, the field value acquisition unit 330 performs inverse diffusion while diffusing based on the field value at the end position of reverse tracking and the result of reverse tracking for each of the calculation lattice points. The field value of the virtual particle when it reaches the calculation grid point from the tracking end position is acquired as the field value at the calculation grid point. More specifically, the field value acquisition unit 330 acquires the field value at each of the calculation lattice points from simultaneous equations obtained based on the locus of the virtual particles and the diffusion equation for each of the calculation lattice points. .
The process performed by the field value acquisition unit 330 corresponds to an example of a process in the field value acquisition step.

The equation acquisition unit 331 acquires simultaneous equations based on the virtual particle trajectory and the diffusion equation for each of the calculation lattice points.
The analysis execution unit 332 acquires the field value at each of the calculation lattice points from the simultaneous equations acquired by the equation acquisition unit 331.

Next, the operation of the analysis apparatus 300 will be described with reference to FIG.
Hereinafter, as in the case of the second embodiment, φ _k (t) indicates the scalar amount of the calculation lattice point k, and φ (tilde) _k (t) indicates the scalar amount of the virtual particle. Other variables used without showing the definition are the same as those in the first embodiment or the second embodiment.

Below, the case where the analysis apparatus 300 performs analysis of an unsteady field will be described first, and then the case where the analysis apparatus 300 performs analysis of a steady field will be described.
Also in the third embodiment, as in the case of the second embodiment, the velocity field and the scalar field are formed in the space to be analyzed, and the scalar quantity constituting the scalar field is advected on the velocity field. Consider the problem of finding the distribution of scalar quantities when diffusing. However, here, when the scalar field is non-stationary, the distribution of the scalar quantity at a new time is obtained.

Also in the third embodiment, as in the case of the second embodiment, a calculation grid having the total number of calculation grid points N is set in the space to be analyzed, and the scalar quantity at each of the calculation grid points is obtained. The calculation grid points are identified by the number k (0 ≦ k ≦ N−1).
Also in the third embodiment, as in the case of the second embodiment, a scalar field is obtained in consideration of the diffusion of virtual particles on the trajectory.
On the other hand, in the third embodiment, a processing loop for updating the scalar field by solving Equation (3) for each calculation grid point is unnecessary, and the processing time is reduced.

FIG. 8 is a flowchart showing a procedure of processing performed by the analysis apparatus 300. Steps S301 and S302 in FIG. 8 are the same as steps S201 and S202 in FIG. Also in the third embodiment, as in the case of the second embodiment, the reverse tracking unit 220 calculates the trajectory x _k (t) of the virtual particle that advects to the calculation lattice point k on the velocity field by the equation (1). ) Is solved in the opposite time direction. Also in the third embodiment, the reverse tracking unit 220 only needs to reverse-track virtual particles over a relatively long time (distance). Therefore, the analysis apparatus 300 can also analyze a field in a system (closed system) without an inflow boundary or a field in a dead water region.

  As a method of solving the equation (1), it is conceivable to solve the equation discretely. At this time, various schemes such as Euler explicit method, complete implicit method, Crank Nicholson method, Runge-Kutta method, Runge-Kutta-Ferberg method or Adams Bash force method can be used as the discretization scheme performed by the inverse tracking unit 220.

Further, when the speed u is defined discretely, it is conceivable to use an interpolation value of the speed u from a nearby speed definition point. At that time, various schemes such as primary accuracy interpolation, bilinear interpolation, trilinear interpolation, secondary accuracy interpolation, or cubic accuracy interpolation can be used as the interpolation scheme performed by the inverse tracking unit 220.
The reverse tracking unit 220 performs the reverse tracking described above for each calculation grid point.

  In addition, when the locus | trajectory obtained by a calculation error turns into an unrealistic locus | trajectory, you may make it the reverse tracking part 220 take a means (for example, correction | amendment of a speed field) which prevents it. For example, when the virtual particle reaches the boundary (wall surface) where the velocity condition is non-slip due to the influence of numerical error, the reverse tracking unit 220 corrects the velocity used for the reverse tracking to avoid stopping the virtual particle. You may make it do.

  After step S302, the equation acquisition unit 331 acquires the simultaneous equations by calculating a coefficient matrix or the like indicating the simultaneous equations related to the scalar amount of each calculation grid point (step S303). Then, the analysis execution unit 332 numerically solves the simultaneous equations acquired by the equation acquisition unit 331, thereby executing analysis and obtaining a scalar field (step S304). Hereinafter, the processing performed by the equation acquisition unit 331 and the analysis execution unit 332 will be described in more detail.

  Also in the third embodiment, as in the case of the second embodiment, the diffusion equation shown in Expression (3) is solved discretely. However, in the second embodiment, the diffusion equation is integrated along the trajectory of the virtual particle, whereas in the third embodiment, the coefficient matrix indicating the simultaneous equations related to the scalar quantities of the respective calculation lattice points from the diffusion equation. Etc. and solve the simultaneous equations. First, the simultaneous equations acquired by the equation acquisition unit 331 will be described.

  As described in the second embodiment, when integrating Equation (3), as a discretization scheme related to time, Euler explicit method, complete implicit method, Crank Nicholson method, Runge-Kutta method, Runge-Kutta-Ferberg method or Adams Bash force method Various schemes can be used. For example, when the complete implicit method is applied, the equation (6) is obtained.

However, δt indicates a time width in discretization with respect to time. Φ indicates a scalar field.
Similarly to the case of the second embodiment, for the Laplacian, the value on the locus is used as the scalar quantity at the center position, and the scalar quantity at the calculation grid point is used as the nearby scalar distribution. As a Laplacian discretization scheme, various schemes such as various Laplacian models, a second-order accurate central difference method, and a fourth-order accurate central difference method can be used. For example, when the secondary accuracy central difference method is applied, the equation (5) is obtained.

  Here, it is assumed that an arbitrary interpolated value φ (bar) (t, x) is given by a linear combination of scalar quantities φ on calculation grid points as shown in equation (7).

However, w _n (x) indicates a weighting coefficient set for each calculation grid point n according to the position x.
Then, equation (5) is discretized as shown in equation (8).

However, W _n ^{k, t} is shown in Equation (9).

  Substituting equation (8) into equation (6) and rearranging results in equation (10).

However, a _{k, t} is expressed as in equation (11).

Further, W ′ _n ^{k, t} is expressed as in Expression (12).

Also, δs _{k, t} is expressed as in equation (13).

When the time at which the virtual particle reaches the calculation lattice point k is t ₀ , Equation (14) is established.

Further, when t → t ₀ and δt → t ₀ -t ₁ (t ₀ -t ₁ > 0) are substituted into the expression (10), the expression (15) is obtained from the expression (14).

Further, when the time t ₂ (t ₂ <t ₁ ) is introduced and the same development is applied to φ (tilde) _k (t ₁ ), the following expression (16) is obtained.

A similar development is repeated until time t _M when the virtual particle leaves the reverse tracking end position to obtain Equation (17).

When the value of φ (tilde) _k (t _M ) is unknown, it is conceivable to approximate equation (17) as equation (18).

Assuming that the value of equation (19) is sufficiently small, the approximation of equation (18) is reasonable. The larger the value of t ₀ -t _M , that is, the longer the reverse tracking time is (ie, the longer the reverse tracking distance), the smaller the value of Equation (19).

On the other hand, if the value of φ (tilde) _k (t _M ) is known from the boundary condition or the like, the value may be simply substituted into equation (17).
When Equation (17) or Equation (18) is considered for all calculation lattice points k (0 ≦ k ≦ N−1), an N- _ary simultaneous linear equation relating to the scalar quantity φ _k (t) of the calculation lattice points can be obtained.
Further, a vector Φ having a scalar quantity φ _k (t ₀ ) of the calculation grid point as a component is defined as in Expression (20).

However, T shows transposition of a vector or a matrix.
The N-ary simultaneous linear equations obtained based on the equation (17) or the equation (18) can be expressed as the equation (21) using Φ.

  However, A is an N-order square matrix having the component shown in Expression (22).

B represents an N-dimensional vector.
If the value of φ (tilde) _k (t _M ) is known, and therefore based on equation (17), the component of b is as in equation (23).

On the other hand, if the value of φ (tilde) _k (t _M ) is unknown, and therefore based on equation (18), the component of b is as in equation (24).

  The equation acquisition unit 331 acquires the simultaneous equations shown in the equation (21) by obtaining the coefficient matrix A and the vector b in the equation (21). Specifically, the equation acquisition unit 331 calculates each element of the coefficient matrix A and the vector b based on the equations (22) to (24).

If Φ (t): [t _M , t ₁ ] is known, it is new (still calculated) using a known solution such as the SOR method (Successive Over-relaxation Method) or the BiCGStab method (Biconjugate Gradient Stabilized Method). No) The scalar quantity Φ of all the calculation lattice points at time t ₀ is obtained. However, (PHI) (t) shows the scalar field in the time t.
Therefore, the analysis execution unit 332 determines a scalar field (scalar amounts of all calculation lattice points) at the new time t ₀ based on the simultaneous equations obtained by the equation acquisition unit 331 and the acquired scalar amount Φ (t). Ask for.

After step S304, the process proceeds to step S305. Step S305 is the same as step S205 in FIG.
Thereafter, the process of FIG.

Next, the case where the scalar field is stationary will be described.
As in the case of the second embodiment, the velocity field and the scalar field are formed in the space to be analyzed, and when the scalar quantity constituting the scalar field diffuses while advancing on the velocity field, the steady state Consider the problem of finding the distribution of scalar quantities.
In the following, the scalar quantity of the calculation grid point _k is expressed as φk.

As in the case where the scalar field is non-stationary, a calculation grid having the total number of calculation grid points N is set in the space to be analyzed, and the scalar quantity at each calculation grid point is obtained. The calculation grid points are identified by the number k (0 ≦ k ≦ N−1).
In addition, a scalar field is calculated in consideration of the diffusion of virtual particles on the trajectory, and a processing loop for updating the scalar field by solving equation (3) for each calculation grid point is unnecessary, thereby reducing processing time. This is the same as when the scalar field is non-stationary.

Similarly to the case where the scalar field is unsteady, the inverse tracking unit 220 reverses the trajectory x _k (t) of the virtual particle that advects the velocity field to the calculation lattice point k and reverses the equation (1) in the time direction. Find by solving. Similar to the case where the scalar field is non-stationary, the reverse tracking unit 220 may perform reverse tracking of the virtual particles over a certain length of time (distance). Therefore, the analysis apparatus 300 can also analyze a field in a system (closed system) without an inflow boundary or a field in a dead water region.
In addition, the solution of equation (1), the point that various schemes can be used, and the means for preventing the trajectory obtained by the calculation error from becoming an unrealistic trajectory are the same as when the scalar field is unsteady. It is the same.

  In addition, the point and solution method for discretely solving the diffusion equation shown in Equation (3) and the point that various schemes can be used are the same as in the case where the scalar field is unsteady. When integrating the formula (3), when the complete implicit method is applied as a discretization scheme for time, the formula (6) is also the same as in the case where the scalar field is unsteady.

  Similarly to the non-stationary case, for the Laplacian, the value on the locus is used as the scalar quantity at the center position, and the scalar quantity at the calculation grid point is used as the scalar distribution in the vicinity. As a Laplacian discretization scheme, various schemes such as various Laplacian models, a second-order accurate central difference method, and a fourth-order accurate central difference method can be used. For example, when the secondary accuracy central difference method is applied, the equation (25) is obtained.

In Expression (25), since a scalar field in a steady state is handled, time t is not included in the direct argument of φ (bar). Other than that is the same as Formula (5).
Here, it is assumed that an arbitrary interpolated value φ (bar) (x) is given by a linear combination of the scalar quantity φ on the calculation grid point as shown in the equation (26).

  Then, equation (25) is discretized as shown in equation (27).

However, W _n ^{k, t} is expressed by Equation (9) as in the case of non-stationary cases.
Substituting equation (27) into equation (6) and rearranging results in equation (28).

However, as in the case of the unsteady state, a _{k, t} , W ′ _n ^{k, t} , and δs _{k, t} are respectively expressed as equations (11), (12), and (13).
If the time at which the virtual particle reaches the calculation lattice point k is t ₀ , Equation (29) is established.

Substituting t → t ₀ and δt → t ₀ -t ₁ (t ₀ -t ₁ > 0) into equation (28) yields equation (30).

Furthermore, when the time t ₂ (t ₂ <t ₁ ) is introduced and the same expansion is applied to φ (tilde) _k (t ₁ ), the equation (31) is obtained.

Similar expansion is repeated until time t _M to obtain equation (32).

When the value of φ (tilde) _k (t _M ) is unknown, it is conceivable to approximate equation (32) as equation (33).

As in the non-stationary case, assuming that the value of Equation (19) is sufficiently small, the approximation of Equation (33) is reasonable.
On the other hand, if the value of φ (tilde) _k (t _M ) is already known due to boundary conditions or the like, the value may be simply substituted into equation (32).
When Equation (32) or Equation (33) is considered for all calculation lattice points k (kε [0, N−1]), an N-ary simultaneous linear equation relating to the scalar quantity φ _k of the calculation lattice points can be obtained.
Further, a vector Φ having a scalar quantity φ _k of the calculation grid point as a component is defined as in Expression (34).

However, T shows transposition of a vector or a matrix.
The N-ary simultaneous linear equations obtained based on the equation (32) or the equation (33) can be expressed as the equation (35) using Φ.

Here, A represents a coefficient matrix of an N-order square matrix, and b represents an N-dimensional vector.
If the value of φ (tilde) _k (t _M ) is known, and therefore based on equation (32), the component of A is as in equation (36).

  Further, when based on Expression (32), the component of b is as shown in Expression (37).

On the other hand, if the value of φ (tilde) _k (t _M ) is unknown, and therefore based on equation (33), the component of A is as in equation (38).

  Further, when based on Expression (33), the component of b is as shown in Expression (39).

  The equation acquisition unit 331 acquires the simultaneous equations shown in Expression (35) by obtaining the coefficient matrix A and the vector b in Expression (35). Specifically, the equation acquisition unit 331 obtains each element of the coefficient matrix A and the vector b based on Expression (36), Expression (23), Expression (38), and Expression (24).

A known solution such as the SOR method or the BiCGStab method is used for the equation (35), and the scalar quantity Φ of all the calculation lattice points is obtained.
Therefore, the analysis execution unit 332 solves the simultaneous equations obtained by the equation acquisition unit 331 and obtains a scalar field (scalar amounts of all calculation lattice points).

As described above, the reverse tracking unit 220 reversely tracks the trajectory of the virtual particles that advect to the calculation grid point, and the field value acquisition unit 330 performs the diffusion when reaching the calculation grid point from the end position of the reverse tracking. The field value of the virtual particle is obtained as the field value at the calculation lattice point.
Thereby, in the analysis apparatus 300, since it is not necessary to interpolate the field value which a virtual particle has in a calculation lattice point, it is hard to produce numerical diffusion and can perform a numerical analysis with higher precision.

The field value acquisition unit 330 acquires the field value at each of the calculation grid points from the simultaneous equations obtained based on the locus of the virtual particles and the diffusion equation for each of the calculation grid points.
Thereby, the field value acquisition unit 330 does not need to repeat updating of the field value. In this respect, the processing time and the storage capacity necessary for the processing can be reduced.

  The inventor actually obtained the field value by the same numerical analysis as in the above embodiment. As a result, the effectiveness of the present invention could be confirmed. Hereinafter, this point will be described.

<First analysis example>
In the first analysis example, the mixing of two liquids in the case where the concentration field is in a steady state in the T-type mixer shown in FIG. 9 was analyzed by the same numerical analysis as in the first embodiment.
FIG. 9 shows an example of a typical T-type mixer used for fluid mixing. In the figure, each dimension is made dimensionless by the ratio to the flow path height h.
The fluid that has entered from the two inlets A31 or A32 passes through the inflow channel and joins in the mixing channel. The two joined fluids are mixed in the mixing channel and discharged from the discharge port A33.
Below, the case where a liquid is mixed is demonstrated to an example. For the incoming two liquids, Newtonian fluid, incompressibility, and two liquid material identity are assumed.
In this case, there are three dimensionless governing equations: a continuity equation, a Navier-Stokes equation, and an advection-diffusion equation relating to concentration, and the continuity equation is expressed as equation (40).

However, u shows speed. In addition, ▽ indicates a Laplacian operator at ▽ ² .
Also, the Navier-Stokes equation is shown as in equation (41).

  However, t shows time and p shows a pressure. Re represents the Reynolds number. The Reynolds number Re is a dimensionless number defined by the equation (42).

Where U represents the average inlet flow velocity, ν represents the kinematic viscosity coefficient, and D represents the diffusion coefficient.
Moreover, the advection diffusion equation regarding the concentration is expressed as shown in Equation (43).

  However, c shows a density | concentration and Sc shows a Schmidt number. The Schmitt number Sc is a dimensionless number defined by the equation (44).

In this analysis, it is assumed that the two liquids are fluids having the same physical property values, and the governing equations of the velocity field (continuous equation and Navier-Stokes equation) do not include concentration. For this reason, the velocity field can be solved separately from the concentration field.
Therefore, the velocity field was analyzed before the concentration field was analyzed. For the analysis of the velocity field, a grid method which is a known analysis method was used, and the total number of calculation grid points was about 2.4 million. The Reynolds number Re was set to 150.
In the analysis of the velocity field, a large lattice dependence on the velocity field was not confirmed. The flow was in a steady state.

In the analysis of the concentration field, a calculation grid different from the calculation grid used for the velocity field analysis was set. The calculation grid for the concentration field was set at the cross section of the mixing channel at the position of y = 7.67 in the y coordinate shown in FIG.
The trajectory x _k (t) of the virtual particle that reaches each calculation grid point is _obtained by converting the equation (1), which is a relational expression related to the advection of the virtual particle, into the equation (45) obtained by discretizing the Euler explicit method in the direction going back in time. Calculate by solving numerically.

The process for solving Expression (45) corresponds to an example of the process for solving Expression (1).
FIG. 10 is an explanatory diagram illustrating an example of reverse tracking of virtual particles in the concentration field analysis. In the figure, initial planes A41 and A42 are set in two inflow channels, respectively. The initial planes A41 and A42 correspond to an example of a region having a known scalar value, and the density c values are 0 and 1, respectively. In the mixing channel, a target plane A43 on which a calculation grid is set is shown.
In the calculation of the trajectory, the trajectory L41 is obtained from the reverse tracking start position P42, which is one of the calculation grid points on the target plane A43, in the direction of going back in time as indicated by the arrow V41 until the initial plane A41 or A42 is reached. .

In obtaining the trajectory, trilinear interpolation from the speed definition point (the position at which the speed was obtained by analyzing the speed field) was applied as an interpolation scheme for the speed u. Further, the Courant number was set to 0.1, and the time width δt was determined so as to satisfy this condition.
Since the obtained trajectory includes the influence of a numerical error, it may reach a boundary (wall surface) where the speed condition is non-slip. As a countermeasure in this case, the velocity used for reverse tracking was corrected so that the velocity in the direction perpendicular to the wall surface becomes zero in the time step before the arrival at the wall surface was detected.

  In reverse tracking, the virtual particles eventually go to one of the two inlets. At that time, among the coordinate values of the positions of the virtual particles, the initial concentration value at the time of inflow is determined based on the value of the x coordinate shown in FIG. Specifically, when x ≦ −2, it is determined that the density c = 0 (initial plane A41 side). On the other hand, when x ≧ 2, it is determined that the density c = 1 (initial plane A42 side).

  By performing backward tracking for all calculation grid points, the scalar value of each calculation grid point is set to 0 or 1. The scalar distribution reflecting the diffusion effect was obtained based on the equation (46) modeling the diffusion equation with the scalar distribution thus obtained as the initial value (ξ = 0).

However, (tau) shows the residence time of a virtual particle. Ξ is a temporarily set variable.
Specifically, the SOR method was applied to the equation (47) obtained by discretizing the equation (46) by the complete implicit method, and the scalar distribution at ξ = 1 was obtained.

The Schmitt number Sc = 3600. Equation (46) corresponds to an example of Equation (2).
FIG. 11 is an explanatory diagram showing a result of comparison of the calculation result of the degree of mixing MI in the target plane A43 (FIG. 10) with the lattice method while changing the interval between calculation lattice points. In FIG. 11, each point indicated by a circle indicates a calculation result in the analysis method of the present embodiment (hereinafter referred to as “the present method” in the present embodiment). On the other hand, each point indicated by a triangle indicates a calculation result in the lattice method.
The degree of mixing MI is defined by equation (48).

However, w shows the perpendicular | vertical speed with respect to the surface A (surface used as the object of integration in Formula (48)). Further, c (bar) indicates a concentration when completely mixed (that is, c (bar) = 0.5).
A commercially available CFD (Computational Fluid Dynamics) code was used for the analysis by the lattice method, and a code generated by the inventor in C ++ was used for the analysis by this method. In addition, the maximum number of calculation lattice points in the lattice method was about 66 million, and the spatial resolution at that time was 120 in points per unit length in one direction. On the other hand, the maximum number of calculation lattice points in this method was about 520,000, and the spatial resolution at that time was 512 in points per unit length in one direction.

FIG. 12 is an explanatory diagram illustrating an example of the analysis result of the density distribution at the maximum spatial resolution between the lattice method and the present method. The upper diagram F11 shows the concentration distribution by the lattice method, and the lower diagram F12 shows the concentration distribution by this method.
Comparing FIGS. F11 and F12, the degree of mixing is greater in FIG. F11. This is thought to be due to overestimation of diffusion due to numerical diffusion in the lattice method.

In the calculation grid used in the grid method, the maximum number of calculation grid points is so large that it can be said to be the limit of a general computer. Nevertheless, the lattice method causes overestimation of diffusion due to numerical diffusion, and the dependence on the computational lattice still remains strong. As a result, a correct analysis result is not obtained by the lattice method.
On the other hand, in this method, since the scalar field of the virtual particle reaching the calculation lattice point is obtained, it is not necessary to interpolate the scalar amount of the virtual particle to the scalar amount of the calculation lattice point. Further, since the scalar quantity is known at the reverse tracking end position, it is not necessary to interpolate the scalar quantity at the calculation grid point with the scalar quantity of the virtual particle. In these respects, the present technique has a feature that it is difficult to cause numerical diffusion fundamentally. As a result, this method shows an overwhelmingly small calculation grid dependency.

  The required analysis time at the maximum spatial resolution was about 10 hours in the lattice method, but about 30 minutes in this method, and the analysis could be performed in about 1/20 time. Also, the storage capacity used was about 45 gigabytes (GB) by the lattice method, but in this method, it was less than 1 gigabyte and analysis was possible with a storage capacity of about 1/45.

<Second analysis example>
In the second analysis example, the mixing of the two liquids when the concentration field is in a steady state in the T-type mixer shown in FIG. 9 was analyzed by the same numerical analysis as in the third embodiment.
Similar to the first analysis example, the case of mixing liquids will be described as an example. For the incoming two liquids, Newtonian fluid, incompressibility, and two liquid material identity are assumed.
Further, as in the first analysis example, the governing equation has three continuations: equation (40), Navier-Stokes equation (41), and advection-diffusion equation relating to concentration shown in equation (43). It is a formula.

  In addition, as in the first analysis example, the velocity field was analyzed before the concentration field was analyzed. As in the first analysis example, for the analysis of the velocity field, a lattice method which is a known analysis method was used, the total number of calculation lattice points was about 2.4 million, and the Reynolds number Re = 150.

  On the other hand, unlike the first analysis example, in this analysis example, the calculation grid for the concentration field is set in the same region as the calculation grid for the velocity field. However, the spatial resolution was changed as compared with the calculation grid for the velocity field, and the calculation of the region existing far from the position of interest was omitted as appropriate.

The trajectory x _k (t) of the virtual particle that reaches each calculation lattice point with respect to the concentration field is similar to the first analysis example. It is obtained by numerically solving Expression (45) discretized in an explicit manner.
In obtaining the trajectory, trilinear interpolation from the speed definition point was applied as an interpolation scheme for the speed u. Further, the Courant number was set to 0.1 and the diffusion number was set to 0.075, and the time width δt was determined so as to satisfy this condition.

As in the first analysis example, as a countermeasure when virtual particles reach the boundary (wall surface) where the speed condition is non-slip, the velocity in the direction perpendicular to the wall surface at the time step before the arrival at the wall surface is detected. The speed used for reverse tracking was corrected so that becomes zero.
Further, the reverse tracking termination condition is that the virtual particles reach a position where the concentration is known or satisfy the formula (49).

Here, the smaller the allowable value ε (ε ≧ 0) is set, the higher the analysis accuracy becomes, but the cost for analysis (analysis time, calculation load, etc.) increases. In this analysis, ε = 0.1 is empirically set.
Analysis was performed under the above settings to obtain a virtual particle trajectory x _k (t): [t _M , t ₀ ]. However, k is a positive integer of 0 ≦ k ≦ N−1, where N is the total number of calculation lattice points. Further, t _M ≦ t ₀ .

In the analysis of the concentration field, the concentration field was determined by numerically solving the N-ary simultaneous linear equation shown in Equation (35). At this time, as shown in the equation (50), as specific examples of the scalar quantities φ ₀ , φ ₁ ,..., Φ _N−1 in the equation (34), the concentrations c ₀ , c _{1 ,. -1} is used.

When the virtual particle concentration c (tilde) _k (t _M ) at the position x _k (t _M ) is unknown, the coefficient matrix A is set as shown in the equation (38).
However, a _{k, tm} is expressed as in equation (51).

Further, W ′ _n ^{k, tm} is expressed as in Expression (52).

W _n ^{k, tm} is expressed as shown in Expression (53).

W _n (x) is a weighting coefficient that gives an interpolation value c (bar) (x) from the scalar distribution on the calculation grid point to the position x. The interpolated value c (bar) (x) is expressed as in equation (54).

Further, when the virtual particle concentration c (tilde) _k (t _M ) at the position x _k (t _M ) is unknown, the vector b is expressed by the equation (55).

On the other hand, if the concentration of the virtual particles c _(tilde) k _{(t M)} is known at position _x k _{(t M),} the coefficient matrix A, and as shown in equation (36).
However, as with the concentration c of the virtual particles _(tilde) k _{(t M)} is unknown at the position _{x k (t M), a} k, tm, W 'n k, tm, W n k, tm , respectively , (50), (51), and (52). As _for the _w n (x), is the same as that concentration of the virtual particles c _(tilde) k _{(t M)} is unknown at the position _x k _{(t M).}
Further, when the virtual particle concentration c (tilde) _k (t _M ) at the position x _k (t _M ) is known, the vector b is expressed by the equation (56).

In addition, Formula (51), Formula (52), Formula (53), Formula (54), Formula (55), and Formula (56) are respectively Formula (11), Formula (12), Formula (9), It corresponds to an example of Expression (26), Expression (24), and Expression (23).
The BiCGStab method was used as a method for solving the N-ary simultaneous linear equations shown in Equation (35). The Schmitt number Sc = 10. In the second analysis example, an example is shown in which the time term of the advection diffusion equation of the concentration field is discretized in a completely implicit manner and the Laplacian is discretized in a second-order central difference method.

  FIG. 13 is an explanatory diagram showing an example of the calculation result of the degree of mixing MI in the target plane A43 (FIG. 10) compared with the lattice method while changing the interval between calculation lattice points. In FIG. 13, each point indicated by a black circle indicates a calculation result in the analysis method of the third embodiment (hereinafter referred to as “this method” in the second analysis example). On the other hand, each point indicated by a white circle indicates a calculation result in the lattice method. The definition of the mixing degree MI is shown in the equation (48).

A commercially available CFD code was used for analysis by the lattice method, and a code generated by the inventor in C ++ was used for analysis by this method.
In addition, the maximum number of lattice points calculated in the lattice method was about 76 million, and the spatial resolution at that time was 126 points per unit length in one direction. In the calculation grid used in the grid method, the maximum number of calculation grid points is so large that it can be said to be the limit of a general computer. Nevertheless, the lattice method causes overestimation of diffusion due to numerical diffusion, and the dependence on the computational lattice still remains strong. As a result, a correct analysis result is not obtained by the lattice method.

On the other hand, in this method, since the scalar field of the virtual particle reaching the calculation lattice point is obtained, it is not necessary to interpolate the scalar amount of the virtual particle to the scalar amount of the calculation lattice point. In addition, when the scalar quantity is known at the reverse tracking end position, it is not necessary to interpolate the scalar quantity of the calculation grid point with the scalar quantity of the virtual particle. In these respects, the present technique has a feature that it is difficult to cause numerical diffusion fundamentally. As a result, this method shows an overwhelmingly small calculation grid dependency, and the analysis has converged. That is, even if the calculation grid point interval is made smaller (shorter), the solution does not change much.
Note that the same analysis was performed in the analysis method in the second embodiment and compared with this method. As a result, in this method, it was possible to perform the analysis with an analysis time of 1/10 or less and a storage capacity of about 1/10 compared with the analysis method in the second embodiment.

FIG. 14 is an explanatory diagram showing an example of the analysis result of the concentration distribution between the lattice method and the present method. The upper diagram F21 shows the concentration distribution obtained at the maximum spatial resolution of the lattice method, and the lower diagram F22 shows the concentration distribution obtained by this method.
Comparing FIGS. F21 and F22, the density distribution of the same level as that of FIG.

  Whereas the required analysis time at the maximum spatial resolution of the grid method is about 10 hours, in this method, the spatial resolution (about 260,000 calculation grid points, 32 per unit length) gives the same degree of accuracy as the degree of mixing. Point) and about 5 minutes. Therefore, with this method, the same level of analysis results could be obtained in about 120 times less than the lattice method. Also, the storage capacity used was about 50 gigabytes in the lattice method, whereas it was less than 2 gigabytes in this method. Therefore, in this method, the same analysis result can be obtained with a storage capacity of about 1/25 of the lattice method.

  FIG. 15 is an explanatory diagram showing a comparative example of the analysis result in the first analysis example, the analysis result in the second analysis example, and the analysis result by the lattice method. Of each figure shown in the figure, FIGS. F311 to F314 show examples of analysis results by the lattice method. Further, FIGS. F321 to F324 show examples of analysis results in the second analysis example. Further, FIGS. F331 to F334 show examples of analysis results in the first analysis example.

  Also, the degree of diffusion is greater in the diagram on the left side. Specifically, in FIGS. F311, F321, and F331, the Schmitt number Sc = 1 is set. Also, Sc = 10 is set in FIGS. F312, F322, and F332. Also, Sc = 100 is set in FIGS. F313, F323, and F333. Also, Sc = 1000 is set in FIGS. F314, F324, and F334.

  FIG. F313 and FIG. F314 have the same analysis results, although the setting of the degree of diffusion is different. This is thought to be because the influence of numerical diffusion in the lattice method is greater than the change in the degree of diffusion. That is, in FIG. F314, the lattice method cannot be analyzed with high accuracy due to the influence of numerical diffusion. On the other hand, in FIGS. F323, F324, F333, and F334, analysis results corresponding to the degree of diffusion are obtained. Thus, in the analysis method in the first embodiment (analysis result in the first analysis example) and the analysis method in the third embodiment (analysis result in the second analysis example), the influence of numerical diffusion is hardly affected. High-precision analysis can be performed without receiving it.

In addition, in FIGS. F321 to F324, good analysis results are obtained. Thus, in the analysis method in the third embodiment, a good analysis result can be obtained regardless of whether the degree of diffusion is large or small.
On the other hand, with respect to FIGS. F331 to F334, the analysis accuracy is slightly lowered in FIG. F331 where the degree of diffusion is large, but good analysis results are obtained in FIGS. F333 and F334 where the degree of diffusion is small. In the analysis method according to the first embodiment, when the degree of diffusion is relatively small, a highly accurate analysis result can be obtained with a short processing time and a small storage capacity.

Note that a program for realizing all or part of the functions of the analysis apparatus 100, 200, or 300 is recorded on a computer-readable recording medium, and the program recorded on the recording medium is read into a computer system and executed. By doing so, the processing of each unit may be performed. Here, the “computer system” includes an OS and hardware such as peripheral devices.
Further, the “computer system” includes a homepage providing environment (or display environment) if a WWW system is used.
The “computer-readable recording medium” refers to a storage device such as a flexible medium, 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 embodiment of the present invention has been described in detail with reference to the drawings. However, the specific configuration is not limited to this embodiment, and includes design changes and the like without departing from the gist of the present invention.

100, 200, 300 Analysis device 110 Precondition acquisition unit 120, 220 Reverse tracking unit 130, 230, 330 Field value acquisition unit 131 Field value setting unit 132 Field value update unit 140 Analysis result output unit 231 Field value calculation Unit 232 convergence determination unit 331 equation acquisition unit 332 analysis execution unit

Claims (6)

For each calculation grid point in the calculation grid set as the field to be analyzed, reverse tracking is performed by tracing the trajectory of a virtual particle having a variable field value and advancing to the calculation grid point back in time. A reverse tracking step;
For each of the calculation grid points, when the calculation grid point is reached from the reverse tracking end position while diffusing, based on the field value at the reverse tracking end position and the result of the reverse tracking, A field value acquisition step of acquiring a field value of the virtual particle as a field value at the calculation lattice point;
The analysis method characterized by comprising. In the reverse tracking step, the trajectory of the virtual particles is obtained for each of the calculation lattice points,
In the field value acquisition step, a field value at each of the calculation grid points is acquired from simultaneous equations obtained based on a locus of the virtual particles and a diffusion equation for each of the calculation grid points.
The analysis method according to claim 1. In the reverse tracking step, the trajectory of the virtual particles is obtained for each of the calculation lattice points,
The field value acquisition step includes:
For each of the calculation grid points, integrating a diffusion equation along the trajectory of the virtual particles, and calculating a field value at each of the calculation grid points,
A convergence determination step for determining whether or not the field value calculated in the field value calculation step has converged;
Comprising
The process in the field value calculation step is repeated until the convergence value determination step determines that the field value calculated in the field value calculation step has converged. Analysis method. The field value acquisition step includes:
For each of the calculation grid points, a field value setting step of setting, in the calculation grid point, a field value that the virtual particle that reaches the calculation grid point has at the end position of the reverse tracking;
For each of the calculation grid points, a field value update step of updating the field value at the calculation grid point according to the residence time required for the virtual particle to reach the calculation grid point from the end position of the reverse tracking When,
The analysis method according to claim 1, further comprising: For each calculation grid point in the calculation grid set as the field to be analyzed, reverse tracking is performed by tracing the trajectory of a virtual particle having a variable field value and advancing to the calculation grid point back in time. A reverse tracking unit;
For each of the calculation grid points, based on the field value at the end position of the reverse tracking and the result of the reverse tracking, when the calculation grid point is reached from the end position of the reverse tracking while diffusing A field value acquisition unit for acquiring a field value in a virtual particle as a field value in the calculation lattice point;
An analysis apparatus comprising: On the computer,
For each calculation grid point in the calculation grid set as the field to be analyzed, reverse tracking is performed by tracing the trajectory of a virtual particle having a variable field value and advancing to the calculation grid point back in time. A reverse tracking step;
For each of the calculation grid points, when the calculation grid point is reached from the reverse tracking end position while diffusing, based on the field value at the reverse tracking end position and the result of the reverse tracking, A field value acquisition step of acquiring a field value of the virtual particle as a field value at the calculation lattice point;
A program for running