Setting Forces and Fluxes

These commands are used within particle-based boundary condition blocks to set loads, tractions, heat fluxes, and concentration fluxes on particles

Load Conditions

The scripted Load command applies a force directly to each particle within the parents's block shape:

Load (dir),(style),(value),<(time)>

In XML files, the command is:

<LoadBC dir='(dir)' style='(style)' load='(value)' time='(time)'
                function='(function)'/>

where

• (dir) is 1, 2, or 3 to specify the loading direction as in the x, y, or z direction (In scripted files, (dir) can be x, y, or z or can be R or Z if axisymmetric).
• (style) specifies the style of the applied load.
• (value) and (time) -each load (style) depends on two parameters specified by (value) and (time). The (value) is required because there is no point in applying zero load. If (time) is not supplied, it is set to zero. The standard units are N for (value) and alt time units for (time) (but the units may change depending on the (style) setting). The final magnitude of the load depends on the LoadType setting. Note that when running axisymmetric calculations, the load is interpreted as N per radian.
• (function) - if the (style) is function (or 6), this attribute in XML files specifies a user defined function. In scripted files, the user defined function should replace the (value) argument and must be entered as quoted text. The function should evaluate to the desired force in N.

LoadType Options

By default, the (value) argument or the result returned by a function gives the load on each particle with that boundary condition (in force units). Alternatively, a load or function can evaluate to the total Newtons applied to all particles in the parents's block shape. To chose to load style in scripted files, include the command:

LoadType (loadType)

within the parents's block shape. If (loadType) is "net", the boundary conditions specifies to the total load spread out over all particles in the parents's block shape. Alternative, if (loadType) is "perParticle" (or omitted), the load for the boundary condition is applied to each particle. You can switch back and forth between net load and per-particle load in a single parents's block shape by using mutliple LoadType commands. Each load boundary condition uses the loading option specified by the most recent LoadType command (or uses "perParticle" is no LoadType command was used).

In XML files, the loading option is set by using the following commands within the parents's block shape:

<net/>
<perParticle/>

As expected, the "net" option is specified with a <net/> command and the "perParticle" option is specified with a <perParticle/> command. These commands can be intersperse to apply different loading options to various load boundary conditions.

These loading-option commands only apply to particle load boundary conditions and have no effect on other traction, heat flux, and concentration flux conditions that happen to be within the same parents's block shape.

Axisymmetric Loads

When applying forces to particles in axisymmetric simulations the loads are interpreted as N per radian. As a consequence, if you intend to load a surface with particle's having different radial positions, the particle loads will have to depend on position to achieve a uniform stress. For this type of loading, it is usually easier (and better) to use traction boundary conditions instead of load boundary conditions. If you need particle loads, here is a sample calculation. Imagine you are applying loads in the Z direction for a surface whose normal is in the Z direction and extends from r_i to r_o (see figure on the right). To get uniform stress, the particle force will have to be proportional to r_p (radial position of the particle) or F_p = k r_p. The goal is to find k to achieve a desired constant stress σ. The total vertical force on all particles (assuming two particles per cell side) will be

      [math]\displaystyle{ F_{tot} = \sum_p 2\pi k r_p = 2\pi k \sum_{i=1}^{n_p} \left( r_i + {2i-1\over 4} \Delta r\right) \approx \pi k n_p(r_o+r_i) }[/math]

where n_p is the total number of particles on the surface and Δr is the cell size in the r direction. The resulting uniform stress is

      [math]\displaystyle{ \sigma = {F_{tot}\over \pi(r_o^2-r_i^2)} = { k n_p \over r_o-r_i} }[/math]

Thus the required function (and the "function" style is needed) for particle force per radian to achieve uniform stress of σ is

      [math]\displaystyle{ F_p = {\sigma(r_o-r_i)x\over n_p} = {\sigma x \Delta r\over 2} }[/math]

where x is the radial position of the particle (you need to use x instead of R for functions in boundary condition commands). This result is simply the stress (σ) times the surface area per radian on top of a particle located at r_p = x when Δr/2 is the radial distance associated with that particle. For a different number of particles per cell side, change the "2" in the denominator to that number.

If the boundary condition is provided using "net" force instead of "per particle" force, the required function for particle force should be:

      [math]\displaystyle{ F_{net} = n_pF_p = \sigma(r_o-r_i)x }[/math]

This result applies for any number of particles per cell side.

Traction Conditions

The scripted Traction command sets traction load on the edge of each particle within the parents's block shape:

Traction (dir),(face),(style),(value),<(time)>

In XML files, the command is:

<TractionBC dir='(dir)' face='(face)' style='(style)' stress='(value)' time='(time)'
                function='(function)'/>

where

• (dir) is 1, 2, or 3 to specify the loading direction as in the x, y, or z direction (In scripted files, (dir) can be x, y, or z or can be R or Z if axisymmetric). For TractionBC, dir can be 11 or 12 to mean normal or shear traction relative to the selected face, respectively. A normal traction is positive if pointing out of the particle domain and a shear traction is positive when oriented in the counter-clockwise direction. Shear tractions (option 12) are not yet allowed in 3D calculations.
• (face) defines which face of the original particle domain is loaded by the traction. In 2D, the faces are the edges numbered 1 through 4 in the counter-clockwise direction with the bottom edge (with normal (0,-1)) being number 1. In 3D, 1 to 4 are the same faces (i.e., faces with normals (0,-1,0), (1,0,0), (0,1,0), and (-1,0,0) on the original particle domain, respectively), 5 is the bottom (with normal (0,0,-1)), and 6 is the top (with normal (0,0,1)).
• (style) specifies the style of the applied traction.
• (value) and (time) - each velocity depends on two parameters specified by arguments (value) and (time). The (value) is required because there is no point in applying zero traction. If (time) is not supplied, it is set to zero. The standard units are pressure units for (value) and alt time units for (time) (but the units may change depending on the (style) setting).
• (function) - if the (style) is function (or 6), this attribute in XML files specifies a user defined function. In scripted files, the user defined function should replace the (value) argument and must be entered as quoted text. The function should evaluate to the desired traction in pressure units.

Note that valid calculation of grid forces resulting from specified traction during 2D simulations requires that you have correctly specified the grid thickness and that the grid thickness matches the particle thickness that is set in Region commands.

Exact Tractions

By using the command

ExactTractions <(option)>

in script files or adding the element

<ExactTractions/>

to the <MPMHeader> element in XML input files, traction loads will be found by exactly integrating shape functions over the deformed particle edge. In script files, set (option) to "yes" to use exact tractions or "no" to revert to standard traction methods. If (option) is omitted, it is treated as "yes".

In theory, this approach is more accurate, especially at large deformations. In practice, it is difficult to confirm an improvement. This option only works for uGIMP and CPDI shape functions and 2D simulations. If used in other simulations, an error will be reported.

Difference Between Loads and Tractions

A Load condition considers the load (in N) to be applied directly to the particle, which in MPM is not really applied to the surface of the body, even when the particle is on the edge. During the MPM algorithm, this load is extrapolated to all nodes around the paraticle resulting in some of a modeled loading being applied internal to the object. It is common, therefore, to see artifacts in particles near edges that are loaded directly by forces. The MPM results away from the particles, however, reflect the edge load better. This form of external loading appears to be common to all MPM codes.

In contrast, a Traction condition applies the load (as a stress) on the actual surface by implementing a surface integral term that appears in MPM derivations. This integral involves integration of the traction times grid shape functions on the particle surface. This surface is found from the current particle deformation.

Traction conditions typically reduce artifacts that occur on particles loaded by Load conditions. Because this method applies forces closer to the actual edge, the potential for artifacts due to low-mass nodes is a concern. The code tries to screen out obvious problem nodes. If instability occurs, try Load conditions instead. Traction conditions will correctly detect one nearby crack and apply force to the appropriate velocity field, but if a Traction condition is interacting with more than one crack, the results will be unpredictable.

Another difference is that a Load condition applies a total force. In large deformation problems, the applied stress will change if the cross-sectional area changes. In contrast, a Traction condition applies a stress to the deformed particle domain. In other words, the stress will be a true stress that remains constant even if the particle cross sectional area changes.

Heat Flux Conditions

The scripted HeatFlux command applies heat fluxes to surfaces near particles within the parents's block shape (the particles are expected to be on the surface):

HeatFlux (mode),(face),(style),(value),<(time)>

In XML files, the command is:

<HeatFluxBC dir='(dir)' face='(face)' style='(style)' value='(value)'
                        time='(time)' function='(function)'/>

where

• (mode) determines the type of heat flux. The two options are:
  1. external - this mode applies a surface flux in heat flux units for transport rate of heat per unit area. A positive value is flux into the material while a negative value is flux out of the material.
  2. coupled - this mode calculates the heat flux from a function of the particle temperature.
• (face) - specifies the particle surface that is on the edge and will have the heat flux. In 2D, imagine a box around the initial particle. Surfaces 1 to 4 are the four edges of the box in the order bottom, right, top, and left (with normals (0,-1), (1,0), (0,1), and (-1,0)). In 3D, imagine a cube around the initial particle. Surfaces 1 to 4 are same as for 2D (now with normals (0,-1,0), (1,0,0), (0,1,0), and (-1,0,0)), 5 in the bottom z surface (normal = (0,0,-1)), and 6 is the top z surface (normal = (0,0,1)).
• (style) specifies the style of the applied heat flux. For coupled heat flux, the (style) must be function (or 6).
• (value) and (time) - for external heat flux, the value depends on two parameters specified by arguments (value) and (time) (if omitted, their values are set to zero). The standard units are heat flux units for (value) and alt time units for (time) (but the units may change depending on the (style) setting). For coupled heat flux, the flux is applied only for time greater than (time) (in alt time units)
• (function) - if the (style) is function (or 6), this attribute in XML files specifies a user defined function. In scripted files, the user defined function should replace the (value) argument and must be entered as quoted text. The function should evaluate to the desired heat flux in heat flux units.

Note that valid calculation of total flux resulting from specified heat flux during 2D simulations requires that you have correctly specified the grid thickness and that the grid thickness matches the particle thickness that is set in Region commands.

Coupled Heat Flux

This heat flux calculates the heat flux from a user-defined function. At each time step, the variable t in the function (which usually means time) will be replaced by the particle temperature. The function should evaluate to the desired heat flux in W/m^2. The function can additionally depend on particle position (x, y, z, D, T, R, or Z in length units), and/or clockwise particle rotation (q in radians). If the (time) parameter is used, the coupled heat flux is calculated only after that time.

Two uses of this flux condition are to implement convective and/or radiative boundary conditions. In convective cooling or heating, the function could be

      [math]\displaystyle{ {\rm function} = h(T_0-t) }[/math]

where h is heat transfer coefficient (in W/(m^2 K)), T₀ is a reservoir temperature, and t is used for surface particle temperature. In radiative cooling or heating, the flux function could be the Stefan-Boltzman law:

      [math]\displaystyle{ {\rm function} = \sigma\varepsilon(T_0^4 - t^4) }[/math]

where σ = 5.6704e-8 W/(m^2 K^4) is the Stefan-Boltzman constant and ε (<1) is the emissivity. In these functions, the physical terms should be replaced by their numerical values. If desired, however, they could also be functions themselves, such as to have h be a function of t or surface particle temperature. They cannot, however, depend on time, because t is used here for particle temperature instead of time.

Concentration or Pore Pressure Flux Conditions

When doing diffusion calculations, these conditions set concentration flux, but when doing poroelasticity calculations, these same conditions set pore pressure flux (entered as pore volume fraction flux). Only one of these options can be active in a single simulation. These boundary conditions are ignored if neither diffusion calculations nor poroelasticity calculations are activated.

When doing diffusion calculations, the scripted ConcentrationFlux command applies concentration fluxes to surfaces of particles within the parents's block shape (the particles are expected to be on the surface):

ConcentrationFlux (mode),(face),(style),<(value)>,<(timeOrRes)>

When doing poroelasticity calculations, the scripted PorePressureFLus command applies concentration fluxes to surfaces near particles within the parents's block shape (the particles are expected to be on the surface):

PorePressureFlux (mode),(face),(style),<(value)>,<(timeOrRes)>

In XML files, the concentration flux and pore pressure flux commands are the same (the simulation type determines which conditions are set):

<ConcFluxBC dir='(mode)' face='(face)' style='(style)' value='(value)'
                        time='(timeOrRes)' function='(function)'/>

In the above commands:

• (mode) determines the type of concentration flux. The two options are:
  1. external - This condition applies a surface flux. For diffusion enter flux in solvent flux units for transport rate of substance per unit area. For poroelasticity, enter flux in pore volume fraction/((length units)²-time units). A positive value is flux into the material while a negative value is flux out of the material
  2. coupled - this mode calculates the concentration flux from a function of the particle concentration potential for diffusion or from a function of the particle pore pressure for poroelasticity.
• (face) - specifies the particle surface that is on the edge and will have the concentration or pore pressure flux. In 2D, imagine a box around the initial particle. Surfaces 1 to 4 are the four edges of the box in the order bottom, right, top, and left (with normals (0,-1), (1,0), (0,1), and (-1,0)). In 3D, imagine a cube around the initial particle. Surfaces 1 to 4 are same as for 2D (now with normals (0,-1,0), (1,0,0), (0,1,0), and (-1,0,0)), 5 in the bottom z surface (normal = (0,0,-1)), and 6 is the top z surface (normal = (0,0,1)).
• (style) specifies the style of the applied concentration or pore pressure flux. For coupled concentration flux or coupled pore pressure flux, the (style) must be function (or 6). For diffusion, the variable (timeOrRes) should specify the solvent reservoir potential (0 to 1). For poroelasticity, the variable (timeOrRes) should specify the reservoir pressure (in pressure units).
• (value) and (timeOrRes) - for external fluxes, the value depends on two parameters specified by arguments (value) and (timeOrRes) (if omitted, their values are set to zero). For diffusion, the standard units are solvent flux units for (value) and alt time units for (timeOrRes). For poroelasticity, the standard units are pore volume fraction/((length units)²-time units) for (value) and alt time units for (timeOrRes). (Note: these units may change depending on the (style) setting).
• (function) - if the (style) is function (or 6), this attribute in XML files specifies a user defined function. In scripted files, the user defined function should replace the (value) argument and must be entered as quoted text. For diffusion, the function should evaluate to the desired concentration flux in solvent flux units. For poroelasticity, the function should evaluate to the desired flux in pore volume fraction/((length units)²-time units).

Note that valid calculation of total flux resulting from specified flux during 2D simulations requires that you have correctly specified the grid thickness and that the grid thickness matches the particle thickness that is set in Region commands. The way that pore volume fraction flux translates to pore pressure increase depends on the poroelasticity properties of materials and is desribed here.

Coupled Concentration Flux

This concentration flux mode calculates the concentration flux from a user-defined function. At each time step, the variable t will be replaced by the difference between the particle concentration potential (0 to 1) and the entered reservoir concentration potential (0 to 1) in parameter (timeOrRes). The function should evaluate to the desired flux in solvent flux units. The function can additionally depend on particle position (x, y, z, D, T, R, or Z in length units), and/or clockwise particle rotation (q in radians). Note that because t is used for particle concentration, a coupled surface flux function cannot depend on time.

Coupled Pore Pressure Flux

This pore pressure flux mode calculates the pore pressure flux from a user-defined function. At each time step, the variable t will be replaced by the difference between the particle pore pressure and the entered reservoir pore pressure in parameter (timeOrRes). The function should evaluate to the desired flux in pore volume fraction/((length units)²-time units). The function can additionally depend on particle position (x, y, z, D, T, R, or Z in length units), and/or clockwise particle rotation (q in radians). Note that because t is used for particle pore pressure, a coupled surface flux function cannot depend on time.

Transport Flux Conditions

When doing other transport calculations, you can set transport fluxes to surfaces of particles within the parents's block shape. The scripted ConcentrationFlux command selects transport fluxes by adding a (phaseStyle) command:

ConcentrationFlux (mode),(face),(style),(value),(timeOrRes),(phaseStyle)

In XML files, the <ConcFluxBC> command adds a phase attribute:

<ConcFluxBC dir='(mode)' face='(face)' style='(style)' value='(value)'
                        time='(timeOrRes)' function='(function)' phase='(phaseStyle)'/>

where (phaseStyle) can be the following options:

• fracture (phaseStyle=3)
• battery (phaseStyle=4)
• conduction (phaseStyle=5)

All other parameters are defined above. To specify (phaseStyle) in scripted files, the (time) parameters is required. It can be set to default zero if not needed.

Some of these options are available only in OSParticulas and only useful when using materials that support them.

Initial Particle Damage

The scripted Damage command adds initial damage state to particles within the parents's block shape. This option only applies if those particles are for a softening material. The format is:

Damage (nx),(ny),<(nz)>,<(dn)>,<(dxy)>,<(dxz)>,<(mode)>

In XML files, the command is:

<Damage nx='0' ny='2' dn='1' dxy='1'/>

where

• (nx, ny, nz) is normal vector for the damage plane. It need not be normalized, but should not be a zero vector (nz only used for 3D).
• (dn, dxy, dxz) defines the three damage parameters (dxz only used for 3D). These should all be from 0 to 1 where 0 is undamaged and 1 is fully damaged. If not specified, the default values are 1 (or fully damaged)
• (mode) damage mode or value for material type in the first damage history variable. This value is only used if all d's are less than 1. The default is 1 for IsoSoftening, but must be set for TransIsoSoftening materials.

For 2D problems, the Damage command must still use up to 7 parameters (for alignment). The two properties entered for 3D problems will be ignoed. In XML files, specify all the parameters that are needed.

Comments on Initial Damage

As a general rule, an thickness of the initial damage zone should span at least one full cell. A single line of particles when using cells with two particles in each direction, will not create the expected damage plane. Instead, the plane will be constrained by other material points connected as a continuum. That line will cause other particles in the same cell to get over stressed and they will eventually become damaged as well.

If damage parameters are equal to 1 (which is most common), the particles are fully damaged and first softening history variable is set to 3 (so visualization can distinguish new damage from predamage). When all damage variables are less than one, the code will need to find maximum strain corresponding to that damage state and first softening history variable is set to provided (mode). See each material type for failure mode settings.

The setting of initial particle damage is not yet implemented for TransIsoSoftening materials.

Initial Particle Phase Field

The scripted PhaseField command adds initial phase field value to particles within the parents's block shape. This option only applies if those particles for a material that uses a phase field (current only for Isotropic Phase Field Softening Material). The format is:

PhaseField (phi)

In XML files, the command is:

<Damage phi='(phi)'/>

where (phi) is a user defined function for the phase field value to set on the particle. The function should evaluate to between 0 and 1 (and values outside that range are truncated to be between 0 and 1).