Two Phase MPM Solver

Overall Structure

The structure of the Two-phase MPM solver is similar to the other MPM schemes. The key differences of the input .JSON file are highlighted below. The two-phase solver implemented in the current code version is the single-point two-phase formulation.

Mesh

The mesh object is generally the same as the explanation given in the explicit scheme with major differences described as follow.

{

	"mesh": {
		"mesh": "mesh.txt",
		"check_duplicates": false,
		"isoparametric": false,
		"entity_sets": "entity_sets.json",
		"cell_type": "ED3H8",
		"io_type": "Ascii3D",
		"node_type": "N3D2P",
		"particle_cells": "particles-cells.txt",
		"boundary_conditions": {
			"velocity_constraints": [
				{
					"nset_id": 0,
					"dir": 2,
					"velocity": 0.00
				},
				{
					"nset_id": 0,
					"dir": 5,
					"velocity": 0.00
				},
				{
					"nset_id": 1,
					"dir": 0,
					"velocity": 0.00
				},
				{
					"nset_id": 1,
					"dir": 3,
					"velocity": 0.00
				},
				{
					"nset_id": 2,
					"dir": 1,
					"velocity": 0.00
				},
				{
					"nset_id": 2,
					"dir": 4,
					"velocity": 0.00
				}
			],
			"pressure_constraints": [
				{
					"phase_id": 1,
					"nset_id": 3,
					"pressure": 0.00
				}
			]
		},
		"particles_pore_pressures": {
			"dir_h": 0,
      			"dir_v": 1,
      			"type": "water_table",
      			"water_tables": [
        		{
          			"h0": 5.0,
          			"position": 0.0
        		},
        		{
          			"h0": 3.0,
          			"position": 1.0
        		}]
    		}
	}

}

• "node_type" is either "N2D2P" or "N3D2P".
• "velocity_constraints" follow the convention described in the explicit scheme with direction "dir" should be set as (0|1|2) for the solid phase constraints and (3|4|5) for the fluid phase constraints in 3D. Meanwhile, in 2D, they should be set as (0|1) for the solid phase and (2|3) for the fluid phase.
• "pressure_constraints" follow the convention described in the Navier-Stokes scheme with "phase_id" to be set at 1.

Particles Pore Pressures

Optional particle pore pressure initializer can be input via a text file (in a similar format to particle stresses). The first row of the text file provides the total number of material points. Each subsequent row provides the particle id and corresponding value of pressure.

Alternatively, if the pore pressure in the entire simulation domain is isotropic, the "isotropic" option can be specified following this format:

{
    
        "particles_pore_pressures": {
            "type": "isotropic",
            "values": 10000000
        }

}

where the "values" are the pore pressure with positive magnitude indicating compression.

One can also define a hydrostatic pressure evolving with depth by using the option "water_table". Here, one needs to specify the following information to help the code to define the initial water table:

• "dir_v" is the vertical direction (Gravity direction) of the water table.
• "dir_h" is the horizontal direction of the water table.
• "water_tables" is the list of reference points indicated by "position" in the horizontal axis defined in "dir_h"-direction and "h0" is the initial water height in the direction of "dir_v". Here, the datum is always measured at the zero coordinate of the "dir_v"-axis.
• It is important to note that to use the "water_table" option, one needs to define "gravity" in the "external_loading_conditions".

Particles

The particle type should be set as "particle_type": "P2D2PHASE" or "P3D2PHASE". Meanwhile, the "material_id" should be set as a list with the first entry indicating the solid-phase material and the second entry indicating the fluid-phase material. Below is the example of the "particles" object:

	"particles": [
		{
			"generator": {
				"check_duplicates": false,
				"location": "particles.txt",
				"io_type": "Ascii3D",
				"particle_type": "P3D2PHASE",
				"material_id": [
					0,
					1
				],
				"pset_id": 0,
				"type": "file"
			}
		}
	],

Materials

The "materials" object should include at least two materials that define both the solid and fluid phases. The parameters defined in the solid-phase material are assumed to follow the effective-stress analysis which solely characterizes the mechanical behavior of the solid skeleton of a porous media.

	"materials": [
		{
			"id": 0,
			"name": "solid",
			"type": "LinearElastic3D",
			"density": 2600.0,
			"poisson_ratio": 0.3,
			"youngs_modulus": 9.80E+7,
			"porosity": 0.3,
			"intrinsic_permeability": false,
			"k_x": 0.001,
			"k_y": 0.001,
			"k_z": 0.001
		},
		{
			"id": 1,
			"name": "water",
			"type": "Newtonian3D",
			"density": 1000,
			"bulk_modulus": 2.2E9,
			"dynamic_viscosity": 8.9E-4
		}
	],

The solid-phase material should include additional information, such as:

• "porosity" which indicates the initial material averaged porosity.
• "k_x","k_y", and "k_z" which indicates the permeability in different axes.
• "intrinsic_permeability" is a boolean that defines the unit of "k_x","k_y", and "k_z", which can be defined as intrinsic permeability (m^2) (true) or soil-mechanics permeability (m/s) (false).

Analysis

The two-phase MPM "analysis" object is described as follow:

{

	"analysis": {
		"type": "MPMSemiImplicitTwoPhase3D",
		"scheme_settings": {
			"beta": 1
		},
		"velocity_update": false,
		"pressure_smoothing": false,
		"pore_pressure_smoothing": false,
		"linear_solver": {
			"assembler_type": "EigenSemiImplicitTwoPhase3D",
			"solver_settings": [
				{
					"dof": "pressure",
					"solver_type": "KrylovPETSC",
					"sub_solver_type": "cg",
					"preconditioner_type": "jacobi",
					"max_iter": 100,
					"tolerance": 1E-5,
					"verbosity": 0
				},
				{
					"dof": "acceleration",
					"solver_type": "KrylovPETSC",
					"sub_solver_type": "bicgstab",
					"preconditioner_type": "none",
					"max_iter": 100,
					"tolerance": 1E-5,
					"verbosity": 0
				}
			]
		},
		"free_surface_detection": {
			"type": "density",
			"volume_tolerance": 0.25
		},
		"uuid": "semi_implicit_3D",
		"dt": 0.0001,
		"nsteps": 2001
	}

}

Type

Supported types:

Analysis	Description
MPMExplicitTwoPhase2D	explicit 2D Two-phase Single-layer MPM
MPMExplicitTwoPhase3D	explicit 3D Two-phase Single-layer MPM
MPMSemiImplicitTwoPhase2D	semi-implicit 2D Two-phase Single-layer MPM
MPMSemiImplicitTwoPhase3D	semi-implicit 3D Two-phase Single-layer MPM

If the semi-implicit scheme is picked

Scheme settings

Users can also specify the "scheme_settings" option to change the default parameter considered in the simulation.

• "beta" is the parameter that defines whether a full (0) or incremental (1) projection scheme is used.

Linear Solver

The "linear_solver" option is used to control the type of assembler and linear solver needed to solve the momentum equation implicitly.

• "assembler_type" should be "EigenSemiImplicitTwoPhase2D" or "EigenSemiImplicitTwoPhase3D" for 2D or 3D simulation.
• "solver_settings" contains all necessary settings to solve the linear equation Ax=b. Two solvers should be defined to solve the predicted acceleration and pressure Poisson equations.
• "dof" is the name of degrees of freedom. Here, "pressure" and "acceleration" should be defined for the semi-implicit scheme.
• "solver_type" is the name of the linear solver library used. Options available are:

Analysis	Description
DirectEigen	Direct solver of Eigen Library
IterativeEigen	Iterative solver of Eigen Library
KrylovPETSC	Iterative solver of PETSc Library

• "sub_solver_type" is the name of the linear solver within the selected library.
• "preconditioner_type" is the name of used preconditioner.
• "max_iter" is the maximum number of iterations for the iterative linear solver.
• "tolerance" is the relative residual tolerance.
• "abs_tolerance" is the absolute tolerance.
• "div_tolerance" is the divergence tolerance.
• "verbosity" is to control how much information is output in the console. The default value is 0.

Free-surface detection

The "free_surface_detection" option is used to detect nodes, particles, and cells nearby the material's free surface where zero pressure BC should be imposed.

• "type" can be "density" or "geometry". Only "density" is available for parallel solvers with MPI.
• "volume_tolerance" defines volume tolerance where all nodes should be considered as free-surface. The default value is 0.25.