Axisymmetry in FLAC2D

Introduction

Axisymmetric modeling is a simplification technique widely used in geotechnical engineering to analyze problems that are rotationally symmetric about a central axis. Instead of modeling a full three-dimensional domain, the problem is reduced to a two-dimensional representation in the radial-axial (\(r\)-\(z\) plane). This approach significantly reduces computational effort while still capturing the essential behavior of the system. A system is axisymmetric when its material properties, geometry, loading, and boundary conditions do not vary with the angular coordinate, \(\theta\). Common applications include analysis of tunnels and shafts, cylindircal boreholes, laboratory-scale tests, underground storage caverns, fluid flow involving injection/extraction wells, among others.

In FLAC2D, axisymmetric logic is enabled by issuing the command model configure axisymmetry. This command must be issued before any zones are created, as it sets the framework for how the model will be interpreted and analyzed.

Figure 1 illustrates the coordinate system and the stresses acting on an elemental volume in an axisymmetric model. The radial direction, \(r\), extends outward from the axis of symmetry, while the axial direction, \(z\), runs parallel to the axis. The circumferential direction, \(\theta\), is implicit in the formulation, with all variables assumed to be independent of this coordinate. As illustrated in Figure 1 below, axisymmetry causes the elemental volume to increase with radius because a unit radian slice in the circumferential directions expands proportionally with radius.

../../../../../_images/axisymmetry.png

Figure 1: Coordinate system and stresses acting on the elemental volume in an axisymmetric model.

General Comments

  • Axisymmetry modeling is supported in mechanical, fluid, and thermal analyses in FLAC2D.

  • \(x=0\) is the axis of symmetry.

  • \(\left( x,y,z \right)\) in cartesian coordinates correspond to \(\left( r,z,\theta \right)\) in cylindrical coordinates. For example, viewing the \(zz\) component of stress in plots is accomplished by plotting the \(\sigma_{yy}\) component of stress. See the section on coordinate system below for more details.

  • Zones must be created in the \(+r\)-\(z\) plane. No negative radial grid point coordinates are allowed. In other words, grid points created to the left of the axis of symmetry are not valid and have no physical meaning.

  • By default, the radial component of velocity is fixed at the axis of symmetry. Explicitly setting it to zero is not required.

  • Many axisymmetric calculations require dividing by the radial coordinate of the grid points. To avoid division by zero at the axis of symmetry, FLAC2D evaluates these calculations at a small distance away from the axis, which is proportional to the zone edge length. FLAC2D uses the following equation: \(r_{e}=0.1 L_{min}\) where \(r_{e}\) is the effective radius employed at the axis of symmetry and \(L_{min}\) is the minimum zone edge length.

  • The area of a zone is weighted by the radial component of its centroid, e.g., the total area of the zone is taken as the area of the zone in the \(r\)-\(z\) plane multiplied by the circumferential length \(r\). Grid point mass is also weighted by the radial coordinate.

Coordinate System Transformation

The transformation from cartesian to cylindrical coordinates is as follows:

  1. Scalar quantities: unchanged, \(s \rightarrow s\).

  2. Vector quantities: \(\left( v_x, v_y, v_z \right) \rightarrow \left( v_r, v_z, 0 \right)\).

  3. Tensor quantities: \(\left( \sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{xz} \right) \rightarrow \left( \sigma_{rr}, \sigma_{zz}, \sigma_{\theta\theta}, \sigma_{rz}, 0, 0 \right)\).

The coordinate transformation is important to keep in mind when setting directional material properties, applying boundary conditions, initializing stresses, plotting results, and using FISH. Internally, FLAC2D handles the necessary coordinate transformations, but users must be aware of the correct components to use in their input and interpretation of modeling results. In general, the commands are identical to cartesian coordinate system. An example of setting boundary conditions and using FISH to initialize stresses is shown below:

model new
model configure axisymmetry

zone create quad size 10 10

zone cmodel assign elastic
zone property density 2e3 bulk 8.33e8 shear 3.85e8

; Set axial velocity to zero at the top and bottom
zone face apply velocity-y 0 range position-y 0 ; axial is y-component
zone face apply velocity-y 0 range position-y 1 ; axial is y-component

; Set radial velocity to zero at the right edge
zone face apply velocity-x 0 range position-x 1 ; radial is x-component

; Tensor components with fish
fish operator inistress(zp)
	local K0  = 0.6
	local sig = zone.stress(zp)
	local sax = sig->yy            ; axial stress is yy-component
	zone.stress(zp)->xx = K0 * sax ; radial stress is xx-component
	zone.stress(zp)->zz = K0 * sax ; hoop stress is zz-component
end

[inistress(::zone.list)]

Mechanical Analysis

In axisymmetric deformation, the problem is formulated in cylindrical coordinates under the assumption that all field variables are independent of the circumferential direction and that the circumferential displacement is zero. As a result, the non-zero stress components are the radial, hoop (circumferential), axial, and radial–axial shear stress. The corresponding strains arise from radial and axial displacements only, with the hoop strain determined kinematically by radial expansion or contraction. Constitutive relations are then applied to link these stresses and strains through an appropriate material law, with the hoop component playing a central role due to geometric confinement. This framework allows inherently three‑dimensional stress states to be captured efficiently within a two‑dimensional formulation.

The strain tensor components in cylindrical coordinates under axisymmetry are:

(1)\[\begin{split} \begin{split} \varepsilon_{rr} &= \frac{\partial u_r}{\partial r} \\ \varepsilon_{zz} &= \frac{\partial u_z}{\partial z} \\ \varepsilon_{\theta\theta} &= \frac{u_r}{r} \\ \varepsilon_{rz} &= \frac{1}{2} \left( \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \right) \\ \varepsilon_{\theta z} &= \varepsilon_{r \theta} = 0 \end{split}\end{split}\]

where the two-dimensional displacement field is defined as \(u = \left( u_r, u_z \right)\). The stress tensor components are related to the strain components through the constitutive model, which can be any material behavior supported by FLAC2D. The stress tensor in axisymmetric form is:

(2)\[\begin{split}\sigma = \begin{bmatrix} \sigma_{rr} & \sigma_{rz} & 0 \\ \sigma_{rz} & \sigma_{zz} & 0 \\ 0 & 0 & \sigma_{\theta\theta} \end{bmatrix}\end{split}\]

Mechanical boundary conditions for axisymmetric models are applied in the same manner as in cartesian coordinates.

Fluid Analysis

In the context of fluid flow in porous media, axisymmetry is particularly useful for modeling problems where flow is radially distributed around a central feature such as an injection or extraction well. Under the axisymmetric assumption, all dependent variables (e.g., pressure, saturation) are functions of radial and axial coordinates only, and are independent of the circumferential coordinate. This reduces the governing equations to a two-dimensional form while still capturing the essential physics of the three-dimensional flow process.

For steady-state single-phase flow in a saturated porous medium, Darcy’s law in axisymmetric coordinates introduces a geometric term reflecting the increasing circumferential area with radius. The continuity equation is:

\[\frac{1}{r} \frac{\partial}{\partial r} \left( r q_r \right) + \frac{\partial q_z}{\partial z} = Q\]

where \(q_r\) and \(q_z\) are the radial and axial components of the Darcy velocity, respectively, and \(Q\) represents any sources or sinks. The presence of the \(1/r\) term accounts for the geometric spreading of flow as it moves radially outward from the axis of symmetry, which is a key feature distinguishing axisymmetric flow from Cartesian flow. When combined with Darcy’s law (\(q=-\lambda \nabla p\), the governing equation becomes:

\[\frac{1}{r} \frac{\partial}{\partial r} \left( r \lambda \frac{\partial p}{\partial r} \right) + \frac{\partial}{\partial z} \left( \lambda \frac{\partial p}{\partial z} \right) = Q\]

The above equation assumes isotropic fluid mobility (\(\lambda\)), but it can be extended to anisotropic conditions by treating the radial and axial mobilities separately. Axisymmetric modeling of fluid processes is also applicable to more complex problems involving partially saturated flow.

Dirichlet boundary conditions for axisymmetric fluid flow are applied in the same manner as in cartesian coordinates. Discharge and leakage boundary conditions, however, cannot be specified directly at the axis of symmetry (\(x=0\)). This is because a radius-weighted area must be associated with the fluid flux. For the bounday conditions to function properly, the location of it must be at a non-zero radial coordinate. In practice, these boundary conditions are often applicable to injection or extraction wells, and the FLAC2D grid can be constructed such that the minimum \(x\) coordinate is the well radius. Volumetric flow rates \(Q\) can be converted to a discharge fluid flux \(q\) with the following equation:

\[q = \frac{Q}{2 \pi r h}\]

where \(r\) is the radial coordinate where the boundary condition is applied and \(h\) is the axial height over which the discharge is applied. This conversion accounts for the circumferential area associated with the flow at that radius, ensuring that the specified flow rate is correctly represented in the axisymmetric model.

Thermal Analysis

Axisymmetric thermal analysis follows an analogous formulation to to fluid flow, with the governing equations taking the same form but involving different physical coefficients and state variables. Where fluid analyis is based on pressure-driven flow described by Darcy’s law, thermal analysis is governed by heat conduction and, where applicable, advection. As a result, the same geometric considerations, including the radial dependence and associated \(1/r\) term, apply to thermal problems. Boundary conditions are also implemented in an equivalent manner, with prescribed temperature and heat flux directly analogous to pressure and fluid discharge.

Axisymmetric Examples