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='(fact)' 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. Furthermore, the use of traction boundary conditions is only allowed when using uGIMP or lCPDI anlaysis methods.

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, which depending on the current MPM method being used, is as follows:

• Dirac: Traction conditions are not allowed in this method; only Load conditions can be used
• uGIMP: Integrate on the surface determined from the initial particle shape but translated to the new particle position.
• lCPDI or qCPDI: Find parallelogram (2D) or parallelopiped (3D) for the deformed particle domain from the particle's deformation gradient and integrate the selected surface terms by expanding its integrand in linear shape functions (i.e., very similar to CPDI methods for shape function integration).

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, which means its effect depends on MPM method being used. When using uGIMP, the deformed particle domain is the same as the initial particle domain and therefore the applied stress will always be an engineering stress based on initial particle area (i.e., the stress will change if particle cross-sectional area changes). But, when using lCPDI or qCPDI, 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. Furthermore, the use of heat flux boundary conditions is only allowed when using uGIMP or lCPDI anlaysis methods.

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 Flux Conditions

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

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

In XML files, the command is:

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

where

• (mode) determines the type of concentration flux. The two options are:
  1. external - This condition applies a surface flux in solvent flux units for transport rate of substance 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 concentration flux from a function of the particle concentration potential.
• (face) - specifies the particle surface that is on the edge and will have the concentration 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 flux. For coupled concentration flux, the (style) must be function (or 6) and the variable (timeOrRes) should specify the solvent reservoir potential (0 to 1).
• (value) and (timeOrRes) - for external heat flux, the value depends on two parameters specified by arguments (value) and (timeOrRes) (if omitted, their values are set to zero). The standard units are solvent flux units for (value) and alt time units for (timeOrRes) (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 concentration flux in solvent flux units.

Note that valid calculation of total flux resulting from specified concentration 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. Furthermore, the use of heat flux boundary conditions is only allowed when using uGIMP or lCPDI anlaysis methods.

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.

Initial Particle Damage

The scripted Damage command adds initial damage state to particles within the parents's block shape. The 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 enter for 3D problems will be ignoed. In XML files, specify all the parameters that are needed.

Comments on Initial Damage

As a general rule, therefore, the 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 you implemented for TransIsoSoftening 1 or 2 materials.