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 kernels calc_fext_follower and fkCfollower_num of every model
are generated with SymPy from this derivation, 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
\(\boldsymbol{u} = u e_x + v e_y + w e_z\), and the position of the deformed
mid-surface is \(\boldsymbol{X} = \boldsymbol{X}_0 + \boldsymbol{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:
where \(\boldsymbol{a}\) is the area vector of the deformed mid-surface per unit undeformed area, \(\boldsymbol{n} \, dA = \boldsymbol{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:
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 panels.models. The Donnell cylinders
neglect \(v/r\) in the rotation, as their strains do. The area vector
\(\boldsymbol{a} = \boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\) is
bilinear in the displacement gradients:
The dead part \(\boldsymbol{a}_0\) is the dead load. The first-order part \(\boldsymbol{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 \(\boldsymbol{a}_0 + \boldsymbol{a}_1\): the terms \(p \, \delta w\) and \(p \, O(\varepsilon) \, \delta w\);
from \(\boldsymbol{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':
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 \(\boldsymbol{c}\), see Bardell’s hierarchical functions,
where \(\boldsymbol{S}^u\), \(\boldsymbol{S}^v\) and \(\boldsymbol{S}^w\) are the rows 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:
where \(\boldsymbol{a}\) is evaluated at \(\boldsymbol{c}\) with Eq.
(5), or Eq. (3) for 'quadratic'. In
the undeformed state, \(\boldsymbol{a} = \boldsymbol{a}_0\) and the follower
load equals the dead load:
which is the reference force vector returned by Shell.calc_fext(),
see shell_fext(). For the first-order area vector \(\boldsymbol{F}_p\) is affine
in \(\boldsymbol{c}\):
with the load stiffness matrix \(\boldsymbol{K}_F\) of Eq. (10).
Load stiffness matrix#
The load stiffness matrix, returned by Shell.calc_kCfollower(), is
the Jacobian of the force vector with the sign of a stiffness, such that
all the stiffness matrices are additive, see Eq. (15):
with \(\boldsymbol{S}_i = (S^u_i, S^v_i, S^w_i)\) the functions of the degree of freedom \(i\) and, from \(\boldsymbol{a} = \boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\):
where \(\partial \boldsymbol{X}_{,x}/\partial c_j = (S^u_{j,x}, S^v_{j,x}, S^w_{j,x})\) and \(\partial \boldsymbol{X}_{,y}/\partial c_j = (S^u_{j,y}, S^v_{j,y} + k_w S^w_j, S^w_{j,y} - k_v S^v_j)\). For the first-order area vector, \(\boldsymbol{X}_{,x}\) and \(\boldsymbol{X}_{,y}\) are replaced by \(e_x\) and \(e_y\) in Eq. (11), \(\boldsymbol{K}_F\) does not depend on \(\boldsymbol{c}\), and its blocks between the fields of the rows (\(\delta u, \delta v, \delta w\)) and of the columns (\(u, v, w\)) are:
For an external pressure on a cylinder with Sanders’ kinematics, \(p < 0\), the block \(-k_v p (\boldsymbol{S}^v)^\top \boldsymbol{S}^v\) of \(\boldsymbol{K}_F\) 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 \(\boldsymbol{F}_p(\boldsymbol{c})\) for the sum over the follower loads with their
multipliers, the residual of the equilibrium is:
where \(\boldsymbol{F}_{ext}(\lambda)\) is Shell.calc_fext(), which includes the
reference part \(\boldsymbol{F}_p(0)\) of the follower loads, and
\(\boldsymbol{F}_{int}^e(\boldsymbol{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:
such that \(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{c}, \lambda) - \boldsymbol{F}_{ext}(\lambda)\),
the form used by structsolve.Analysis. The tangent stiffness
matrix, Shell.calc_kT(), is its exact Jacobian, \(\partial \boldsymbol{R}/\partial \boldsymbol{c}\):
and minus 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\):
From version 0.5.0, structsolve passes the load factor to the callables
that accept the keyword argument inc and uses Eq. (16) 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 \(\boldsymbol{K}_{C0}\boldsymbol{c} = \boldsymbol{F}_p(\boldsymbol{c})\), i.e., from Eq. (9):
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 pre-buckling state \(\hat{\boldsymbol{c}}\) of the reference load, the geometric stiffness \(\boldsymbol{K}_G(\hat{\boldsymbol{c}})\) and the load stiffness both scale with the load factor:
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. (18) does not detect: the eigenvalues \(\lambda^2
= -\omega^2\) of \((\boldsymbol{K}_T + \lambda^2 \boldsymbol{M})\boldsymbol{c} = \boldsymbol{0}\), with the load
stiffness in \(\boldsymbol{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, the block of the rows \(\delta u\) and of the columns \(w\) of Eq. (12) is \((\boldsymbol{S}^u)^\top \boldsymbol{S}^w_{,x}\). A symmetric matrix would have its transpose, \((\boldsymbol{S}^w_{,x})^\top \boldsymbol{S}^u\), as the block of the rows \(\delta w\) and of the columns \(u\), which is instead \(-(\boldsymbol{S}^w)^\top \boldsymbol{S}^u_{,x}\), and similarly for the coupling of \(v\) and \(w\). The difference between the first block and the transpose of the second one is a derivative,
whereas the blocks of the curvature, \(-k_v (\boldsymbol{S}^v)^\top \boldsymbol{S}^v\) and \(-k_w (\boldsymbol{S}^w)^\top \boldsymbol{S}^w\), are symmetric. Integrating by parts over the loaded region \(A_p\), whose contour \(\partial A_p\) has the outward normal \(\{n_x, n_y\}\), the skew part of the load stiffness matrix of Eq. (12) is:
The curvature terms having cancelled, Eq. (20) holds for the plates and for the Donnell and Sanders cylinders alike. Hence, \(\boldsymbol{K}_F\) is symmetric, and the follower pressure conservative, if and only if:
\(p\) is uniform over the loaded region, such that the area integral vanishes; and
at each point of the contour of the loaded region, either \(w = 0\) or the in-plane displacement normal to the contour is zero, i.e. \(u = 0\) on the edges \(x = const\) and \(v = 0\) on the edges \(y = const\), such that the contour integral vanishes.
For a load over the whole domain, the second condition is met on the edge
\(x = 0\) with the flag x1w = 0 or x1u = 0 of Shell, see
Bardell’s hierarchical functions, and likewise on the other edges. These
conditions are met by simply supported and clamped edges, by symmetry
planes, and by closed surfaces such as complete rings and cylinders, for
which no boundary exists; in these cases the load potential is the pressure
times the change of the enclosed volume. Conversely, \(\boldsymbol{K}_F\)
is unsymmetric for free or partially restrained edges, for a pressure
applied on a patch whose edges lie inside the shell, and for a pressure that
varies over the surface. This agrees with [hibbitt1979Follower], p. 938,
for whom a uniform pressure on a fully enclosed volume, or on a surface
whose boundaries have prescribed displacements, gives a symmetric load
stiffness, unlike free edges and a non-uniform pressure, p. 939, and, for
body-attached loads, with Table 2 of [schweizerhof1984Pressure], whose
boundary terms vanish under the conditions of their Sec. 3.2.1, Eqs.
(32)-(35), e.g. when at least two displacement components are prescribed
at the edge, or when the edge is restrained normal to a continuous
surface, such as the undeformed surface of the shell, \(w = 0\).
A closed cylinder of panels is an assembly of domains, see
create_cylinder(), whose load stiffness matrices are each
unsymmetric: the contour integrals along a connected edge cancel only when
the displacements of the two sides are equal, which holds for the matrices
reduced by the null-space method, symmetric to round-off, but not for the
assembled matrix of the penalty method, which remains unsymmetric. When it
is symmetric the load is conservative, with the potential
For the quadratic area vector \(\boldsymbol{K}_F\) depends on \(\boldsymbol{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
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
with \(\beta_y = -w_{,y} + k_v v\). In the inextensional limit its singularity gives \(p_{cr} r^3/D\), derived with SymPy:
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] and [brush1975Buckling]. For a ring, with \(D\) replaced by the bending stiffness \(EI\) of the ring and the extension of the ring axis neglected, they are the critical loads \(3EI/r^3\), \(4EI/r^3\) and \(4.5EI/r^3\) of [pflueger1975Stability], p. 403 (case I, D, b, 1 of the appendix), which Table 3 of [schweizerhof1984Pressure] lists as reference for their finite element results. The follower value is also the hydrostatic critical pressure \(3D/r^3\) of the long cylinders, Eq. (66) of [nasa2020SP8007] and Eq. (97b), p. 150, of [pflueger1975Stability]. 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].
The Donnell values differ because the Donnell kinematics neglect \(v/r\) in the rotation \(\beta_y\), in the change of curvature \(\kappa_{yy}\) and in the area vector, \(k_v = 0\). [kardomateas1993Orthotropic] obtains the classical value \((n^2 - 1) D/r^3\), i.e. \(3D/r^3\) for \(n = 2\), Eq. (29a), p. 199, from the Donnell non-linear shell equations of [brush1975Buckling], but with a circumferential displacement that varies through the thickness with the slope \((v - w_{,\theta})/r\), in the notation of panels, which is the rotation \(\beta_y = -w_{,y} + v/r\) of the Sanders kinematics.
The centrally directed load has no kernel in panels: its column of the
table is checked only with the load stiffness central_load_stiffness
of tests/tests_shell/test_follower_pressure_validation.py, which is
integrated independently of the kernels.
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#
Eqs. (3), (7) and (10) are derived with SymPy in the notation of the kernels, which also verifies the symmetry identity Eq. (20) and the ring closed forms, and writes the kernels
calc_fext_follower(Eq. (7)) andfkCfollower_num(Eq. (10), assembled in full, not as an upper triangle) of the modulespanels.models.<model>_num, integrated withnx*nyGauss-Legendre points over each loaded patch, the same points as Eq. (8) inshell_fext().tests/tests_shell/test_follower_pressure.py: Eq. (8), Taylor tests of Eq. (7) against Eq. (10), finite-difference Jacobian of Eq. (14) against Eq. (15), quadratic convergence of Newton-Raphson and Riks, the symmetry conditions of Eq. (20) and Eq. (17).tests/multidomain/test_follower_pressure_multidomain.py: the follower loads in aMultiDomainassembly, and the symmetry of the load stiffness of a closed cylinder reduced by the null-space method.tests/tests_shell/test_follower_pressure_validation.pyand the notebooks of Follower (hydrostatic) pressure: the ring closed forms, finite cylinders, [nasa2020SP8007], [kardomateas1993Orthotropic], [han2004SandwichCylinder], [kardomateas2003Sandwich] and [schweizerhof1984Pressure].