Abstract

The geometrically exact nonlinear deflection of a beamshell is considered here as an extension of the formulation derived by Libai and Simmonds (1998, The Nonlinear Theory of Elastic Shells, Cambridge University Press, Cambridge, UK) to include deformation through the thickness of the beam, as might arise from transverse squeezing loads. In particular, this effect can lead to receding contact for a uniform beamshell resting on a smooth, flat, rigid surface; traditional shell theory cannot adequately such behavior. The formulation is developed from the weak form of the local equations for linear momentum balance, weighted by an appropriate tensor. Different choices for this tensor lead to both the traditional shell equations corresponding to linear and angular momentum balance, as well as the additional higher-order representation for the squeezing deformation. In addition, conjugate strains for the shell forces are derived from the deformation power, as presented by Libai and Simmonds. Finally, the predictions from this approach are compared against predictions from the finite element code abaqus for a uniform beam subject to transverse applied loads. The current geometrically exact shell model correctly predicts the transverse shell force through the thickness of the beamshell and is able to describe problems that admit receding contact.

1 Introduction

When objects are in contact, the extent of the contact area can vary with the applied compressive load between the bodies. In particular, as described by Dundurs and Stippes [1] and illustrated in Fig. 1, advancing contact describes situations in which the contact area under load ΓC is greater than that in the unloaded configuration, defined as Γ0, so that Γ0 ∈ ΓC. Stationary contact refers to examples where the contact area remains unchanged with the increasing load, so that ΓC = Γ0. In contrast, receding contact describes cases where the contact area decreases with the increasing load, that is, ΓC ∈ Γ0 [24]. Receding contact was first predicted by Filon [5] during his study of an elastic continuum squeezed by two equal and opposite forces, and it was found that if the continuum was sufficiently long compared to its height that at some point along its length, the transverse stress along the midplane changes from compressive to tensile. Later, Coker et al. [6] experimentally found that for an elastic strip on an elastic foundation subject to a point load, the ratio of the contact length to the height of the beam is constant, with a value of approximately 1.35. Dundurs and Stippes [1] found that for cases of receding contact with rectangular elastic blocks, the shape of the contact patch is discontinuous, forming immediately at load application and unchanging with an increase of load [7], while the magnitude of the stress field in the body is proportional to the load. Ahn and Barber [3] examined Dundur’s postulations including Coulomb friction. They found that under monotonic loading, the original conclusion holds. However, if the load is oscillatory, the shape and the extent of the contact patch changes.

Fig. 1
Classification of contact areas with the increasing load
Fig. 1
Classification of contact areas with the increasing load
Close modal

Receding contact also arises in the study of the nonlinear response of jointed structures, where a joint is defined as a connecting region for two or more structural components whose tangential motion relative to one another is restricted by friction [8]. These mechanical interfaces can play a significant role in the overall dissipation observed in mechanical systems, so that understanding the dynamics of joints and interfaces is an important component in the prediction of the overall structural response. In particular, localized regions of slip near the boundaries of the contact area, referred to as microslip, can play a significant role in the overall damping observed in the structure [8]. Heinstein and Segalman [9] investigated the transient geometry of the contact area in a joint and showed that the receding contact observed in the dynamic response area plays a significant role in the resulting energy dissipation.

Discrete Iwan models, formed by collections of spring-slider elements, have a natural continuum limit represented by a beam sliding on a rough foundation [1012] and can be incorporated into larger structural dynamics models [13]. Recent research has shown that hysteretic Iwan models are capable of describing the dissipation arising from microslip [1417]. More generally, shell theory provides a framework for the development of reduced-order models capable of representing the elastic and dissipative effects arising from mechanical joints. However, the traditional shell theory is unable to account for all effects in the response of jointed structures and in particular cannot capture receding contact. Standard shell theories provide no mechanism to introduce strains in the transverse direction through the thickness of the shell [18] although higher-order theories can incorporate such effects by making assumptions on the distribution of the strain [19]. By assuming a through thickness stress state, Essenberg [20] accurately predicts receding contact in a shell resting on a smooth rigid foundation.

This study introduces transverse strains within the geometrically exact nonlinear shell theory formulated by Libai and Simmonds [21] and generalizes the formulation of higher-order shell theories within this framework. Therefore, the predictive capability of the shell for problems in microslip can be improved by incorporating transverse strains and allowing for receding contact within the shell formulation. In particular, the work of Libai and Simmonds [21], and by extension the work herein, makes no specific assumptions on the distribution of strains throughout the three-dimensional continuum; instead, the shell reduction occurs by way of appropriately weighted through-the-thickness averaging. Thus, receding contact as well as the transition from compressive to tensile strain observed by Filon [5] is directly predicted by the theory without additional assumptions on the stress state.

2 Shell Formulation

In the following, bold underlined characters such as r and L represent vector quantities and bold characters with overhats denote unit vectors (e.g., e^), while characters in normal type (e.g., s, σ, ℓ) indicate scalars. Furthermore, bold capitalized quantities such as P or Z represent general tensor quantities. Finally, overdots (e.g., x˙) describe derivatives with respect to time, while primes (e.g., r′) represent derivatives with respect to the (to be specified) shell coordinate σ. Spatial derivatives are also represented as r,s denoting the derivative of r with respect to the spatial coordinate s, so that r′ ≡ r,σ.

2.1 General Equations.

This formulation begins with the geometrically exact shell theory defined and developed by Libai and Simmonds [21], and much of the initial description is repeated from that work for convenience. Consider a shell generated by translating a planar region normal to itself in the e^z direction, defined by Libai and Simmonds [21] as a beamshell. In particular, beamshells are assumed to undergo plane strain deformation in the e^z direction. The reference shape of the beamshell can be specified as follows:
(1)
with re^z=0. Here, r(σ, ζ) defines the generating planar region. We note that σ and ζ are the shell coordinates and are associated with the distance along the beamshell and normal through the thickness, respectively, although they need not be an actual arclength of material points. The curve defined by r(σ, 0) is known as the base curve, so that one can take
(2)
Finally, the deformed position of a particle that is initially located at r(σ, ζ) is denoted as r¯(σ,ζ) as illustrated in Fig. 2. A covariant basis for the beamshell can be defined as (r,σ,r,ζ,e^z), while a corresponding contravariant basis can be identified as (nσ,nζ,e^z), with
(3)
where μ=e^z(nσ×r,ζ).
Fig. 2
Shell coordinates and directions
Fig. 2
Shell coordinates and directions
Close modal
Consider a subvolume of the body defined by
(4)
For this control volume, linear momentum balance reduces to
(5)
where Sds is the traction acting on the boundary of A over the interval ds, and fdA represents the body force acting on the elements of the body within dA.
By virtue of the invariance of the beamshell in the e^z direction, the net force and motion in this direction vanish, so that the balance laws are only considered in the coordinates (σ, ζ). As a result,
(6)
With n representing the outward normal, on the edges of the volume σ = (a, b), so that nds=r,ζ×e^zdζ=nσμdζ. Likewise, on the faces, ζ = ±h and nds=e^z×r,σ=nζμdζ.2 Defining the first Piola–Kirchhoff stress tensor as P = r,σSσ + r,ζSζ, on the edges, Sds = Sσμ, while on the faces, Sds = Sζμ. Therefore,
(7)
With this, linear momentum balance from Eq. (5) can be written as follows:
(8)
and using the divergence theorem this reduces to
(9)
Thus, the local form of the linear momentum balance can be recovered as follows:
(10)

2.2 Generalized Shell Equations.

Now consider an arbitrary tensor T operating on the local form of linear momentum balance given in Eq. (10). When integrated over the region A, this yields
(11)
Note that
(12a)
(12b)
so that Eq. (11) can be written as follows:
(13)
Finally, using the divergence theorem, this reduces to
(14)
(15)
In terms of the beamshell coordinates (σ, ζ) over the volume A, Eq. (15) can be written as follows:
(16)

2.3 Linear and Angular Momentum Balance.

The integral forms of linear and angular momentum balance can be recovered by a proper choice of the tensor T and reproduce the shell equations found in Libai and Simmonds [21]. In particular, for T = Inir,i, then T=0 and d/dt(T) = 0. Moreover, the deformed position of material points can be expressed as r¯=y¯+z¯, where
(17)
Thus, y¯(σ,t) can be identified with the mass averaged position of material points located for constant σ, and
(18)
defined by Libai and Simmonds [21] as the dynamic consistency condition. As a result, Eq. (16) reduces to
(19)
with v¯y¯˙, and F and p defined as follows:
(20)
Here, F represents the internal force acting on the shell as illustrated in Fig. 2, while p is the external shell loading. This expression is identical to that identified in Eq. (IV.B.13) in Libai and Simmonds [21].
Angular momentum balance can be obtained by choosing T=Z¯, where Z¯ represents the skew-symmetric operator of the cross-product with z¯, so that aZ¯a×z¯. As a result,
(21)
Furthermore,
(22)
where aR¯a×r¯ and aY¯a×y¯. Note, however, that the local balance of angular momentum can be expressed as [22] follows:
(23)
so that
(24)
Making use of the dynamic consistency condition, angular momentum balance for the beamshell can be written as follows:
(25)
with M, l, and m¯ defined as follows:
(26)
and finally I is defined as follows:
(27)
As with linear momentum balance, this equation is identical to that obtained in Eq. (IV.B.17) by Libai and Simmonds [21], although Eq. (25) was obtained through a slightly different procedure, namely, with TZ¯ while those authors effectively use TR¯.

2.4 A Transverse Shell Equation.

Building on the previous work of Libai and Simmonds [21], higher-order shell theories can be developed by way of additional choices for T. In particular, the effect of compressive forces across the faces of the beamshell, which vanish from p in Eq. (20), can be included with the choice Tz¯. Therefore, returning to Eq. (16),
(28)
so that
(29)
with
(30)
and
(31)
with I defined in Eq. (27). The choice Tz¯ leads to a scalar equation in the direction of z¯, where η describes a weighted distribution of the stress over the cross section of the beam. Although this cannot be expressed as a conservation principle, it nonetheless describes the transverse shear distribution through the variable η, similar to the identification of the equivalent moment M in the equation for angular momentum balance given in Eq. (25).
Note that the term P:z¯ does not vanish and instead with z¯=z¯σr,σ+z¯ζr,ζ can be written as follows:
(32)
For linearized deformations, z¯/σ=0 and z¯/ζ=r,ζ so that P:z¯=σζζ, and this term arises from the squeezing stress transverse to the shell.3 Likewise, the linearization of Sσz¯ reduces to ζσσζ. Therefore, for small deformations, the internal states η and χ can be related to the continuum state of stress as follows:
(33)

2.5 Constitutive Equations.

The resulting differential equations for the beamshell can be written as follows:
(34)
These balance laws still require constitutive equations appropriate for the beamshell formulation, here based on the deformation power (identified below). For classical solutions, these equations can be multiplied by (V, Ω, U), respectively, combined, and integrated over the subegion (a, b) × (t1, t2). Finally, removing spatial derivatives on F, M, and η yields the weak form of the equations of motion
(35)
Taking V=v¯, Ω = ω, and U = u, this reduces to
(36)
known as the mechanical work identity, where
(37)
is the apparent external mechanical power,
(38)
is the reduced kinetic energy, and
(39)
is the deformation power. We note that these terms are identical to those found in Sec. IV.B in Libai and Simmonds [21] with the exception of the terms that depend on u, which arise from the additional higher-order equation given in Eq. (29). The contribution of the present work lies precisely in these additional terms together with the transverse beamshell equation given in Eq. (29).
Define a set of basis directions (T^,B^), which spin with angular speed ω. Thus, defining an angle β such that β˙=ω, the bending strain can be identified as follows:
(40)
so that ω=k˙. Further, y¯ can be expressed in terms of this basis as follows:
(41)
where e and g are defined as the extensional and bending strains, respectively. Then, the internal force F (compare with Eq. (20)) can be expressed in terms of these directions as follows:
(42)
Finally, the transverse strain ψ can be identified as ψ˙u, so that the deformation power reduces to
(43)
For an elastic system, the deformation power can be related to the strain energy density Φ as follows:
(44)
As a result, if Φ is a differentiable function of the strains, then
(45)
Therefore, since (e˙,g˙,k˙,ψ˙,ψ˙) can be chosen arbitrarily [23], Eq. (45) implies the constitutive equations:
(46)
The strains (e, g, k) describe the deformation of the curve describing the beamshell, while (ψ, ψ′) represent deformations transverse to the beamshell. Therefore, coupling in the strain energy density between (e, g, k) and either (ψ, ψ′) can couple the squeezing load defined by q to the deformation of the beamshell.

3 Isotropic Homogeneous Beamshell

In terms of the conjugate forces and strains identified above, the equilibrium response for Eq. (34) can be decomposed into the (T^,B^) directions to yield
(47)
where p=pTT^+pBB^. Note that these equations are geometrically exact and can be used to describe large deflections subject to appropriate constitutive equations. Linearizing these results in
(48)

3.1 Strain Energy Density.

An appropriate shell strain energy density Φ can be obtained by a descent from three to two dimensions, with specific kinematic assumptions on the three-dimensional deformation profile. The general strain energy is as follows:
(49)
where ϵ is the Lagrangian strain, defined as ϵ=12((r¯)T(r¯)I). In particular, the constitutive equation for isotropic linear elastic materials is as follows:
(50)
where E is the Young’s modulus and ν is the Poisson ratio of the material. Therefore, the strain energy reduces to
(51a)
(51b)
Determination of an appropriate strain energy density requires specific assumptions on the deformation field in the continuum and its relationship to the shell strain field. This procedure is described in Appendix  A, where the strain energy density is identified in Eq. (A11), and is given as follows:
(52)
Here, the parameters c1 and c2 describe correction coefficients that arise from assumed warping of material planes, similar to Timoshenko shear correction coefficients [24,25]. In the results given later in this study, these will be chosen based on comparison between predictions of the resulting shell model and equivalent predictions from finite element analysis (FEA).

With the kinematic assumptions described in Appendix  A, the shell strain energy density is therefore coupled between e and ψ, the extensional and transverse strains. Although a linearly elastic material model has been used to develop this strain energy density, this formulation retains geometric nonlinearities in the deformation arising from Eq. (34), and nonlinear material models could be easily used in place of Eq. (50).

3.2 Nondimensionalization.

Nondimensional and/or scaled variables are denoted by overhats of the corresponding dimensional quantities, for example, N^ represents the nondimensional normal shell force, while N is its dimensional counterpart. As described in Appendix  B, the strain energy density can be nondimensionalized to yield
(53)
where k^ is the nondimensional bending strain. The dimensional linearized equations are transformed to
(54)
where the nondimensional forces applied to the beamshell are defined in Appendix  B. In these equations, the nondimensional conjugate forces and strains are related as follows:
(55)
Note that the nondimensional parameters c1 and δ are correction coefficients related to the influence of the transverse strain energy while α is related to the transverse loading, defined in Appendix  A as follows:
(56)
As a result, the linearized nondimensional equations describing the deformation reduce to
(57)
If the beamshell is only subject to transverse loading, so that p^T=l^=0, the normal force N^ is constant, so that
(58)
Moreover, if N^=0, this implies that
(59)
and the differential shell equations can be written as follows:
(60a)
(60b)
with
(61)
Note that Eq. (60a) is identical to the standard transverse shell equation as described by an Euler–Bernoulli shell in the absence of shear or a Timoshenko shell with the inclusion of shear deformation. In contrast, Eq. (60b) describes the deformation arising from compressive external loading identified through q^. These equations can be coupled together through the external loadings p^ and q^. The dimensional external loads and shell forces are related to their nondimensional counterparts as given in Eqs. (B10) and (B14).

The resulting deformation of the beamshell is illustrated in the following three examples, where different loadings and boundary conditions are applied to a beamshell of length ℓ = 200 mm with an externally applied loading over a width 2 d, where d = 50 mm, vanishing outside this interval. The thickness of the beamshell is h = 1 mm, and the material properties are assumed to be E = 209 GPa and ν = 0.30. Each of the three example problems for the beamshell formulation was numerically verified using comparable continuum finite element problems developed in abaqus. In all three cases, the beamshell was meshed using uniform, square CPE4R elements, four-node linear plane strain quadrilateral continuum elements, with a constant side length of 0.025 mm. The mesh size was chosen to provide enough points near the edge of the load to resolve the phenomena studied in each example and provide comparable results to the beamshell model.

For algebraic convenience, in the nondimensional system, the spatial origin s = 0 is shifted so that it coincides with the edge of the loading. Therefore, the point of symmetry for the physical system is shifted to s=d^, and the deformation is symmetric about this location. For contact problems, the beamshell is assumed to rest on a smooth rigid surface, so that the normal contact pressure, defined as n^(s), is unknown and must be determined as part of the overall solution.

3.3 Uniform Squeezing.

Consider an isotropic homogeneous beamshell subjected to a uniform squeezing force of magnitude f^0, as shown in Fig. 3. Due to symmetry, s is only considered in the interval sd^. As a result, the external loading on the beamshell reduces to
(62)
Within traditional shell theory, the strains that result from this loading are identically zero, so that k^0—the beamshell undergoes no deformation and there are no internal forces developed in the shell. However, incorporating the squeezing strain ψ^ and the strain energy density given in Eq. (52), the squeeze strain couples with the extensional strain to create a Poisson effect in the direction transverse to the applied load.
Fig. 3
Uniform beamshell with squeezing force f^0. The spatial variable s is measured from the edge of the loading and the nondimensional width of the loading interval is 2d^, while the overall nondimensional beam length is 2ℓ^.
Fig. 3
Uniform beamshell with squeezing force f^0. The spatial variable s is measured from the edge of the loading and the nondimensional width of the loading interval is 2d^, while the overall nondimensional beam length is 2ℓ^.
Close modal
If the length of the beamshell is sufficiently long compared to its height, the beamshell can be considered in the limit ℓ → ∞. The response is symmetric about s=d^ and bounded, so that the boundary conditions reduce to
(63)
while ψ(s) is continuous and smooth at the edge of the loading, that is,
(64)
Equation (60a) vanishes identically, while Eq. (60b) reduces to a single linear differential equation with a piecewise constant forcing that can be solved subject to the above boundary conditions for ψ(s). Consequently, the transverse force χ^(s)κ2ψ(s) is expressed as follows:
(65)
As shown in Fig. 4(a) with κ/δ = 1, the beamshell allows the transverse beamshell stress χ^(s) to extend past the edge of the load application at s = 0 [5]. Further, the shape and the distribution of the squeezing stress as obtained from Eq. (65) are independent of f^0, the magnitude of the uniform external applied squeezing pressure, which is consistent with Dundurs and Stippes [1].
Fig. 4
Transverse beamshell force χ (the external load q0 is shown for reference as the dotted line): (a)χ^(s) (nondimensional; q^0/α=1.00, κ = 1, δ = 1) and (b) χ(σ). FEA, solid; beamshell (κ = 1.119, δ = 0.894), dashed.
Fig. 4
Transverse beamshell force χ (the external load q0 is shown for reference as the dotted line): (a)χ^(s) (nondimensional; q^0/α=1.00, κ = 1, δ = 1) and (b) χ(σ). FEA, solid; beamshell (κ = 1.119, δ = 0.894), dashed.
Close modal

Consider a uniform squeezing pressure of f0 = 100 MPa applied to both the top and bottoms surfaces along the mid-span of the beamshell, so that the nondimensional shell loading reduces to f^0=2.375×104. Note that the values of κ and δ depend on c1 and c2, respectively, the correction coefficients associated with the squeezing stress and identified in Appendix  A. In the limit κ/δ → ∞, the distribution of the squeezing stress χ(σ) approaches a step function, consistent with the traditional shell theory. In Fig. 4(b), the dimensional stress χ(σ) is shown for (κ, δ) = (1.119, 0.894), chosen based on the contact problem described in Sec. 3.4.1. Although these values slightly overpredict the extension of the transverse shell force χ beyond the loading interval, the predictions nevertheless compare favorably with the FEA results, evaluating the stress σyy integrated over the thickness as suggested by Eq. (33). Note that for these parameter values, the dimensional length scale is cs = 3.819 × 10−4, so that the nondimensional width of the loading interval is d^=130.93, while α = 0.764.

3.4 Contact.

As illustrated in the previous example, with the strain energy density given in Eq. (52), the transverse strain ψ can couple with the longitudinal strain e. In addition, the transverse strain can couple with the deformation of the beamshell through the loading and boundary conditions. It is this latter instance that gives rise to receding contact as the transverse stress extends beyond the edge of the load application. For this example, the beamshell presented in the previous example is set on a flat, frictionless, rigid surface assumed to be in the e^x direction, as shown in Fig. 5. The (nondimensional) external loading on the top of the beamshell is assumed to be f^(s)e^y, while the normal pressure from the foundation is unknown, but described as +n^(s)e^y, so that the shell loading is given as follows:
(66)
Therefore, the equations of motion for the shell reduce to
(67)
which can be combined to yield
(68)
However, the impenetrability of the foundation introduces a kinematic constraint on the continuum, so that 0=r¯,σ(σ,h/2)e^y, which is written in terms of the shell variables as follows:
(69)
Unilateral contact also requires that the normal load from the surface n(σ) is in the +e^y direction. This constraint can be linearized to yield
(70)
so that in terms of nondimensional variables,
(71)
Note that the parameter α was given in Eq. (56). Together with the linearized equations given in Eq. (60), the deformation equations can be written as follows:
(72)
These equations can be written in terms of the bending strain k^ alone as follows:
(73)
with f^0 for s > 0. In terms of k^, the transverse strain ψ can be written as follows:
(74)
Fig. 5
Contact with a smooth, rigid, flat surface
Fig. 5
Contact with a smooth, rigid, flat surface
Close modal
If the external force is discontinuous at s = 0, then d2k^/ds2 jumps as well, so that
(75)
where f^0 is the jump in the external force at the boundary of the loading, while
(76)
Also, the normal contact pressure can be discontinuous at the loading edge, so that
(77)
In the limit d^, the symmetry boundary conditions at s = −d are replaced by limsk^(s)=0.

These equations represent the analog of the traditional shell equations for contact on a rigid surface. Note that the general form is fourth order in the bending strain, similar to the expression obtained by Couchaux et al. [26], which was based on the theory of Baluch et al. [27]. These works enhance the stress state to account for transverse normal deformation. However, these theories fail to satisfy local equilibrium conditions and Hooke’s law, in contrast to the present approach that provides a consistent description of the shell equations.

The general solution to Eq. (73) can be written as follows:
(78)
The four values of λj are obtained by combinations of the two ± values in the aforementioned expression for λ. Finally, it is convenient to identify the parameter ϕ as follows:
(79)
Note that the parameters γ and ϕ are written in terms of κ and δ, which are directly related to the correction coefficients c1 and c2 as identified in Appendix  A. The parameters κ and δ are chosen to fit equivalent shell forces as obtained from the finite element analysis. For notational convenience, all four of these dependent parameters are retained in the results explained later.

3.4.1 Persistent Contact.

If contact is assumed to persist for the entire interval, then n^ can take both positive and negative values, with n^<0 implying an adhesive surface. In addition, as s → ∞, limsk^(s)=0. As a result, Eqs. (73)(76) can be solved for the bending strain as follows:
(80)
As a result, the normal pressure from the rigid surface becomes
(81)
In particular, consistent with the predictions of Filon [5], the normal contact traction decays exponentially outside the loading interval but alternates between a compressive and adhesive characteristic, so that the normal force vanishes at spatial locations s, where
(82)
The results predicted from the aforementioned beamshell model are compared against FEA results for the same physical parameters as used in the previous example. Note that the domain for the FEA calculations is σ ∈ (46, 52) mm. The components of the stress distribution from the FEA analysis are shown in Fig. 6 near the edge of the loading interval. Note that the shell model requires that the correction coefficients c1 and c2 be determined. In examining the FEA results, the normal pressure at the rigid surface is continuous at the edge of the loading, so that based on Eq. (77), (αδ)2 = 1/2, which implies that c2 = 3. The correction coefficient c1 can then be chosen to fit the shell predictions to those of the FEA. Minimizing the error between the moment and shear resultants, M and Q, suggests a value c1 = 0.885, so that
(83)
With these coefficients, the dimensional shell resultants are shown in Fig. 7 along with the equivalent stress resultants as obtained from FEA. The agreement between the moment and shear resultants shown in Figs. 7(a) and 7(b) is excellent for the identified shear correction coefficients, and the shell resultants χ and η are shown in Figs. 7(c) and 7(d), respectively. Note that the agreement for the transverse stress resultant χ is again good. In contrast, the agreement between the predicted and FEA quantities for η is not as good, but the comparison of qualitative results is favorable. Note that this quantity depends on the distribution of the shear stress across the beam and therefore is sensitive to the kinematic assumptions made in the development of the constitutive equations (compare with Appendix  A). Finally, the predicted contact pressure is shown in Fig. 8 together with the FEA results and again the agreement is excellent. In particular, both models predict an overshoot of the normal pressure on each side of the loading boundary compared to the external loading and the transition between compressive and adhesive contact is well approximated by Eq. (82).
Fig. 6
Persistent contact, stress field as predicted by FEA: (a) σxx, (b) σxy, and (c) σyy
Fig. 6
Persistent contact, stress field as predicted by FEA: (a) σxx, (b) σxy, and (c) σyy
Close modal
Fig. 7
Persistent contact, stress resultants (FEA, solid; beamshell, dashed): (a) Q, (b) M, (c) χ, and (d) η
Fig. 7
Persistent contact, stress resultants (FEA, solid; beamshell, dashed): (a) Q, (b) M, (c) χ, and (d) η
Close modal
Fig. 8
Persistent contact, contact pressure (FEA, solid; beamshell, dashed)
Fig. 8
Persistent contact, contact pressure (FEA, solid; beamshell, dashed)
Close modal

3.4.2 Receding Contact.

If the contact assumed to be unilateral then adhesion is prohibited so that contact only occurs for n^>0. If the contact interval extends to s=w^, then the spatial interval can be represented in terms of three domains, over which
(84)
The boundary conditions at s = 0 are still described by Eqs. (75) and (76), and the general solution for s<w^ remains as in Eq. (78). For s>w^, the strain field vanishes as s → ∞, which implies that
(85)
while ψ(s) decays exponentially, so that
(86)
Therefore, at s=w^,
(87)
corresponding to zero bending moment and shear force at s=w^. Using the continuity of ψ(s) and /ds(s), at the edge of the contact interval dψ/ds(w^)=(κ/δ)ψ(w^). Making use of Eq. (74), this boundary condition can be written as follows:
(88)
where w^ is the limit as sw^ with s<w^. Note that w^, the interval of contact, is also unknown and must be determined as part of the overall solution. The resulting boundary value problem is solved numerically to determine k^(s) and w^. However, note that these boundary conditions do not imply that n^(w^)=0. Rather, the normal pressure on the shell can jump at s=w^, so that
(89)
and then vanishes for s=w^+.

The components of the stress tensor obtained from FEA are shown in Fig. 9, which are similar to those obtained in Sec. 3.4.1 for persistent contact (compare with Fig. 6). However, note that the normal stress at the rigid surface (y = 0) is identically zero for σ > w, indicating a loss of contact. This is in contrast to the distribution seen in Fig. 6(c), where the normal stress only vanishes at isolated points as discussed in Sec. 3.4.1. The stress resultants and shell variables are compared in Fig. 10 for the same physical parameters as used in the previous example, while the shell coefficients κ and δ are again chosen based on the persistent contact problem. The resultant contact pressure is shown in Fig. 11 from both the FEA and the shell analyses. The predicted contact length from the FEA is wFEA = 50.40 mm, while for the shell analysis, wshell = 50.31 mm. Note that as expected, the normal load predicted by the shell jumps at σ = w, so that n(w) = 22.07 MPa before vanishing for σ > w+. Finally, in Fig. 12, the rotation angle β of the shell is shown. At the edge of the contact interval β(w) = 1.752 × 10−4 rad, so that the beamshell lifts off the surface with a small but finite angle of rotation. As predicted by Filon [5], the contact length is independent of the magnitude of the applied load, while the rotation angle at the edge of the contact interval is scales linearly with the applied load.

Fig. 9
Receding contact, stress field as predicted by FEA: (a) σxx, (b) σxy, and (c) σyy
Fig. 9
Receding contact, stress field as predicted by FEA: (a) σxx, (b) σxy, and (c) σyy
Close modal
Fig. 10
Receding contact, stress resultants; the edge of the contact interval at σ = w is indicated with the dotted line at w = 50.31 mm (FEA, solid; beamshell, dashed): (a) M, (b) Q, (c) χ, and (d) η
Fig. 10
Receding contact, stress resultants; the edge of the contact interval at σ = w is indicated with the dotted line at w = 50.31 mm (FEA, solid; beamshell, dashed): (a) M, (b) Q, (c) χ, and (d) η
Close modal
Fig. 11
Receding contact, contact pressure; the edge of the contact interval at σ = w is indicated with the dotted line at w = 50.31 mm (FEA, solid; beamshell, dashed)
Fig. 11
Receding contact, contact pressure; the edge of the contact interval at σ = w is indicated with the dotted line at w = 50.31 mm (FEA, solid; beamshell, dashed)
Close modal
Fig. 12
Receding contact, shell rotation angle β; the edge of the contact interval at σ = w is indicated with the dotted line at w = 50.31 mm
Fig. 12
Receding contact, shell rotation angle β; the edge of the contact interval at σ = w is indicated with the dotted line at w = 50.31 mm
Close modal

4 Conclusions

This work has derived a higher-order shell theory that can incorporate transverse squeezing stress coupled with the longitudinal and bending deformation of the shell. Following the geometrically exact theory presented by Libai and Simmonds [21], the current work also generalizes the development of higher-order shell theories through the choice of appropriate weighting operators T (compare with Eq. (16)). Finally, several examples are presented including uniform squeezing as well as receding contact. In each case, the inclusion of the coupling between the transverse squeezing load and the deformation of the shell cannot be captured by traditional shell theories, but allows for the higher-order shell theory developed here to agree with finite element simulations. Future efforts could extend the present formulation to consider multilayer laminates [28] through an appropriate description of the strain energy function, as well as the development of finite elements based on the proposed shell model.

Footnotes

2

“Edges” can be considered as internal “cross sections,” while “faces” are synonymous with the external surface of the body.

3

Note that the stress component σij represents the ij component of the stress tensor P in contrast to the arclength variable σ in the shell formulation.

Acknowledgment

This research was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the US Department of Energy’s National Nuclear Security Administration under contract DE–NA–0003525.

Appendix A: Strain Energy Density—Descent From Three Dimensions for a Uniform Beamshell

From Eq. (51), the continuum strain energy density is
(A1)
The deformed position of material points can be expressed in terms of the identified basis directions (T^,B^) as follows:
(A2)
so that the right Cauchy-Green strain tensor can be identified as follows:
(A3)
With the shell extensional and bending strains identified in Eqs. (40) and (41),
(A4a)
(A4b)
The Lagrangian strain tensor is defined as ε=12(GI), and the lowest order terms of its components can be expressed as follows:
(A5a)
(A5b)
(A5c)
Identification of a shell strain energy requires an assumption on the deformation of the three-dimensional continuum, here taken to be
(A6)
The spatial function f(ζ) describes the out-of-plane warping of material planes, while the material points are allowed to deform in the B^ direction for q(ζ) ≠ 0; taken together, these deformation profiles imply that Kirchhoff’s kinematic assumptions do not hold. Moreover, f(ζ) is assumed to be even, so that it is symmetric about the centerline of the beam, although the dynamic consistency condition given in Eq. (18) still holds. Although not formally specified, an example of the spatial shape meeting these requirements is as follows:
(A7)
Therefore, the strain field reduces to
(A8a)
(A8b)
(A8c)
As a result, the strain energy may be written as follows:
(A9)
with
(A10a)
(A10b)
In the expression for the strain energy, the constants c1 and c2 describe correction coefficients similar to the classic shear correction coefficients in the Timoshenko beam theory. These are typically chosen by either assuming a specific spatial distribution or comparison with predictions from the finite element analysis [24,25]. In addition, the final term outside the integral is constant and does not contribute to the continuum strain energy density, which is identified as follows:
(A11)

Appendix B: Nondimensionalization

In terms of the constitutive equations derived from the strain energy density, the equations of motion can be written as follows:
(B1)
The strain energy density and the spatial variable can be scaled as follows:
(B2)
In addition, both the bending strain k and ψ′ have units of length, so that
(B3)
where k^ is represents the nondimensional bending strain. As a result, the equations of motion reduce to
(B4)
With these scales, the strain energy density becomes
(B5)
The spatial length scale is therefore chosen as follows:
(B6)
and the nondimensional strain energy density can be identified as follows:
(B7)
where
(B8)
and
(B9)
Finally, the nondimensional loading on the shell is defined as follows:
(B10)
Note that the nondimensional loads, in particular pB and q that both arise from loading transverse to the shell, are scaled differently. Therefore, these differences are explicitly identified so that
(B11)
with
(B12)
The nondimensional shell resultants can be identified as follows:
(B13)
so that the dimensional and nondimensional shell forces are related as follows:
(B14)

References

1.
Dundurs
,
J.
, and
Stippes
,
M.
,
1970
, “
Role of Elastic Constants in Certain Contact Problems
,”
ASME J. Appl. Mech.
,
37
(
4
), pp.
965
970
. 10.1115/1.3408725
2.
Keer
,
L. M.
,
Dundurs
,
J.
, and
Tsai
,
K. C.
,
1972
, “
Problems Involving a Receding Contact Between a Layer and a Half Space
,”
ASME J. Appl. Mech.
,
39
(
4
), pp.
1115
1120
. 10.1115/1.3422839
3.
Ahn
,
Y. J.
, and
Barber
,
J. R.
,
2008
, “
Response of Frictional Receding Contact Problems
,”
Int. J. Mech. Sci.
,
50
(
10–11
), pp.
1519
1525
. 10.1016/j.ijmecsci.2008.08.003
4.
Barber
,
J. R.
,
2018
,
Contact Mechanics (Volume 250 of Solid Mechanics and Its Applications)
.
Springer Nature
,
New York
.
5.
Filon
,
L. N. G.
,
1903
, “
On an Approximate Solution for the Bending of a Beam of Rectangular Cross-Section Under Any System of Load
,”
Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.
,
206A
(
331–345
), pp.
63
155
.
6.
Coker
,
E. G.
,
Filon
,
L. N. G.
, and
Turner Jessop
,
H.
,
1957
,
A Treatise on Photo-Elasticity
, 2nd ed.,
Cambridge University Press
,
Cambridge, UK
.
7.
Tsai
,
K. C.
,
Dundurs
,
J.
, and
Keer
,
L. M.
,
1974
, “
Elastic Layer Pressed Against a Half Space
,”
ASME J. Appl. Mech.
,
41
(
3
), pp.
703
707
. 10.1115/1.3423375
8.
Gaul
,
L.
, and
Lenz
,
J.
,
1997
, “
Nonlinear Dynamics of Structures Assembled by Bolted Joints
,”
Acta Mech.
,
125
(
1–4
), pp.
169
181
. 10.1007/BF01177306
9.
Heinstein
,
M. W.
, and
Segalman
,
D. J.
,
2002
, “
Bending Effects in the Frictional Energy Dissipation in Lap Joints
,”
Sandia National Laboratories
,
Albuquerque, NM
,
Technical Report No. SAND2002–0083
.
10.
Menq
,
C.-H.
,
Bielak
,
J.
, and
Griffin
,
J. H.
,
1986
, “
The Influence of Microslip on Vibratory Response, Part I: A New Microslip Model
,”
J. Sound Vib.
,
107
(
2
), pp.
279
293
. 10.1016/0022-460X(86)90238-5
11.
Dane Quinn
,
D.
, and
Segalman
,
Daniel J.
,
2005
, “
Using Series-Series Iwan-Type Models for Understanding Joint Dynamics
,”
ASME J. Appl. Mech.
,
72
(
5
), pp.
778
784
. 10.1115/1.1875552
12.
Miller
,
J. D.
, and
Dane Quinn
,
D.
,
2009
, “
A Two-Sided Interface Model for Dissipation in Structural Systems With Frictional Joints
,”
J. Sound Vib.
,
321
(
1–2
), pp.
201
219
. 10.1016/j.jsv.2008.09.037
13.
Dane Quinn
,
D.
,
2012
, “
Modal Analysis of Jointed Structures
,”
J. Sound Vib.
,
331
(
1
), pp.
81
93
. 10.1016/j.jsv.2011.08.017
14.
Segalman
,
D. J.
,
2001
, “
An Initial Overview of Iwan Modeling for Mechanical Joints
,”
Sandia National Laboratories
,
Albuquerque, NM
,
Technical Report No. SAND2001–0811
.
15.
Song
,
Y.
,
Hartwigsen
,
C. J.
,
McFarland
,
D. M.
,
Vakakis
,
A. F.
, and
Bergman
,
L. A.
,
2004
, “
Simulation of Dynamics of Beam Structures With Bolted Joints Using Adjusted Iwan Beam Elements
,”
J. Sound Vib.
,
273
(
1–2
), pp.
249
276
. 10.1016/S0022-460X(03)00499-1
16.
Segalman
,
D. J.
,
2006
, “
Modelling Joint Friction in Structural Dynamics
,”
Struct. Control Health Monitor.
,
13
(
1
), pp.
430
453
. 10.1002/stc.119
17.
Deaner
,
B. J.
,
Allen
,
M. S.
,
Starr
,
M. J.
,
Segalman
,
D. J.
, and
Sumali
,
H.
,
2015
, “
Application of Viscous and Iwan Modal Damping Methods to Experimental Measurements From Bolted Structures
,”
ASME J. Vib. Acoust.
,
137
(
2
), p.
021012
. 10.1115/1.4029074
18.
Brink
,
A. R.
, and
Dane Quinn
,
D.
,
2016
, “
Shear Effects on Energy Dissipation From an Elastic Beam on a Rigid Foundation
,”
ASME J. Appl. Mech.
,
83
(
1
), p.
011004
. 10.1115/1.4031764
19.
Alijani
,
F.
, and
Amabili
,
M.
,
2014
, “
Non-Linear Static Bending and Forced Vibration of Rectangular Plates Retaining Non-Linearities in Rotations and Thickness Deformation
,”
Int. J. Non-Linear Mech.
,
67
, pp.
394
404
. 10.1016/j.ijnonlinmec.2014.10.003
20.
Essenberg
,
F.
,
1962
, “
Shear Deformation in Beams on Elastic Foundations
,”
ASME J. Appl. Mech.
,
29
(
2
), pp.
313
317
. 10.1115/1.3640547
21.
Libai
,
A.
, and
Simmonds
,
J. G.
,
1998
,
The Nonlinear Theory of Elastic Shells
, 2nd ed.,
Cambridge University Press
,
Cambridge, UK
.
22.
Antman
,
S. S.
,
1995
,
Nonlinear Problems of Elasticity
(
Number 107 in Applied Mathematical Sciences
),
Springer-Verlag
,
New York
.
23.
Coleman
,
B. D.
, and
Noll
,
W.
,
1963
, “
The Thermodynamics of Elastic Materials With Heat Conduction and Viscosity
,”
Arch. Rational Mech. Anal.
,
13
(
1
), pp.
167
178
. 10.1007/BF01262690
24.
Cowper
,
G. R.
,
1966
, “
The Shear Coefficient in Timoshenko’s Beam Theory
,”
ASME J. Appl. Mech.
,
33
(
2
), pp.
335
340
. 10.1115/1.3625046
25.
Hutchinson
,
J. R.
,
2001
, “
Shear Coefficients for Timoshenko Beam Theory
,”
ASME J. Appl. Mech.
,
68
(
1
), pp.
87
92
. 10.1115/1.1349417
26.
Couchaux
,
M.
,
Hjiaj
,
M.
, and
Ryan
,
I.
,
2015
, “
Enriched Beam Model for Slender Prismatic Solids in Contact With a Rigid Foundation
,”
Int. J. Mech. Sci.
,
93
, pp.
181
190
. 10.1016/j.ijmecsci.2014.12.012
27.
Baluch
,
M. H.
,
Azad
,
A. K.
, and
Khidir
,
M. A.
,
1984
, “
Technical Theory of Beams With Normal Strain
,”
J. Eng. Mech.
,
110
(
8
), pp.
1233
1237
. 10.1061/(ASCE)0733-9399(1984)110:8(1233)
28.
Reddy
,
J. N.
,
1984
, “
A Simple Higher-Order Theory for Laminated Composite Plates
,”
ASME J. Appl. Mech.
,
51
(
4
), pp.
745
752
. 10.1115/1.3167719