calculate stress command

Purpose

Calculates the stress of a collection of particles

Synatx

calculate stress keyword/value pairs

The following keyword/value pairs are allowed:

Keywords	Description
id	user-assigned name for the calculate stress command
region	region in which the stress tensor is calculated
particle_group	particle group for which the stress tensor is calculated
kernel_radius	radius of the smoothing kernel default: 4*maximum radius; range: (0,∞); units: [length]
kernel_type	type of the smoothing kernel, must be one of: `Gaussian`, `Top_Hat`, `Wendland` default: Gaussian;

Examples

calculate stress id cs
calculate stress id cs kernel_radius 1e-2 kernel_type Wendland region outlet
# calculate average xx stress:
calculate average id cs_xx quantity id_stressTensor_cs.xx

Description

The stress tensor is calculated for each particle according to Goldhirsch. The formula is given by:

$\sigma_{i,ab} = -1/2 \sum_{j,k} f_{jk,a} r_{jk,b} \int_0^1 \phi(r_i - r_j + s r_{jk}) ds - \sum_j m_j v'_{ij,a} v'_{ij,b} \phi(r_i - r_j) + \sigma_{wall,i,ab}$

where

$v'_{ij} = v_j - V_i V_i = p_i / rho_i p_i = \sum_j m_j v_j \phi(r_i - r_j) \rho_i = \sum_j m_j \phi(r_i - r_j) \phi(r) = H(radius - |r|) / \Omega(radius) \Omega(radius) = 4/3 \pi radius^3$

and

 is the velocity vector
 force acting from j onto i
 vector from center of j to center of i
 mass of particle i
 is the Heaviside function

In case solid boundaries are present the last term is given according to Weinhart et al. by

$\sigma_{wall,i,ab} = - \sum_{j,k} f_{jk,a} a_{jk,b} \int_0^1 \phi(r_i - r_j + s a_{jk}) ds$

where

$a_{jk} = r_j - c_{jk}$

and $c_{jk}$ is the contact point of particle j with wall k and the sum runs over all particles j and walls k.

The following three kernel_type values are implemented at the moment:

Top_Hat - Top hat kernel

$w(r) = a_t \quad if \, q < 1 \quad (q = r / kernel\_radius) = 0 \quad otherwise$

Gaussian - Gaussian kernel

$w(r) = a_g exp(-q^2) \quad (q = 3 r / kernel\_radius)$

Wendland - Quintic radial polynomial

$w(r) = a_w (1-q/2)^4 (1+2q) \quad (q = 2 r / kernel\_radius)$

Note that all kernels are equal to zero if $r$ is greater than kernel_radius (this implies a cut-off for the Gaussian). The constants $a$ (different for each kernel) are chosen such that the integral of $w$ over the ball of radius kernel_radius is equal to one. In case of the top hat kernel $a_t$ is equal to the volume of this sphere.

The command can be restricted to work on a specific group of particles or a specific region by using the particle_group or region keywords, respectively.

The stress tensor is written to a per-particle property which can be accessed using id_stressTensor_ID, where ID is the id of the calculate command. It has nine components which can be accessed using .xx, .xy, .xz, .yx, .yy, .yz, .zx, .zy and .zz for the respective components.

References

(Goldhirsch) Goldhirsch; Stress, stress asymmetry and couple stress: from discrete particles to continuous fields, Granular Matter (2010)

(Weinhart) Weinhart, Thornton, Luding, Bokhove; From discrete particles to continuum fields near a boundary (2012)