JP2012026785A - Ground deformation analysis device, ground deformation analysis method, program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a ground deformation analysis device and the like capable of performing a three dimensional large deformation analysis and dynamic analysis of ground.SOLUTION: The ground deformation analysis device 1 for analyzing a ground deformation using a particle method includes: a storage means for storing element information on an element region having a solid shape and particle information on a particle; an initial element region calculation means for calculating an initial element region where the particle locates; an initial stress calculation means for calculating the initial stress of the particle; a node velocity calculation means for calculating a first velocity of a node; a position/strain increment calculation means for calculating the position and the strain increment of the particle; a strain/stress/density calculation means for calculating the strain, stress, and density of the particle; and an element region calculation means for calculating the element region where the particle locates, and if the particle is outside of the element region, translating the particle position to a predetermined position on the outer surface of the element region according to the particle position.

Description

  The present invention relates to a ground deformation analysis apparatus, a ground deformation analysis method, and a program for analyzing ground deformation.

  For example, during an earthquake or the like, disasters accompanied by large ground deformation such as slope failure and liquefaction frequently occur due to dynamic action on the ground. At that time, in order to ensure the safety of the surrounding environment, it is important to predict the degree of damage caused by those disasters. For these damage predictions, ground deformation analysis is indispensable, and it is desirable that three-dimensional deformation analysis and dynamic analysis are possible.

  In the case of such a ground deformation analysis, a deformation analysis by a finite element method has been conventionally performed. Such an example is described in Patent Document 1, and the action of water that resists the volume deformation of the constitutive law of soil is simply modeled, and an earthquake response analysis is performed in consideration of the effect of liquefaction on sand ground and the like. It is described.

JP 2006-038455 A

  However, the conventional finite element method has a problem that it cannot cope with large deformation analysis of the ground. For example, in the finite element method, the deformation of the ground is expressed by the movement of the nodes that form the element. However, due to problems such as the failure of the calculation due to the large deformation of the element, the ground is greatly displaced as described above. Difficult to use for analysis.

  On the other hand, in recent years, a method called a particle method that can cope with large deformation analysis of an object has been developed. Here, an outline of the particle method will be briefly described. FIG. 14 is a diagram schematically showing an outline of MPM (Material Point Method) which is a particle method.

  The particle method is an explicit time history analysis method. As shown in FIG. 14A, individual particles 100 (Lagrange particles) having substance information (Lagrange variable) as particle information move on or between elements 105 in the space. Material information which is a Lagrangian variable of the particle 100 is aggregated through an interpolation function or the like at a node 120 of the element 105 where the particle 100 exists at every fixed time interval. Then, the equation of motion is solved for the node 120 to determine the speed increment of the node 120 in the next step.

  At this time, as shown in FIG. 14B, the element 105 is deformed while the particle 100 is placed by the movement of the node 120, and the particle 100 moves. Further, the substance information of the particle 100 is updated from the substance information of the node 120 through an interpolation function or the like. Then, as shown in FIG. 14C, the deformed element 105 is returned to the original position again in preparation for the next step, leaving the particles 100 that have moved. By repeating the above process, the particle 100 moves its position and its substance information is updated.

  However, the particle method cannot be used as it is in consideration of the above-mentioned problems of three-dimensional deformation of the ground and dynamic analysis. In other words, it is necessary to deal with three-dimensionalization, in addition to dealing with dynamic analysis and groundwater behavior.

  The present invention has been made in view of the above-described problems, and an object of the present invention is to provide a ground deformation analysis apparatus and the like that enable three-dimensional large deformation analysis and dynamic analysis of the ground. For the MPM, which is a particle method that was devised and devised by the scientists, it was not possible with the conventional FEM analysis by using a method that realized the correspondence to three-dimensionalization and the dynamic analysis and the behavior of groundwater. 3D large deformation analysis of ground is possible.

  A first invention for achieving the above-described object is a ground deformation analysis apparatus that performs a ground deformation analysis using a particle method, and includes at least identification information of the element region as element information relating to an element region having a three-dimensional shape. Storage means for storing the position of the node forming the element region, and storing at least particle identification information, identification information of the element region in which the particle is located, position, mass, density, and stress as particle information about the particle; Based on the position of the particle and element information, the particle is located adjacent to the nodal point, with the nodal point forming the one surface as the origin and the one side forming the one surface with respect to the one side of the element region having the three-dimensional shape. Determine whether the particle is located in the direction indicated by the outer product of two vectors toward each of the two nodes, and the particle is indicated by the outer product for all faces of the element region An initial element region calculating means for calculating the number of the particles in the element region, and calculating the mass of the particle using an interpolation function. The initial nodal force of the nodal point is calculated based on the initial nodal force of the nodal point based on the initial nodal force of the nodal point, the mass of the particle, and the density using a gradient function. An initial stress calculating means for calculating an initial stress of the particle based on a directional stress; and for the element region, the mass of the node is calculated based on the mass of the particle using an interpolation function, and an interpolation function and The nodal force of the node is calculated based on the mass, density, and stress of the particle using a gradient function, the acceleration of the nodal point is calculated based on the mass of the nodal point and the nodal force of the nodal point, and the acceleration of the nodal point Calculate the first velocity of the node based on A nodal velocity calculating means for calculating a position of the particle based on a first velocity of the nodal point using an interpolation function and an acceleration of the nodal point using an interpolation function. Position / strain increment calculating means for calculating the particle velocity and calculating a strain increment of the particle based on the velocity of the particle using a gradient function, and based on the particle strain increment in the element region Based on the strain / stress / density calculating means for calculating strain, stress and density of the particles, and the position of the particles and the element information, the bottom surface of the element region located at the lowest position corresponding to the planar position of the particles From the vertical coordinate and the vertical length of the element area, an element area where the particle is located is calculated, and when the particle is outside the element area, the position of the particle is Depending on the position, the element region It is an element region calculation means for converting into a predetermined position on the outer surface of the region.

  The ground deformation analysis apparatus according to the first aspect of the present invention may further include a water pressure calculating means for calculating a water pressure using the strain of the particles or the strain increase of the particles.

  The nodal velocity calculating means may calculate the acceleration of the nodal point using a motion equation simulating a dashpot for some of the nodal points.

  The second invention for achieving the object described above stores at least the identification information of the element area and the position of the node forming the element area as element information relating to the element area having a three-dimensional shape, and as particle information relating to particles. An information processing apparatus having storage means for storing at least particle identification information, element region identification information, position, mass, density, and stress, and performing ground deformation analysis using a particle method, Based on the position of the particle and the element information, the particle is located adjacent to the node by forming the one surface with respect to one surface of the element region having the three-dimensional shape, with the predetermined node forming the one surface as the origin. It is determined whether it is located in the direction indicated by the outer product of two vectors toward each of the two nodes, and the particles are positioned in the direction indicated by the outer product for all the faces of the element region. An initial element region calculating step of calculating the element region as an element region in which the particle is located and calculating the number of the particles in the element region; and based on the mass of the particle using an interpolation function. The initial nodal force of the particle is calculated based on the initial nodal force of the nodal point, the mass and density of the particle using a gradient function, and the initial vertical stress of the particle is calculated. An initial stress calculating step for calculating the initial stress of the particle based on the above, and for the inside of the element region, the mass of the node is calculated based on the mass of the particle using an interpolation function, and the interpolation function and the gradient function To calculate the nodal force of the node based on the mass, density, and stress of the particle, calculate the nodal acceleration based on the mass of the nodal point and the nodal force of the nodal point, and based on the acceleration of the nodal point To calculate the first velocity of the node A nodal velocity calculating step for calculating the position of the particle based on a first velocity of the nodal point using an interpolation function and an acceleration of the nodal point using an interpolation function. A position / strain increment calculating step for calculating the particle velocity and calculating a strain increment of the particle based on the particle velocity using a gradient function; and based on the particle strain increment in the element region. Based on the strain / stress / density calculating step for calculating strain, stress, and density of the particle, and the position of the particle and the element information, the bottom surface of the element region located at the lowest position corresponding to the planar position of the particle From the vertical coordinate and the vertical length of the element area, an element area where the particle is located is calculated, and when the particle is outside the element area, the position of the particle is Depending on position An element region calculating step for converting the element region into a predetermined position on the outer surface of the element region.

  In the ground deformation analysis method of the second invention, the information processing apparatus may further perform a water pressure calculation step of calculating a water pressure using the strain of the particles or the strain increment of the particles.

  In the nodal velocity calculation step, the acceleration of the nodal point may be calculated for some of the nodal points using an equation of motion simulating a dashpot.

  A third invention for achieving the above-described object is a program for causing an information processing device to function as the ground deformation analysis device of the first invention.

  With the above configuration, the ground deformation analysis is performed using the particle method. At this time, an element region of an arbitrary three-dimensional shape is determined, and the calculation of the element region where the particle is located and the initial stress are set correspondingly. . By using the particle method, large deformation analysis of the ground becomes possible, by defining the element area of any three-dimensional shape, by calculating the element area where the particle is located and setting the initial stress correspondingly, Three-dimensional deformation analysis becomes possible. In addition, by calculating water pressure based on particle strain and strain increment, it is possible to analyze groundwater. It is also possible to perform dynamic analysis using an equation of motion simulating a dashpot.

  According to the present invention, it is possible to provide a ground deformation analysis apparatus and the like that enable three-dimensional large deformation analysis and dynamic analysis of the ground.

The figure which shows the hardware constitutions of the ground deformation analysis apparatus 1 The figure shown about the particle | grains 100 and the element area | region 105 Flow chart showing the overall flow of the ground deformation analysis method by the ground deformation analysis apparatus 1 Flow chart showing the flow of the ground deformation analysis method by the ground deformation analysis apparatus 1 The figure shown about the process of the ground deformation analysis method by the ground deformation analysis apparatus 1 The figure shown about the process of the ground deformation analysis method by the ground deformation analysis apparatus 1 The figure shown about the process of the ground deformation analysis method by the ground deformation analysis apparatus 1 The figure shown about the process of the ground deformation analysis method by the ground deformation analysis apparatus 1 Flow chart showing the overall flow of the ground deformation analysis method by the ground deformation analysis apparatus 1 The figure shown about the process of the ground deformation analysis method by the ground deformation analysis apparatus 1 The figure which shows the example of the ground deformation analysis by the ground deformation analysis apparatus 1 The figure which shows the example of the ground deformation analysis by the ground deformation analysis apparatus 1 The figure which shows the example of the ground deformation analysis by the ground deformation analysis apparatus 1 Diagram showing the outline of the particle method

Hereinafter, embodiments of the ground analysis device and the like of the present invention will be described with reference to the drawings.
First, the hardware configuration of the ground deformation analysis apparatus 1 of the present embodiment will be described with reference to FIG. In the ground deformation analysis apparatus 1, for example, a control unit 11, a storage unit 13, a media input / output unit 15, a peripheral device I / F unit 17, a communication unit 21, an input unit 23, a display unit 25, and the like are connected via a bus 19. Configured.

  The control unit 11 includes a CPU, ROM, RAM, and the like. The CPU calls and executes a program stored in the storage unit 13, ROM, recording medium, or the like to a work memory area on the RAM, and drives and controls each unit connected via the bus 19. The ROM permanently holds a computer boot program, a program such as BIOS, data, and the like. The RAM temporarily holds the loaded program and data, and includes a work area used by the control unit 11 to perform various processes.

  The storage unit 13 is a hard disk drive, and stores a program executed by the control unit 11, data necessary for program execution, an operating system, and the like. These program codes are read by the control unit 11 as necessary, transferred to the RAM, and read and executed by the CPU. The storage unit 13 stores programs, data, and the like related to ground deformation analysis processing described later.

  The media input / output unit 15 (drive device) is a media input / output device such as a floppy (registered trademark) disk drive, CD drive, DVD drive, and MO drive, and performs data input / output.

  The peripheral device I / F (interface) unit 17 is a port for connecting a peripheral device, and transmits and receives data to and from the peripheral device via the peripheral device I / F unit 17. The peripheral device I / F unit 17 is configured by a USB or the like, and usually includes a plurality of peripheral devices I / F. The connection form with the peripheral device may be wired or wireless.

  The communication unit 21 includes a communication control device, a communication port, and the like, and is a communication interface that mediates communication with a network or the like, and performs communication control.

  The input unit 23 is an input device such as a keyboard, a pointing device such as a mouse, or a numeric keypad, and outputs input data to the control unit 11.

  The display unit 25 includes a display device such as a liquid crystal panel or a CRT monitor, and a logic circuit (video adapter or the like) for executing display processing in cooperation with the display device, and is input under the control of the control unit 11. Display information is displayed on a display device.

  The bus 19 is a path that mediates transmission / reception of control signals, data signals, and the like between the devices.

  As described above, during an earthquake or the like, disasters with large ground deformation such as slope failure and liquefaction frequently occur due to dynamic action on the ground. For these damage predictions, ground deformation analysis is indispensable, and it is desirable to be able to analyze three-dimensional behavior. However, it is difficult to cope with large deformation analysis of the ground by a finite element method conventionally used for ground deformation analysis.

  The ground analysis apparatus 1 of the present invention solves the above problems, adapts the MPM, which is a particle method, to large ground deformation analysis considering the three-dimensional behavior, and analyzes the three-dimensional behavior of the ground and its This is a process that can be applied to dynamic analysis and analysis of groundwater behavior. Hereinafter, the processing will be described. Each process shown below is executed by the ground deformation analysis apparatus 1 when the CPU of the control unit 11 reads a program, data, or the like stored in advance in the storage unit 13 or the like.

In the processing performed by the ground deformation analysis apparatus 1 according to the present embodiment, a three-dimensional model of the ground is obtained by the planar distribution of the earth pillars 30 including the pillars of the elements 105 formed by the nodes 120 as shown in FIG. To do. By performing the calculation process shown in the flow of FIG. 3 in the model, the three-dimensional analysis of the ground using MPM is enabled. Each particle 100 is associated with a particle number that is particle identification information as particle information 110 that is a Lagrange variable, as shown in FIG. 2 (b-1), and the element number, position, and speed at which the particle is located. , Stress, strain, density, mass, and other information, and the particle information 110 is stored in the RAM or the storage unit 13.
FIG. 3 shows an algorithm flow in the ground deformation analysis method by the ground deformation analysis apparatus 1.

  As shown in FIG. 3, first, in step S100, initial conditions related to calculation, such as input of initial values required for calculation, are set.

That is, processing before starting the time history calculation is required first. Here, the position X _p (k) = {x _p (k), y _p (k), z _p (k)} (k represents the number of time steps), velocity v _p (k) = {v _xp ( k), v _yp (k), v _zp (k)} ^T , strain ε _p (k) = ε _{ij p} (k), and initial values such as density ρ _p (k), mass M _p (k = 0) ) In association with a particle number or the like, which is identification information indicating each particle 100, and is stored in the RAM or the storage unit 13.

For mass M _p ,
There is a relationship. V _p = L ³ / n _s ³ , L ³ is the volume of the element region, and n _s is the number of particles per element side. In the MPM, the object is modeled as a particle having a constant mass, so that the mass M _p becomes an invariable value during time history calculation.

  In addition, if there are constants necessary for calculation processing described later, or values that require setting of initial values when using models with history dependence such as elasto-plastic models, those constants and initial values are also set in the same way. Set in S100 and stored in RAM or the like.

In addition, an element area of each element 105 is set. As described above, each element region is formed by the nodes 120, and each element region has a flat earth column 30 formed of columns of elements (cells) 105 constituted by the nodes 120, as shown in FIG. To be distributed.
For this reason, as shown in FIG. 2B-2, for example, the node number for identifying the node 120 forming this element region is the element number for identifying each element, and the local node number to be described later for each element region. It is memorized in association with. Each node number is stored in association with its position as shown in FIG.
The above information is stored as element information 111 in the RAM, the storage unit 13 or the like. In this embodiment, each node forms an Euler lattice having a rectangular element region in the XY plane shown in FIG. 2A, the element region of the element 105 forms an Euler lattice also in the XZ plane or the YZ plane in FIG. 2A. However, the present embodiment is not limited to this.

  Next, in step S110, various calculations relating to the position of the particles are performed. FIG. 4 shows the processing in step S110.

  As shown in FIG. 4, in step 110, first, before starting the time calendar calculation (k = 0) (Yes in step S <b> 111), in step S <b> 112, corresponding to each particle based on the particle position and element information. An element region (element number) to be calculated is calculated and stored as particle information 110. Thereby, the particle | grains located in each element area | region are pinpointed, for example, it can link with an element number and can also record a particle number.

That is, in order to define the presence of particles as an initial setting, it is necessary to extract particles in each element region. Since the element region has an arbitrary three-dimensional shape, particles are extracted by using the normal of each surface constituting the element region using the element information and the position of the particle.
Specifically, for each particle, with respect to one surface of the element region, the particle is defined as the origin of a predetermined node that forms the one surface, and the two nodes that are adjacent to the predetermined node that form the one surface. It is determined whether it is located in the direction indicated by the outer product of the two vectors toward each. When the above determination is made for all the faces of the element area and the particle is positioned in the direction indicated by the outer product with respect to all the faces, the element information of the element area (initial element area) as the particle information of the particle Write a number to identify that the particle is inside the element area. This will be described below.

  Since the element region is composed of a hexahedron, a method of extracting particles inside the hexahedron will be described with reference to FIG.

First, the local coordinates of the quadrilateral with nodes 1 to 4 that are one side of the hexahedron are calculated using the outer product. For example, the vectors m ₁ , m ₂ , and m ₃ in FIG. 5A showing the local coordinate system are respectively expressed as vectors a and b in FIG.
Is calculated.

Also, the position of the particles represented by a vector _x = shown in FIG. _{5 (b) (x 1,} x 2, x 3), and the local coordinate system by a vector _{m i,} 2 two of the global coordinate system by a vector _{e i} When displayed in the coordinate system,
It becomes. By taking an inner product vector m _i, For = vector after local coordinate changes _{X (X 1, X 2,} X 3),
Therefore, the following formula
Thus, the particle position represented by coordinates x = (x ₁ , x ₂ , x ₃ ) can be calculated as position X = (X ₁ , X ₂ , X ₃ ) in the local coordinate system. _{Here, m ij} indicates a _{e j} component of the _{m i.}

Here, as shown in FIG. 5C, it is assumed that the position of the particle p is represented by a vector S in the global coordinate system, and one surface 105b of an element region which is a hexahedron is constituted by the nodes 1 to 4. . At this time, using the vector T of the global coordinate system of the point 1 which is one of the points constituting the one surface 105b, the matrix M which is a transformation matrix expressed by the elements m _{11 to} m ₃₃ in the above equation (5). ,
When the component L ₃ corresponding to the component X ₃ in the above equation (5) is L ₃ ≧ 0, an arbitrary point represented by the vector S is shown on the plane 105b. It is on the vector m ₃ side indicated by 5 (a). In FIG. 5C, the point 1 is the origin, and the vector of the point 1 → 2 from the point 1 to the point 2 adjacent to the point 1 is adjacent to the vector a and the point 1 to the point 1 in FIG. The point 1 → 4 toward the matching point 4 corresponds to the vector b, and using these, the elements m _{11 to} m ₃₃ of the matrix M are calculated by the equation (2). That is, the vector m ₃ side is a direction indicated by the outer product of the vector of the point 1 → 2 and the vector of the point 1 → 4.

If there are particles on the m ₃ side under all the following conditions, the particles exist inside the hexahedron. Note that (1 → 2, 3, 4) means that the point 1 is the start point of the vectors a and b described above. Further, m ₁ (1 → 2) indicates that the direction of the vector m ₁ shown in FIG. 5A is the direction from point 1 to point 2 shown in FIG. Others are the same.
For example, in the case of 2), the origin is moved to the point 5 and it is determined whether particles exist on the m ₃ side by the method described above. If there is a particle on the m ₃ side under all the conditions 1) to 6) in which the above determination is made for all the faces of the element region, the particle exists inside the hexahedron. The above determination is performed based on the position of each particle and the element information of each element region, the element number of the initial element region including each particle is calculated, and stored in the RAM or the like as particle information of each particle.

  After step S112, the number of particles included in each element region is calculated in step S114 of FIG. For example, the number of particles in each element region is calculated by counting the number of particles having the element number of the element region as particle information.

  Here, an interpolation function and a gradient function for obtaining physical information of a node based on the particle information of each particle, or performing the inverse calculation, used in the subsequent processing will be described.

Now, assume that the particle p exists in an element region. As shown in FIG. 6, in order to calculate an interpolation function or the like, the global coordinate system (x, y, z axes) is converted to the local coordinate system (r ₁ , r ₂ , r ₃ axes), Consider that the element after conversion is an isoparametric square element. The local node numbers of the element regions converted into the local coordinate system are determined in advance for each element region as indicated by numbers 1 to 8 in FIG. 6 and FIG. 2 (b-2), and the coordinates X _p (k of the particle p ) As a function of shape function S _i (X _p (k)) (= S _ip (k) = {S _xip (k), S _yip (k) S _zip (k)} ^T ) is calculated by the following equation. .
Further, a gradient function ∇ S _ip (k) (= G _ip (k) = {G _xip (k), G _yip (k) G _zip (k)} ^T ) is calculated by the following equation.
Here, i is a local node number. R ₁ , r ₂ , and r ₃ represent three axes of the local coordinate system based on isoparametric square elements, and have a relationship represented by the following expression with respect to the overall coordinate system. Especially in the XY plane
There is a relationship. x _c and y _c respectively indicate the positions of the origin of the local coordinate system of the isoparametric square element in the global coordinate system. l _x and l _y are values of ½ of the length in the X-axis direction and the Y-axis direction in the global coordinate system of isoparametric square elements, respectively.
The above interpolation function and gradient function can be obtained based on the particle position of each particle and the element information of the element region where each particle is located. The shape function calculated by the isoparametric square element can be used as an interpolation function in the calculation process described later. For this reason, the symbols S _ip and N ⁱ _pj indicating the shape function are sometimes referred to as interpolation functions.

Returning to FIG. 4, before starting the time calendar calculation (k = 0), after step S114 described above, in step S115, the initial value (k = 0, initial value) of the stress σ _p (k) of each particle Stress) is calculated and stored as particle information 110 in a RAM or the like.

  In the present embodiment, in step S115, the particle stress is obtained from the stress distribution in the element region of the element 105 of the earth pillar 30. In the element region of the earth pillar 30, the stress is constant at the top and bottom inside, and the stress of the particles in the element region is determined by the stress of the upper four nodes. In the initial stage, the stress placed at the upper four nodes is based on the gravity of the particles, and the initial stress of the particles in the element region 105 is determined thereby.

  Specifically, the initial nodal force of a node is calculated based on the particle mass using an interpolation function, and then the initial particle force based on the initial nodal force, particle mass, and density of the node using a gradient function. The vertical stress is calculated, and the initial stress of the particle is calculated from the initial vertical stress. This will be described below.

In step S115, first, the nodal force associated with the weight of the particle is obtained using the particle mass and the interpolation function. That is, when attention is paid to one element region, the influence (force) of each particle on each node 120 is expressed by the following equation for the node 120 above the element region:
For the lower node 120:
Each is represented. Here, N ⁱ _pj is an influence function of the particle j with respect to the node i (i is a local node number), and a displacement shape function of the node i is used. As described above, the shape function calculated by the isoparametric square element is equal to the interpolation function.

_{Further, M Pj} is the mass of the particle j, _{g z} are weight acceleration. Both of the above formulas are equivalent to formulas used in a commonly used FEM. Superimposing the nodal force of nodes determined by the above equation between vertical element regions (plus match) soil column by its own weight of the particles by the initial nodal force f _{i (i} = 1 to 4) is calculated.

On the other hand, the formula usually used in FEM
Thus, the relationship between stress and nodal force can be described by the following equation using MPM. In the following equation, the differential symbol ∂ is the side differentiation in the z direction. In the above equation, ε _vp is the particle strain, and σ _pi is the particle stress.
For the relationship between the stress at the node and the stress at the particle,
For nodal forces and particle stresses
There is a relationship.

When Equation (15) and Equation (16) are connected,
It becomes. The matrix A represents a matrix for converting stress into nodal force. Since the matrix A is a square matrix, an inverse matrix of A can be obtained, and by solving this, the stress σ _{i at the} node position can be calculated from the nodal force f. In step S115, it obtains the initial vertical stress sigma _pi particles by substituting this manner (initial) nodal force f _i from the obtained node positions (initial) stress sigma _i in the above equation (14).

In order to obtain the inverse matrix of A, it is necessary that nx ≧ 2 and ny ≧ 2. Here, nx and ny are the numbers of particles in the x direction and y direction in one element, respectively, and can be obtained using the number of particles in the element region and the x direction length and y direction length of the element region. . If the inverse matrix cannot be obtained, the components of the A matrix are diagonalized as shown in the following equation to obtain an inverse matrix.

Thereafter, the initial stress of the particles in the element region is σ _z = σ _Pi , σ _X = σ from the initial vertical stress σ _pi and the lateral pressure coefficient K ₀ which is a value input in advance to the ground deformation analysis apparatus 1 as an input condition. _Pi K ₀ and σ _Y = σ _Pi K ₀ are calculated and stored in the RAM or the like as particle information. Here, the shear components are all zero.

Incidentally, the lateral pressure coefficient K ₀ is preferentially utilized if set non-zero. When the input is zero, K ₀ = ν / (1−ν) determined from the Poisson's ratio ν in the case of elasticity calculation. In the case of elasto-plastic analysis, the static earth pressure coefficient K ₀ = 1−sin φ determined from the internal friction angle φ of the soil is used.

Next, returning to FIG. 4, the calculation of the element number during the time calendar calculation will be described.
As will be described later, in the MPM, the position of the particle is updated at each time step. Therefore, during the time calendar calculation in which k ≠ 0 (No in step S111), in step S113, the element region including the particle and the local coordinates in the element region are updated based on the particle position.
In this embodiment, using the particle position and element information, the vertical coordinate of the bottom surface of the element region corresponding to the XY plane position of the particle and the length of the element region in the Z direction (vertical direction) To calculate the element region in which the particle is located. When the particle is actually outside the element region, the particle position is converted into a predetermined position on the outer surface of the element region corresponding to the particle position. This will be described below with reference to FIGS.

In this embodiment, the position of the element region on the XY plane of the particle p can be determined relatively easily because the element region forms an Euler lattice on the XY plane. As shown in FIG. 7A, for example, the positions of the element regions in the X and Y directions of the particles at the positions x and y are (x−X _s ) / L _x +1,
It can be calculated as (y−Y _s ) / L _y +1. Here, X _s and Y _s are the coordinates of the lower left node of the element number (1, 1, Nz), and L _x and L _y are the lengths of the element regions in the X and Y directions, respectively.

  On the other hand, as described above, in the present embodiment, as shown in FIG. 7B, the element region including the particles is an Euler lattice on the XY plane and the element regions having the same shape overlap in the Z direction. May not form an Euler lattice in the Z direction. Therefore, the position Nz of the element region in the Z direction of the particle in the Z direction position z is calculated by the following method.

  That is, since the XY plane is an orthogonal lattice, the element region on the XY plane can be determined by the above-described method. However, in the Z direction, since the undulation of the bottom surface is taken into consideration, the same method as the XY plane cannot be used. For the Z direction, an element region including particles and local coordinates in the element region are specified by the following procedure.

First, using the above-described (x−Xs) / Lx + 1, (y−Ys) / Ly + 1 and element information, among the element regions at the bottom of each earth pillar, the bottom XY region is the planar position of the particle p ( An element region 105a-1 including x, y) is calculated.
Next, using the element information of the element region 105a-1, corresponding to the planar position (x, y) of the particle p on the bottom surface of the element region 105a-1 indicated by a cross in FIG. 7B. The Z direction coordinate Ze to be calculated is calculated. And, using the vertical position Z of the particle p, the following equation
N that satisfies the condition shown in FIG. Lz is the height in the Z direction of the element region 105a-1 (105a-2). Thereby, it can be specified that particles are included in the nth element region from the lowest element region 105a-1. In this way, the element region where the particle is located is calculated and updated based on the particle position and the element information. For example, when n = 2, it is assumed that the particle p is included in the second element region 105a-2 from the lowermost element region 105a-1.

  As described above, the element region may not form an Euler lattice in the Z direction. Therefore, as shown in FIG. 7B, the particle p does not actually exist in the element region obtained as described above. There is a case.

  Therefore, in the present embodiment, after calculating the element region including (assuming) the particle p, the element region is an isoparametric square element by Newton-Raphson method as described with reference to FIG. In this case, the coordinate transformation is performed for each region viewed from the Z axis positive side with respect to the bottom surface 130 (region AA) of the element region shown in FIG. For convenience, the following position conversion is performed.

  For example, when the particle p is in the AA plane on the XY plane in FIG. 8A, the particle p is present when the particle p is above the element region 105 after conversion, as indicated by P1 in FIG. 8B. It is assumed that the point is on a point on the upper surface of the element region 105 that is perpendicular to the upper surface. On the contrary, when it is below the element region 105, it is assumed that it is at a point on the lower surface where a perpendicular is dropped from the particle p toward the lower surface (bottom surface 130) of the element region 105. As indicated by P <b> 2 in FIG. 8B, when the height of the particle p is in the element region, the particle p is included in the element region 105, and the position is left as it is.

  On the other hand, when the particle p is not in the AA plane but in the A1 to A4 plane in FIG. 8A, the particle p is higher than the isoparametric element after conversion, as indicated by P3 in FIG. 8B, for example. Is located at a point on the upper side obtained by dropping a perpendicular from the particle p to the corresponding side (upper side) of the upper surface of the element region 105. Conversely, when it is below the element region 105, it is assumed that it is at a point on the lower side where a perpendicular line is drawn from the particle p toward the corresponding side (lower side) of the lower surface of the element region 105. As indicated by P4 in FIG. 8 (b), when the height of the particle p is in the element region, it is at a point on the side surface that is perpendicular to the corresponding side surface of the element region 105. To do.

  Further, in FIG. 8A, when not in the AA plane but in the B1 to B4 plane, for example, as shown by P5 in FIG. 8B, when the particle p is above the isoparametric element after conversion. The particle p is assumed to be at a corresponding corner (upper corner) of the upper surface of the element region 105. Conversely, when it is below the element region 105, the particle p is assumed to be at a corresponding corner (lower corner) of the lower surface of the element region 105. As indicated by P6 in FIG. 8B, when the height of the particle p is in the element region, the vertical line is drawn toward the corresponding side (upper and lower sides) of the element region 105. Suppose that

  In this manner, the position is converted so that the position of the particle is on the outer surface of the element area in accordance with the position of the particle with respect to the element area after conversion. The above conversion can be performed using the particle position and element information. Thereafter, in step S114, the number of particles in each element region is calculated and updated as described above.

  Return to the explanation of the flow of ground deformation analysis. In step S120, the nodal mass and nodal force are calculated in the element region, the nodal acceleration is calculated from the nodal mass and nodal force, and the nodal first velocity is calculated from the nodal acceleration.

In step S120, first, the mass m _i (k) of the node is calculated by the following equation using the mass M _{p of the} particle _p and the interpolation function S _ip .
Here, N _p is the total number of particles included in the element region.

  Next, as a nodal force, an internal force and an external force are calculated using an interpolation function and a gradient function based on the particle mass, stress, density, volume force, and contact force.

First, the nodal force f _i ^int (k) = {f _xi ^int (k), f _y ^int (k), f _zi ^int (k)} ^T as the internal force is defined as the mass of the particle, the stress σ _{ij p} (k), Using the gradient function G _nip based on the density ρ _p (k),
Calculate with

On the other hand, nodal force as the external force _f ⁱ ext _{^(k) =} body force _b i (k) in the MPM and ^{{f xi ext (k),} f yi ext (k), f zi ext (k)} T = {b _xi (k), b _yi (k)} ^T and contact force τ _i (k) = {σ _p (k) ^T · n},
Calculate with here,
It is. n represents a normal vector of dS. The position vector of X _p the above particles, ∂Omega is a region in which the action of contact stress, dS represents a unit area.

Next, the acceleration of the node is obtained based on the mass of the node and the nodal force. That is, the node acceleration a _i (k) = {a _xi (k), a _yi (k), a _zi (k)} ^T is expressed by the following equation of motion of the node:
Calculate from

Finally, based on the acceleration of the node, the first velocity v _i (k + 1) of the node at the next step k + 1 is
To calculate and update. Δt is the length of the time step k and can be determined in advance as appropriate.

Next, in step S130, the particle position X _p (k + 1), velocity v _p (k + 1), and strain increment Δε _{ij p} (k + 1) are calculated.

In step S130, first, the position of the particle is calculated by the following equation using an interpolation function based on the first velocity of the node, and updated.

In addition, the velocity of the particle is calculated by the following equation using an interpolation function based on the acceleration of the node, and updated.

Further, smoothing of the velocity v _i (k + 1) proposed by Sulsky is performed, and the node velocity (second velocity) v _i (k + 1) is obtained again by the following equation. That is, the velocity of the node (second velocity) is obtained by the following equation using the interpolation function based on the mass of the particle, the velocity, and the mass of the node.

Finally, the particle strain increment Δε _{ij p} (k + 1) is calculated from the following equation using the gradient function based on the nodal velocity. Note that the subscripts m and n represent velocity components.
However, regarding engineering strain Δγ _mn ,
It is.

Subsequently, in step S140, the particle strain ε _{ij p} (k + 1), stress σ _{ij p} (k + 1), and density ρ _p (k + 1) are calculated and updated based on the particle strain increment.

In step S140, first, a particle stress increment Δσ _{ij p} (k + 1) is calculated based on the particle strain increment using the object stiffness matrix T _ijkl by the following equation. The subscripts k and l represent distortion components.

Further, the strain of the particle is calculated by the following formula based on the increment of the strain of the particle and updated.

The particle stress is calculated and updated based on the particle stress increment.
It is known that it is preferable to consider the influence of rotation when performing large deformation analysis. Here, the Jamann rate is used to take this effect into account. At this time, the stress σ _{ij p} (k + 1) of the particle p at the step k + 1 is calculated by the following equation using the particle velocity, the interpolation function, etc. in addition to the particle stress increment and updated.
here
It is.

Further, the density of the particles is calculated by the following equation using the particle strain increment and updated.

  Thereafter, until the calculation for all time steps is completed (No in step S150), the value of the time step number k is incremented (step S160), and the process returns to step S110 to perform the next time calendar calculation. In the course of this processing, particle information that is a Lagrangian variable of each particle is updated. For example, the particle strain is updated in step S140. In this way, when the time calendar calculation is repeated by the predetermined number of time steps to finish the time calendar calculation (Yes in step S150), the ground analysis is finished, and the position, strain, stress, etc. of each particle such as position, strain, stress, etc. are displayed on the display unit. Various values are displayed in different colors. In addition, when the above series of processes is completed in one time step, the process returns to step S110 as described above, and the element number where each particle is located is updated in step S113 to identify the particle in each element region. At this time, using the (original) element number of each particle, for each element region, only the particle having the element number of the neighboring element region can be processed, thereby reducing the amount of calculation.

  In the present embodiment, as MPM, instead of the processing in steps S120 to S140 described above, the time increment in one time step called USAVG is halved, and the strain, stress, and density of the particle are calculated in two steps. You can also take a method. This will be briefly described below with reference to FIG.

  The initial condition setting in step S200 shown in FIG. 9 and various calculation processes related to the particle position in step S210 are the same as the processes in steps S100 and S110 described above. In USAVG, first, in step S220, the velocity of the node is calculated.

In step S220, first, the mass m _i (k) of the node is calculated by the above equation (20), and then the velocity of the node is calculated by the following equation.
This is the same processing as that for obtaining the node speed (second speed) in step S130 described above.

Next, in step 230, the particle strain increment Δε _{ij p} (k) is calculated from the following equation.
This is the same processing as that for obtaining the particle strain increment in step S130 described above.

  Next, in step S240, the particle strain, stress, and density at time step k + 1/2 are calculated from the particle strain increment.

First, the particle stress increment Δσ _{ij p} (k) is calculated by the following equation.
Also, the strain of the particle is calculated by the following formula and updated.
And the stress of particle | grains is computed by following Formula.
here,
It is.
Further, the density of the particles is calculated by the following formula.
The above process is the same as the process in step S140 described above.

  In step S250, the nodal force is calculated in the element region, the nodal acceleration is calculated from the nodal mass and the nodal force, and the first velocity of the nodal is calculated from the nodal acceleration.

First, the mass of nodal force _f ⁱ int (k + 1) particles as an internal force, stress, by using a gradient function _{G an nip} based on the density equation
The nodal force f _i ^ext (k + 1) as an external force is calculated by using the bulk force b _i (k + 1) and the contact force τ _i (k + 1) in the MPM as follows:
Calculate with here,
It is. Next, based on the mass of the node and the nodal force,
To obtain the acceleration a _i (k + 1) of the node. Finally, based on the acceleration of the node, the first velocity v _i (k + 1) of the node at the next step k + 1 is
Calculated by The above processing corresponds to the processing in step S120 described above.

  Next, in step S260, the position, velocity, and strain increment of the particle are calculated.

In step S260, first, the particle position X _p (k + 1) is calculated by the following equation using the interpolation function based on the first velocity of the node and updated.
Further, the velocity v _p (k + 1) of the particle is calculated by the following equation using the interpolation function based on the acceleration of the node, and is updated.
Further, smoothing of the velocity v _i (k + 1) proposed by Sulsky is performed, and the node velocity (second velocity) v _i (k + 1) is obtained again by the following equation. That is, the velocity of the node (second velocity) is calculated and updated by the following equation using an interpolation function based on the mass, velocity, and mass of the particle.
Finally, the particle strain increment Δε _{ij p} (k + 1) is calculated from the following equation using the gradient function based on the nodal velocity.
The above processing corresponds to the processing in step S130 described above.

  Subsequently, in step S270, the strain, stress, and density of the particle are calculated again based on the increment of the strain of the particle and updated.

In step S270, first, the particle stress increment Δσ _{ij p} (k + 1) is calculated by the following equation based on the particle strain increment using the object stiffness matrix T _ijkl .
Also, the particle strain ε _{ij p} (k + 1) is calculated by the following formula based on the particle strain increment and updated.
The particle stress σ _{ij p} (k + 1) is calculated and updated based on the particle stress increment. The stress σ _{ij p} (k + 1) of the particle p is calculated by the following equation and updated as in the previous case.
here
It is. Further, the particle density ρ _p (k + 1) is calculated by the following equation using the particle strain increment and updated.
The above processing corresponds to the processing in step S140 described above.

  Thereafter, until the calculation for all time steps is completed (No in step S280), the value of the time step number k is incremented (step S290), and the process returns to step S210 to perform the next time calendar calculation. In the course of this processing, particle information that is a Lagrangian variable of each particle is updated. In this way, when the time calendar calculation is repeated by the predetermined number of time steps to finish the time calendar calculation (Yes in step S280), the ground analysis is ended, and the position, strain, stress, etc. of each particle such as position, strain, stress, etc. are displayed on the display unit. Various values are displayed in different colors.

  Next, consideration of water pressure when using the MPM analysis method of the present embodiment for ground analysis will be described. In order to use MPM for ground analysis, considering the water pressure, it becomes possible to incorporate groundwater behavior such as liquefaction and osmotic flow into the analysis.

In the Darcy's law equation expressed in the u-P format shown below, it is said that the calculation is more stable when the acceleration term is removed from the term of the flow velocity used.
Here, k is permeability, g is the gravitational acceleration, n represents porosity, k _f is the bulk modulus of the water, the [rho _f shows the density of water. Assuming non-drainage, there is no acceleration related to water permeability,
It becomes. Furthermore, considering that water is an elastic body, equation (60) is integrated with respect to time,
It can be expressed. However, t = 0, εν = 0, and P = 0. εν is a volume strain and can be obtained from the strain of particles.

Here, in the MPM, the water pressure can be obtained from the particle strain or the particle strain increment. The water pressure can be included as an internal force when determining the nodal force. For example, in the above-mentioned USAVG, after obtaining the particle strain at the time step k + 1/2 in the above step S240,
The water pressure P (k + 1/2) at time step k + 1/2 is calculated by the following equation, or the particle strain increment calculated in step S230 is used to calculate
To calculate and update the water pressure P (k + 1/2). These water pressures are stored in a RAM or the like. In step S250, the generated water pressure may be included in the internal force of the node by the following equation.
That is, the term P (k + 1/2) is added to the normal USAVG described above. Here, δ _nm represents the Kronecker delta.

Modulus k _f is sufficiently large and the volume trε soil water becomes almost zero, the calculation does not occur in the volume change. In FEM, the average principal stress and water pressure are evaluated at the center point of the element to match the stress and water pressure evaluation points used for nonlinear analysis. In the particle method, the stress evaluation point is on the particle. However, by setting the evaluation point of the water pressure to the center position of the element region regardless of the position of the particles, it is possible to maintain the state represented by the formula (61) in all the particles. Specifically, in the strain calculation and the internal force calculation, only the volume strain ε _V is a selective integration that is calculated at the center of the element region.

Next, the extension to the dynamic analysis when the analysis method by the particle method is used for the ground analysis will be described. For example, in the ground deformation analysis method of the present embodiment, a downward viscous boundary can be set. The downward viscous boundary is expressed by attaching a dashpot as shown in FIG. 10A to a node, that is, using the following equation as the equation of motion of the node. In the dashpot, the equation of motion is as follows, assuming that the earthquake acceleration a ″, the gravitational acceleration g, the relative acceleration u ″, and the spring constant k.
In the present embodiment, for the necessary nodes, the above equation (65) can be used as the equation of motion of the nodes. In the above-described flow, the first velocity of the nodal point is obtained based on the above-described relative acceleration.

As for the nodes of the viscous boundary, as shown in FIG.
The internal force of the node is calculated by In FIG. 10B, R _i is a reaction force, and the reaction force at the node 1 is R ₁ . At the initial stress described above, f _i ^int and R _{i are} balanced, so that the acting force becomes zero and the node does not move.

Therefore, as shown in FIG. 10C, the acting force due to the viscous boundary is obtained from the dynamic stress change (σ−σ _i ). Therefore, at the nodes where the viscous boundary is required, the initial stress and the gravity are taken. The viscous force calculated using the stress change as the applied force is added to the internal force of the node in the equation of motion.

In addition, the fixed boundary of the ground can be expressed by setting the velocity v _i and acceleration a _i in the above-mentioned flow to 0 at the node of the bottom surface of the predetermined element region.

  Furthermore, an element region can be specified as a side boundary condition, and a direction in which particles cannot move within the element region can be specified. As for the movement of particles in each element region, for example, if the movement in the Z direction (up and down direction) is set and the particle position in the Z direction is fixed, the particles do not move up and down. Since the movement of the particles is forcibly set to zero, the same strain as that when moving is generated and a stress is generated. However, the particles remain, and a reaction force to the surrounding particles due to the contact force is generated.

Also, Rayleigh attenuation can be applied to each node. In the case of Rayleigh damping, an equation of motion simulating a dashpot having a damping matrix with values obtained by multiplying the mass matrix and the stiffness matrix by coefficients α and β is used for calculating the acceleration of the node. For example, in the case of the above-mentioned USAVG, the equation of motion can be solved by using any one of the following equations as the internal force of the node in k + 1/2 steps.
Here, D represents a stiffness matrix.

  Next, an example of analysis using a three-dimensional MPM by the ground analysis apparatus will be shown. FIG. 11 shows an example in which a three-dimensional analysis of the flow of dry sand is performed by the ground analysis apparatus. FIG. 11 (a) shows a state before flow, FIG. 11 (b) shows a state after 5 seconds of flow generation, and FIG. 11 (c) shows a state after flow. It can be seen that the flow of dry sand flowing out of the slit can be reproduced in three dimensions.

  On the other hand, Fig. 12 shows an example of three-dimensional dynamic analysis of large deformation of the slope using this ground deformation analysis device, and shows the maximum shear strain as well as the ground deformation for the large deformation behavior when the slope receives earthquake motion. It is a thing. 12 (a) to 12 (c) show cases where the slope gradient is gentle, FIG. 12 (a) shows the state before the collapse, FIG. 12 (b) shows the state after 0.5 seconds from the start of the collapse, and FIG. ) Shows the state after collapse. FIGS. 13A to 13C show a case where the slope is steep, FIG. 13A is a state before the collapse, FIG. 13B is a state where 0.5 seconds have passed since the collapse, FIG. (C) shows the state after collapse. It is shown that the large deformation behavior of the slope differs depending on the slope gradient. It can be seen that when the slope gradient is gentle, the slope mass is deformed slowly, and when the slope slope is steep, it behaves rapidly.

  The preferred embodiments of the ground deformation analysis apparatus and the like according to the present invention have been described above with reference to the accompanying drawings, but the present invention is not limited to such examples. It will be apparent to those skilled in the art that various changes or modifications can be conceived within the scope of the technical idea disclosed in the present application, and these naturally belong to the technical scope of the present invention. Understood.

1 ... Ground deformation analysis device 100 ... Particle 105 ... Element 110 ... Particle information 120 ... Node

Claims (7)

A ground deformation analysis device for performing ground deformation analysis using a particle method,
At least the identification information of the element region and the position of the node forming the element region are stored as element information related to the element region having a three-dimensional shape, and at least the particle identification information and the element region where the particle is positioned as particle information about the particle Storage means for storing the identification information, position, mass, density, stress of
Based on the position of the particle and element information, the particle is located adjacent to the nodal point, with the nodal point forming the one surface as the origin and the one side forming the one surface with respect to the one side of the element region having the three-dimensional shape. A determination is made as to whether the particle is located in the direction indicated by the outer product of two vectors going to each of the two nodes, and if the particle is located in the direction indicated by the outer product for all faces of the element region, An initial element region calculating means for calculating the element region as a positioned element region, and calculating the number of the particles in the element region;
The initial nodal force of the node is calculated based on the mass of the particle using an interpolation function, and the initial vertical direction of the particle based on the initial nodal force of the node, the mass of the particle, and the density using a gradient function. Initial stress calculating means for calculating stress and calculating initial stress of the particle based on initial vertical stress of the particle;
For the element region, the interpolation function is used to calculate the mass of the node based on the mass of the particle, and the interpolation function and the gradient function are used to calculate the mass of the node based on the mass, density, and stress of the particle. A nodal velocity calculating means for calculating a nodal force, calculating an acceleration of the nodal based on the mass of the nodal point and the nodal force of the nodal point, and calculating a first velocity of the nodal based on the acceleration of the nodal point;
For the element region, the position of the particle is calculated based on the first velocity of the node using an interpolation function, and the velocity of the particle is calculated based on the acceleration of the node using an interpolation function. A position / strain increment calculating means for calculating a strain increment of the particle based on the velocity of the particle using a gradient function;
For the element region, strain / stress / density calculating means for calculating strain, stress, density of the particle based on strain increment of the particle;
Based on the position of the particle and the element information, when the particle is located from the vertical coordinate of the bottom surface of the element region corresponding to the planar position of the particle and the vertical length of the element region An element region to be calculated, and when the particle is outside the element region, the position of the particle is converted to a predetermined position on the outer surface of the element region according to the position of the particle Area calculation means;
A ground deformation analysis apparatus comprising:   The ground deformation analysis apparatus according to claim 1, further comprising a water pressure calculation unit configured to calculate a water pressure using the strain of the particles or the strain increment of the particles.   The ground deformation analysis apparatus according to claim 1 or 2, wherein the nodal velocity calculating means calculates an acceleration of the nodal point using a motion equation simulating a dashpot for some of the nodal points. At least the identification information of the element region and the position of the node forming the element region are stored as element information related to the element region having a three-dimensional shape, and at least the particle identification information and the element region where the particle is positioned as particle information about the particle An information processing apparatus having storage means for storing identification information, position, mass, density, and stress, and performing ground deformation analysis using a particle method,
Based on the position of the particle and element information, the particle is located adjacent to the nodal point, with the nodal point forming the one surface as the origin and the one side forming the one surface with respect to the one side of the element region having the three-dimensional shape. A determination is made as to whether the particle is located in the direction indicated by the outer product of two vectors going to each of the two nodes, and if the particle is located in the direction indicated by the outer product for all faces of the element region, An initial element region calculating step of calculating the element region as a positioned element region and calculating the number of the particles in the element region;
The initial nodal force of the node is calculated based on the mass of the particle using an interpolation function, and the initial vertical direction of the particle based on the initial nodal force of the node, the mass of the particle, and the density using a gradient function. Calculating an initial stress of the particle based on the initial vertical stress of the particle;
For the element region, the interpolation function is used to calculate the mass of the node based on the mass of the particle, and the interpolation function and the gradient function are used to calculate the mass of the node based on the mass, density, and stress of the particle. Calculating a nodal force, calculating an acceleration of the nodal based on the mass of the nodal point and the nodal force of the nodal point, and calculating a first velocity of the nodal based on the acceleration of the nodal point; and
For the element region, the position of the particle is calculated based on the first velocity of the node using an interpolation function, and the velocity of the particle is calculated based on the acceleration of the node using an interpolation function. A position / strain increment calculating step of calculating a strain increment of the particle based on the velocity of the particle using a gradient function;
For the element region, a strain / stress / density calculation step for calculating strain, stress, and density of the particle based on strain increment of the particle;
Based on the position of the particle and the element information, when the particle is located from the vertical coordinate of the bottom surface of the element region corresponding to the planar position of the particle and the vertical length of the element region An element region to be calculated, and when the particle is outside the element region, the position of the particle is converted to a predetermined position on the outer surface of the element region according to the position of the particle Region calculation step;
A ground deformation analysis method comprising:   The ground deformation analysis method according to claim 4, wherein the information processing apparatus further performs a water pressure calculating step of calculating a water pressure using the strain of the particles or the strain increase of the particles.   6. The ground deformation analysis method according to claim 4, wherein, in the nodal velocity calculating step, acceleration of the nodal point is calculated for a part of the nodal points using an equation of motion simulating a dashpot.   A program for causing an information processing device to function as the ground deformation analysis device according to any one of claims 1 to 3.