fix continuum/weighted command

Warning

GPU support for this command has not been tested and may not work as expected.

Warning

This command is deprecated, please use calculate stress and/or the calculate strain commandsinstead

Syntax

fix ID group-ID continuum/weighted keyword value

ID, group-ID are documented in fix command
continuum/weighted = style name of this fix command
zero or one keyword/value pairs may be appended

keyword = {kernel_radius, kernel_type, compute}
  kernel_radius value = radius
    radius = Radius of the smoothing kernel
  kernel_type value = type
    type = Type of kernel {Top_Hat, Gaussian, Wendland}
  compute value = compute-type
    compute-type = Which tensor(s) to compute {stress, strain, stress_strain}

Examples

fix 1 all continuum/weighted kernel_radius 0.01 compute stress
fix 1 all continuum/weighted kernel_radius 0.2 kernel_type Wendland compute stress_strain

Description

If the compute keyword is set to either stress or stress_strain this fix calculates the complete stress tensor at 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.

If the compute keyword is set to either strain or stress_strain this fix calculates the incremental strain tensor at each particle according to Zhang et al. The formula is given by

$\epsilon_{i,ab} = \frac{1}{2 rho_i} \sum_{j,k} m_j m_k \phi(r_{ij}) dt (v_{jk,a} \nabla\phi(r_{ik},b) + v_{jk,b} \nabla\phi(r_{ik},a))$

where most of the variables are given as above and additionally

 is the a-th component of the velocity difference between i and j
 is the a-th component of the gradient of phi with respect to 
dt is the time-step size

The following three kernel types 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 - Quinitic 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 > 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.

Restart, fix_modify, output, run start/stop

No information about this fix is written to binary restart files. None of the fix_modify options are relevant to this fix.

The values can be dumped by using the f_stressTensor_[i] and/or f_strainTensor_[i] (0 <= i <= 8) values in dump commands

No parameter of this fix can be used with the start/stop keywords of the run command.

Restrictions

Strain computation does not work with the default TOP_HAT kernel as its derivative is zero.

In order to ensure that all particles in the kernel radius are considered make use of the neigh_modify command. In particular the contact_distance_factor which should be set such that

2*min(radius)*contact_distance_factor >= kernel_radius