Release nodal force boundary condition

This MR presents a new boundary condition for mechanical process: ReleaseNodalForce, particularly for excavation modelling.

This boundary condition applies a time-dependent release nodal force to nodes on an exposed surface to simulate excavation processes.

The initial state assumes a non-equilibrium stress denoted by \sigma_0. In the finite element method (FEM), the nodal force is given by:

\mathbf{b} = \int \left( \mathbf{B}^\top (\boldsymbol{\sigma} - \boldsymbol{\sigma}_0) + (\mathbf{f} - \mathbf{f}_0) \right) N \, \mathrm{d}\Omega + \int_{\Gamma_q} (\boldsymbol{\sigma}-\boldsymbol{\sigma}_0)\cdot \mathbf n \mathrm{d}\Gamma

Where:

• \mathbf{B} is the strain-displacement matrix,
• \boldsymbol{\sigma} is the current total stress,
• \boldsymbol{\sigma}_0 is the initial total stress,
• \mathbf{f} is the current body force,
• \mathbf{f}_0 is the initial body force,
• \int_{\Gamma_q} is the boundary where the traction condition is applied, and \mathbf n is the outward normal of the boundary,
• N is the shape function,
• \Omega is the domain of integration.

After excavation, the stress and body force inside the excavated domain vanish.
This leads to non-zero nodal forces at the exposed surface nodes, computed as:

\mathbf{b}_0 = -\int \left( \mathbf{B}^\top \boldsymbol{\sigma}_0 + \mathbf{f}_0 N \right) \, \mathrm{d}\Omega - \int_{\Gamma_q} \boldsymbol{\sigma}_0 \cdot \mathbf n \mathrm{d}\Gamma

The components of \mathbf{b}_0 corresponding to the exposed surface nodes define the release nodal force vector:

\mathbf{f}_\text{r} := (\mathbf{b}_0)_i, \quad i \in \text{exposed surface nodes}

Note: If \Omega is the excavated domain, this leads to a negative \mathbf{f}_\text{r}.

To simulate gradual excavation, the release nodal force vector is applied to the global right-hand side (RHS) vector \mathbf{b}, scaled by a time- and position-dependent release parameter g(t, \mathbf{x}):

\mathbf{b} = \mathbf{b} + \mathbf{f}_\text{r} \cdot g(t, \mathbf{x})

Parameter requirements:

• Monotonically decreasing: \frac{\partial g}{\partial t} < 0
• g(0, \mathbf{x}) = 1
• g(t_e, \mathbf{x}) = 0, where t_e is the end time of excavation

1. Feature description was added to the changelog
2. Tests covering your feature were added?

Kirsch's problem, a classic example, is used to verify the implementation of the ReleaseNodalForce boundary condition.

Kirsch's problem defines an infinite elastic plate with a circular hole of radius a subjected to uniaxial tension \sigma. If the tension is along y-axis,

the stress solutions of this problem in polar coordinates (r, \theta) are given by [1]:

\begin{align} \sigma_r &= \frac{\sigma}{2} \left( 1 - \frac{a^2}{r^2} \right) - \frac{\sigma}{2} \left( 1 - \frac{4a^2}{r^2} + \frac{3a^4}{r^4} \right) \cos 2\theta \\ \sigma_\theta &= \frac{\sigma}{2} \left( 1 + \frac{a^2}{r^2} \right) + \frac{\sigma}{2} \left( 1 + \frac{3a^4}{r^4} \right) \cos 2\theta \\ \tau_{r\theta} &= \frac{\sigma}{2} \left( 1 + \frac{2a^2}{r^2} - \frac{3a^4}{r^4} \right) \sin 2\theta \end{align}

• At the hole's boundary r = a: \sigma_r = 0, \tau_{r\theta} = 0.
• At \theta = 0, \pi: \sigma_\theta = 3\sigma (maximum tensile stress).
• At \theta = \pm \pi/2: \sigma_\theta = -\sigma (compressive stress).

For the finite element analysis, the geometrical domain is bounded to 140 m \times 140 m with a circular hole of radius 6.5 m located at the domain center. By leveraging the problem's symmetry, we use a quarter of the domain, resulting in the mesh shown below:

The problem is analyzed using the release nodal force approach and a standard approach. The initial and boundary conditions for both approaches are provided in the table below:

Item	Standard approach	Release Nodal Force	Unit
Initial stress	(0.0, 0.0, 0.0, 0.0)	(0.0, -20.0, 0.0, 0.0)	MPa
Bottom	u_y=0	u_y=0	m
Left	u_x=0	u_x=0	m
Top	\sigma_y=-20	\sigma_y=-20 or 0	MPa
Arc	Traction \mathbf\tau=\mathbf 0	ReleaseNodalForce	N/A

These conditions result in \sigma=-20 MPa and a=6.5 m for the analytical solutions.

Note: The release nodal force approach requires the compensation of non equilibrium initial residuum.

The two approaches give the identical stress solutions:

Note: The Jupytext benchmark script includes mesh generation. However, the mesh files are still attached to provide a complete test for the user.

Reference

[1]: Kirsch, E.G., 1898. Die Theorie der Elastizität und die Bedürfnisse der Festigkeitslehre. Zeitschrift des Vereines Dtsch. Ingenieure 42, 797–807.

3. Any new feature or behaviour change was documented?