Follower pressure loads: theory

A pressure created by a gas or a liquid acts normal to the current surface and on its current area. With Shell.add_pressure_load() and follower=True the pressure \(p\) follows the deformed mid-surface of the shell, as opposed to the default dead load, follower=False, which keeps the direction and the area of the undeformed mid-surface. The pressure \(p(x, y)\) is a function of the undeformed coordinates, i.e. a body-attached follower load in the classification of Schweizerhof and Ramm [schweizerhof1984Pressure]. It is the load of the hydrostatic buckling of rings and cylinders, see e.g. [timoshenko1961Stability], [brush1975Buckling] and [simitses1976Stability], and of the elasticity solutions of thick orthotropic and sandwich cylinders under external pressure [kardomateas1993Orthotropic], [kardomateas1995Orthotropic], [kardomateas2003Sandwich], [han2004SandwichCylinder], [kardomateas2005Sandwich], in which the pressure “remains normal to the deflected surface during the buckling process”.

This section derives the force vector of the follower pressure, its contribution to the internal force vector and its load stiffness matrix Shell.calc_kCfollower(), and discusses when this matrix is symmetric. The derivation is implemented with sympy in theory/shells/follower_pressure/follower_pressure.py, which also generates the kernels calc_fext_follower and fkCfollower_num of every model, see Implementation and verification.

Virtual work of the pressure

Let \(x\) and \(y\) be the coordinates of the undeformed mid-surface, \(y\) being the arc length for the cylindrical shells, and \(\{e_x, e_y, e_z\}\) the local basis of the undeformed mid-surface, with \(e_z\) the normal, positive outwards for the cylinders of radius \(r\), such that \(e_{y,y} = -e_z/r\) and \(e_{z,y} = e_y/r\). The displacement of the mid-surface is \(\mathbf{u} = u e_x + v e_y + w e_z\), and the position of the deformed mid-surface is \(\mathbf{X} = \mathbf{X}_0 + \mathbf{u}\). A pressure \(p\), positive along \(+e_z\) as the force fz of Shell.add_point_load() (an external pressure on a cylinder is negative), acting normal to the deformed mid-surface and on its deformed area, does the virtual work:

(1)\[\delta W_p = \int_{A} p \, \mathbf{n} \cdot \delta\mathbf{u} \, dA = \int_{A_0} p \, \mathbf{a} \cdot \delta\mathbf{u} \, dx \, dy, \qquad \mathbf{a} = \mathbf{X}_{,x} \times \mathbf{X}_{,y}\]

where \(\mathbf{a}\) is the area vector of the deformed mid-surface per unit undeformed area, \(\mathbf{n} \, dA = \mathbf{a} \, dx \, dy\). For the models based on the shear deformation theories, FSDT and TSDT, the pressure acts on the mid-surface, such that only \(u\), \(v\) and \(w\) enter Eq. (1), never the rotations \(\phi_x\) and \(\phi_y\).

In the local basis:

(2)\[\begin{split}\mathbf{X}_{,x} = \mathbf{A} = \begin{Bmatrix} 1 + u_{,x} \\ v_{,x} \\ w_{,x} \end{Bmatrix}, \qquad \mathbf{X}_{,y} = \mathbf{B} = \begin{Bmatrix} u_{,y} \\ 1 + v_{,y} + k_w w \\ w_{,y} - k_v v \end{Bmatrix}\end{split}\]

with the factors of the kinematics of each model:

kinematics	\(k_w\)	\(k_v\)
plate	\(0\)	\(0\)
cylinder, Donnell	\(1/r\)	\(0\)
cylinder, Sanders	\(1/r\)	\(1/r\)

For the plates and the Sanders cylinders Eq. (2) is the exact geometry, and \(w_{,y} - v/r\) is the rotation of the Sanders kinematics, \(\phi_y = -w_{,y} + v/r\), see theory/shells/cylshell_clpt_sanders and theory/shells/fsdt_tsdt. The Donnell cylinders neglect \(v/r\) in the rotation, as their strains do. The area vector \(\mathbf{a} = \mathbf{A} \times \mathbf{B}\) is bilinear in the displacement gradients:

(3)\[\begin{split}\mathbf{a} = \mathbf{a}_0 + \mathbf{a}_1 + \mathbf{a}_2, \qquad \mathbf{a}_0 = \begin{Bmatrix} 0 \\ 0 \\ 1 \end{Bmatrix}, \qquad \mathbf{a}_1 = \begin{Bmatrix} -w_{,x} \\ -w_{,y} + k_v v \\ u_{,x} + v_{,y} + k_w w \end{Bmatrix},\end{split}\]

(4)\[\begin{split}\mathbf{a}_2 = \begin{Bmatrix} v_{,x} (w_{,y} - k_v v) - w_{,x} (v_{,y} + k_w w) \\ w_{,x} u_{,y} - u_{,x} (w_{,y} - k_v v) \\ u_{,x} (v_{,y} + k_w w) - v_{,x} u_{,y} \end{Bmatrix}\end{split}\]

The dead part \(\mathbf{a}_0\) is the dead load. The first-order part \(\mathbf{a}_1\) contains the rotation of the normal, \((-w_{,x}, -w_{,y} + k_v v)\), and the change of area, \(u_{,x} + v_{,y} + k_w w\), which is the sum of the linear membrane strains \(\varepsilon_{xx} + \varepsilon_{yy}\).

Truncation consistent with the strains

The strains of all models are of moderate rotations, von Kármán, Donnell or Sanders, i.e. with the in-plane gradients \(u_{,x} \sim \varepsilon\) and the rotations \(w_{,x} \sim \varepsilon^{1/2}\). The virtual displacements are ordered alike, \(\delta u \sim \varepsilon^{1/2} \delta w\), such that the integrand of Eq. (1), \(p (a_z \delta w + a_x \delta u + a_y \delta v)\), has:

• from \(\mathbf{a}_0 + \mathbf{a}_1\): the terms \(p \, \delta w\) and \(p \, O(\varepsilon) \, \delta w\);

• from \(\mathbf{a}_2\): the terms \(p \, O(\varepsilon^2) \, \delta w\), of the order of the terms neglected by the strains themselves.

The consistent truncation is therefore the first-order area vector, follower=True or 'linear':

(5)\[\mathbf{a} \approx \mathbf{a}_0 + \mathbf{a}_1\]

With follower='quadratic' the complete area vector of Eq. (3) is used instead, which is useful to check the sensitivity of a result to the truncation.

Force vector

With the approximation of the displacements by the Ritz constants \(\{c\}\), see Bardell’s Hierarchical Functions,

(6)\[u = \{N_u\}^T \{c\}, \qquad v = \{N_v\}^T \{c\}, \qquad w = \{N_w\}^T \{c\}\]

where \(\{N_u\}\), \(\{N_v\}\) and \(\{N_w\}\) are the vectors of the approximation functions of each field, null at the positions of the other degrees of freedom, Eq. (1) gives the force vector of the follower pressure:

(7)\[\{F_p(c)\} = \int_{A_0} p \left( \{N_u\} a_x + \{N_v\} a_y + \{N_w\} a_z \right) dx \, dy\]

where \(\mathbf{a}\) is evaluated at \(\{c\}\) with Eq. (5), or Eq. (3) for 'quadratic'. In the undeformed state, \(\mathbf{a} = \mathbf{a}_0\) and the follower load equals the dead load:

(8)\[\{F_p(0)\} = \int_{A_0} p \{N_w\} \, dx \, dy\]

which is the reference force vector returned by Shell.calc_fext(), see shell_fext(). For the first-order area vector \(\{F_p\}\) is affine in \(\{c\}\):

(9)\[\{F_p(c)\} = \{F_p(0)\} + [K_p] \{c\}\]

Load stiffness matrix

The Jacobian of the force vector, the load stiffness, is:

(10)\[[K_p] = \frac{\partial \{F_p\}}{\partial \{c\}}, \qquad K_{p_{ij}} = \int_{A_0} p \, \mathbf{N}_i \cdot \frac{\partial \mathbf{a}}{\partial c_j} \, dx \, dy\]

with \(\mathbf{N}_i = (N_{u_i}, N_{v_i}, N_{w_i})\) and, from \(\mathbf{a} = \mathbf{A} \times \mathbf{B}\):

(11)\[\frac{\partial \mathbf{a}}{\partial c_j} = \frac{\partial \mathbf{A}}{\partial c_j} \times \mathbf{B} + \mathbf{A} \times \frac{\partial \mathbf{B}}{\partial c_j}\]

where \(\partial \mathbf{A}/\partial c_j = (N_{u_j,x}, N_{v_j,x}, N_{w_j,x})\) and \(\partial \mathbf{B}/\partial c_j = (N_{u_j,y}, N_{v_j,y} + k_w N_{w_j}, N_{w_j,y} - k_v N_{v_j})\). For the first-order area vector, \(\mathbf{A}\) and \(\mathbf{B}\) are replaced by \(\mathbf{A}_0 = e_x\) and \(\mathbf{B}_0 = e_y\) in Eq. (11), \([K_p]\) does not depend on \(\{c\}\), and its blocks between the fields of the rows (\(\delta u, \delta v, \delta w\)) and of the columns (\(u, v, w\)) are:

(12)\[\begin{split}[K_p] = \int_{A_0} p \begin{bmatrix} 0 & 0 & -\{N_u\}\{N_{w,x}\}^T \\ 0 & k_v \{N_v\}\{N_v\}^T & -\{N_v\}\{N_{w,y}\}^T \\ \{N_w\}\{N_{u,x}\}^T & \{N_w\}\{N_{v,y}\}^T & k_w \{N_w\}\{N_w\}^T \end{bmatrix} dx \, dy\end{split}\]

The matrix returned by Shell.calc_kCfollower() is the contribution to the tangent stiffness matrix:

(13)\[[K_{C_{follower}}] = - \frac{\partial \{F_p\}}{\partial \{c\}} = -[K_p]\]

such that all the stiffness matrices are additive, see Eq. (16). For an external pressure on a cylinder with Sanders’ kinematics, \(p < 0\), the block \(-k_v p \{N_v\}\{N_v\}^T\) of \([K_{C_{follower}}]\) is positive definite, the stiffening that cancels the destabilizing hoop prestress of the rigid rotation, see Ring under uniform pressure.

Internal force vector, residual and tangent stiffness

The pressures of the loads added with cte=False are multiplied by the load factor \(\lambda\), inc, and those with cte=True are constant. Writing \(\{F_p(c)\}\) for the sum over the follower loads with their multipliers, the residual of the equilibrium is:

(14)\[\{R(c, \lambda)\} = \{F_{ext}(\lambda)\} + \left( \{F_p(c)\} - \{F_p(0)\} \right) - \{F_{int}^e(c)\}\]

where \(\{F_{ext}(\lambda)\}\) is Shell.calc_fext(), which includes the reference part \(\{F_p(0)\}\) of the follower loads, and \(\{F_{int}^e(c)\}\) is the internal force vector of the strains. Shell.calc_fint() returns the internal force vector with the configuration-dependent part of the follower loads:

(15)\[\{F_{int}(c, \lambda)\} = \{F_{int}^e(c)\} - \left( \{F_p(c)\} - \{F_p(0)\} \right)\]

such that \(\{R\} = \{F_{ext}(\lambda)\} - \{F_{int}(c, \lambda)\}\), the form used by structsolve.Analysis. The tangent stiffness matrix, Shell.calc_kT(), is its exact Jacobian:

(16)\[[K_T] = \frac{\partial \{F_{int}\}}{\partial \{c\}} = [K_0] + [K_{0L}] + [K_{L0}] + [K_{LL}] + [K_{G_{NL}}] + [K_G(N_0 + N_L)] + [K_{C_{follower}}]\]

and the derivative of the residual with respect to the load factor, used by the arc-length methods, is the load vector of the current configuration, calc_fext(inc=1., c=c) for the loads proportional to \(\lambda\):

(17)\[\frac{\partial \{R\}}{\partial \lambda} = \{F_{ext}(1)\} + \{F_p(c)\} - \{F_p(0)\}\]

From version 0.5.0, structsolve passes the load factor to the callables that accept the keyword argument inc and uses Eq. (17) in the arc-length methods.

Linear analyses

Linear static analysis. A geometrically linear structure under a follower load with the first-order area vector is in equilibrium when \([K_0]\{c\} = \{F_p(c)\}\), i.e., from Eqs. (9) and (13):

(18)\[\left( [K_0] + [K_{C_{follower}}] \right) \{c\} = \{F_p(0)\}\]

an unsymmetric linear system in general, solved by structsolve.Analysis.static(NLgeom=False) when calc_fext accepts the configuration c.

Linear buckling, static criterion. With the prestate \(\{c_0\}\) of the reference load, the geometric stiffness \([K_G(c_0)]\) and the load stiffness both scale with the load factor:

(19)\[\left( [K_0] + \lambda \left( [K_G(c_0)] + [K_{C_{follower}}] \right) \right) \{c\} = \{0\}\]

solved by structsolve.lb(kC0, kG + kCfollower), which uses the solvers of unsymmetric matrices when needed.

Kinetic criterion. A non-conservative system may lose stability by flutter, which Eq. (19) does not detect: the eigenvalues \(\lambda^2 = -\omega^2\) of \(([K_T] + \lambda^2 [M])\{c\} = \{0\}\), with the load stiffness in \([K_T]\), are real below the critical load and become complex conjugate pairs at a flutter load, structsolve.freq.

Symmetry of the load stiffness matrix

For the first-order area vector, \([K_p]\) is the matrix of the bilinear form

(20)\[P(\delta, \Delta) = \int_{A_0} p \, \delta\mathbf{u} \cdot \mathbf{a}_1(\Delta\mathbf{u}) \, dx \, dy\]

and integrating by parts, the skew part is (verified symbolically by check_symmetry of follower_pressure.py):

(21)\[P(\delta, \Delta) - P(\Delta, \delta) = - \oint_{\partial A_p} p \left[ (\delta u \, \Delta w - \delta w \, \Delta u) n_x + (\delta v \, \Delta w - \delta w \, \Delta v) n_y \right] ds + \int_{A_p} \left[ p_{,x} (\delta u \, \Delta w - \delta w \, \Delta u) + p_{,y} (\delta v \, \Delta w - \delta w \, \Delta v) \right] dx \, dy\]

over the loaded region \(A_p\), independently of \(k_w\) and \(k_v\). Therefore \([K_{C_{follower}}]\) is symmetric if and only if \(p\) is uniform over the whole integration domain and, on each edge, either \(w = 0\) or the displacement normal to the edge is zero, \(u\) on the edges \(x = const\) and \(v\) on the edges \(y = const\). This is the case of simply supported or clamped panels, of symmetry planes and of closed rings. Free or partially restrained edges, patches, whose edges are inside the shell, and a non-uniform \(p(x, y)\) give an unsymmetric load stiffness, a non-conservative load. This is Table 2 of [schweizerhof1984Pressure] for body-attached loads, see also [hibbitt1979Follower]. When it is symmetric the load is conservative, with the potential

(22)\[\Pi_p = - \lambda \left( \{F_p(0)\}^T \{c\} + \frac{1}{2} \{c\}^T [K_p] \{c\} \right)\]

For the quadratic area vector \([K_p]\) depends on \(\{c\}\) and is not symmetric in general. A hydrostatic head, a pressure that depends on the current position, would add a further term, not implemented.

Load surface

A pressure \(p\) on the surface \(z = z_p\) of a cylinder, e.g. the outer face \(z_p = h/2\) of a thick sandwich, has the resultant \(p (r + z_p)\) per unit length of mid-surface arc. With add_pressure_load(..., zp=zp) it is applied as

(23)\[p \left( 1 + \frac{z_p}{r} \right)\]

per unit mid-surface area, with the direction and the change of area of the mid-surface, as in a shallow shell, see pressure_factor(). The terms of relative order \(z_p \kappa\), \(\kappa\) the change of curvature, in the change of area and the moment of the in-plane components about the mid-surface are neglected; for a ring mode \(n\) they are of relative order \(n^2 z_p / r\) in the load stiffness. Plates are not affected.

Ring under uniform pressure

For a ring in plane strain, mode \(w = W \cos(n y/r)\), \(v = V \sin(n y/r)\), membrane prestress \(N_{yy} = p r\) and the CLPT strains of the kernels, the second variation of the total potential is

(24)\[U_2 = \frac{1}{2} \oint \left( A \varepsilon_{yy}^2 + D \kappa_{yy}^2 + N_{yy} \beta_y^2 \right) dy - \frac{1}{2} \oint p \, \mathbf{u} \cdot \mathbf{a}_1(\mathbf{u}) \, dy\]

with \(\beta_y = w_{,y} - k_v v\). In the inextensional limit its singularity gives \(p_{cr} r^3/D\) (ring_buckling in follower_pressure.py):

kinematics	follower (hydrostatic)	dead (constant direction)	centrally directed
Sanders	\(-3\)	\(-4\)	\(-9/2\)
Donnell	\(-16/5\)	\(-4\)	\(-64/15\)

The Sanders values are the classical ones, see [timoshenko1961Stability], [brush1975Buckling] and Table 3 of [schweizerhof1984Pressure], and the follower value is the hydrostatic critical pressure \(3D/r^3\) of the long cylinders, Eq. (66) of [nasa2020SP8007]. With the FSDT and the transverse shear stiffness \(S\), the follower load gives \(p_{cr} = 3D/(r^3 (1 + 4D/(S r^2)))\), the shell formula of [han2004SandwichCylinder].

With Sanders’ kinematics and a dead pressure, the rigid rotation \(v = const\), \(w = 0\), strain free, has the energy \(\oint N_{yy} (v/r)^2 dy < 0\) under an external pressure, a spurious mode at a vanishing critical load whenever \(v\) is free at an edge. The follower load stiffness adds \(-\oint p \, v^2/r \, dy\), which cancels it exactly.

Implementation and verification

• theory/shells/follower_pressure/follower_pressure.py: the sympy derivation of Eqs. (3), (7) and (10) in the notation of the kernels, the symmetry identity Eq. (21) and the ring closed forms; write_pyx.py writes the kernels calc_fext_follower (Eq. (7)) and fkCfollower_num (Eq. (13), assembled in full, not as an upper triangle) into every panels/models/<model>_num.pyx, integrated with nx*ny Gauss-Legendre points over each loaded patch, the same points as Eq. (8) in shell_fext().

• tests/tests_shell/test_follower_pressure.py: Eq. (8), Taylor tests of Eq. (7) against Eq. (10), finite-difference Jacobian of Eq. (15) against Eq. (16), quadratic convergence of Newton-Raphson and Riks, the symmetry conditions of Eq. (21) and Eq. (18).

• tests/tests_shell/test_follower_pressure_validation.py and the notebooks of Follower (hydrostatic) pressure: the ring closed forms, finite cylinders, [nasa2020SP8007], [kardomateas1993Orthotropic], [han2004SandwichCylinder], [kardomateas2003Sandwich] and [schweizerhof1984Pressure].