On Spencer’s displacement function approach for problems in second-order elasticity theory

Abstract

The paper describes the displacement function approach first proposed by AJM Spencer for the formulation and solution of problems in second-order elasticity theory. The displacement function approach for the second-order problem results in a single inhomogeneous partial differential equation of the form $E^4\psi (\mathbf{x})=f(\mathbf{x})$ , where $E^2$ is Stokes’ operator and $f(\mathbf{x})$ depends only on the first-order or the classical elasticity solution. The second-order isotropic stress $p(\mathbf{x})$ is governed by an inhomogeneous partial differential equation of the form ${\nabla }^{2}p(\mathbf{x})=g(\mathbf{x})$ , where ${\nabla }^2$ is Laplace’s operator and $g(\mathbf{x})$ depends only on the first-order or classical elasticity solution. The introduction of the displacement function enables the evaluation of the second-order displacement field purely through its derivatives and avoids the introduction of arbitrary rigid body terms normally associated with formulations where the strains need to be integrated. In principle, the displacement function approach can be systematically applied to examine higher-order effects, but such formulations entail considerable algebraic manipulations, which can be facilitated through the use of computer-aided symbolic mathematical operations. The paper describes the advances that have been made in the application of Spencer’s fundamental contribution and applies it to the solution of Kelvin’s concentrated force, Love’s doublet, and Boussinesq’s problems in second-order elasticity theory.

1. Introduction

Rubber-like elastic materials are characterized by their ability to sustain large strains without fracture and damage. The seminal works of RS Rivlin and co-workers presented a comprehensive and systematic approach to the description of hyperelastic materials, commencing with the formulation of constitutive relations by appeal to the theory of invariants, the solution of benchmark problems involving homogeneous strains, inflation and eversion of annular regions, torsion and bending of prismatic bodies, and performing experiments to validate the mathematical approach. References to these accomplishments are also contained in the collected works of RS Rivlin edited by Barenblatt and Joseph [1] and well documented [2–20] in terms of the contribution of these works to the overall advancement of non-linear continuum mechanics and, in particular, the development of constitutive equations for rubber-like elastic materials [21–35].

When dealing with hyperelastic materials, it is essential to consider the mathematical formulations and the relevant strain energy functions that can accommodate large strains in a consistent fashion. In addition, the realm of application of the hyperelasticity theory is also important. For example, the form of a strain energy function needed to model very large strains in biological tissues [36–43], and may have little relevance to applications dealing with the mechanics of industrial rubber mountings that experience moderately large strains and vice versa [44–51]. A case in point relates to the experiments reported by Rivlin [52] on the torsion of a pure gum rubber cylinder measuring 53.85 mm in diameter and 25.4 mm in length. With these dimensions, the problem is that of the torsion of a disk made of gum rubber and the torque required to induce measurable large strains can be substantial, to the extent that tensile tearing or even surface instabilities [53] can materialize.

The equations governing large strain behavior of rubber-like materials are highly non-linear and unless recourse is made to computational approaches [54–59], the possibility of obtaining analytical solutions to boundary value problems involving the complete non-linear hyperelasticity formulation with generalized forms of the strain energy function is perhaps unrealistic. For this reason, the early research applications in the area of finite elasticity also focused on developing solutions to problems through the method of successive approximations, where the first-order solution corresponded to the relevant classical elasticity solution. The mathematical developments in this area commence with the studies by Signorini [60–62], Mişicu [63], Stoppelli [64], Rivlin [65], Sheng [66], and Grioli [67]. The notion of higher-order theories of elasticity that can be obtained by the suitable expansion of the strains in terms of a small parameter is also discussed by Murnaghan [68,69] but without the development of a formal approach for the solution of boundary value problems. Also, it should be noted that the second-order effects of interest to this paper result from the manifestations of large strains in the hyperelastic medium as opposed to non-linearities that result from purely large deflections and large rotations of a medium that maintains small strain elastic behavior (e.g., the elastica). An appraisal of the work of Signorini is also provided by Capriz and Podio-Guidugli [70]. Reviews of the topic of second-order elasticity are given in Truesdell and Noll [7], Green and Adkins [8], Spencer [9] and Rivlin [65], and extensive developments and applications of the theory of second-order elasticity to boundary value problems are documented in the literature cited. For example, Adkins et al. [71] presented a complete development of the plane strain problem in finite elasticity in generalized curvilinear coordinates and with no restriction on the form of the strain energy function. Complex variable techniques [72,73] were used to solve certain benchmark plane strain problems. Adkins et al. [71] also discuss the application of the method of successive approximations, where the first-order approximation relates to the plane strain problem in classical elasticity and the second-order elasticity ensues. Green and Shield [74] continued the developments in the Adkins et al. [71] generalized tensor formulation of the second-order problem applicable to arbitrary forms of the strain energy function. A complex potential function approach was then used to arrive at the relevant integral equations that can be used to examine problems of axial tension and torsion in prismatic bodies that included elliptical cylinders. Adkins and Green [75] presented a complex variable approach to the study of two-dimensional problems in second-order elasticity theory. Blackburn and Green [76] also used two complex potential functions for the analysis of the second-order torsion problem and extended the studies to include bending of a cylinder. Blackburn [77] subsequently extended the work to include the torsion of a region bounded by a single closed curve of a transversely isotropic incompressible elastic material. The second-order bending of an isotropic incompressible cylinder by terminal couples was also examined by Blackburn [78] who reduced the problem to the solution of a single boundary value for two complex potential functions and the first-order problem corresponded to the classical boundary value problem. The axial extension that occurs during the twisting of cylinder was observed by Rivlin [52]. The problem was extended by Green [79] to include arbitrary forms of the strain energy function. The torsion problem was also examined by Sheng [66] using a stress function technique. Green and Spratt [80] investigated the second-order effects in the deformation of both compressible and incompressible elastic bodies and considered the torsion of a solid of revolution. A stress function is used to satisfy the incompressibility condition and equations of equilibrium for second-order elasticity in the axial and radial directions give rise to an inhomogeneous bi-harmonic equation for the stress function. Specific results for the torsion of a cone were presented in exact closed form and the second-order normal force generated during first-order torsion was evaluated explicitly. Application of complex variable methods to the study of two-dimensional problems in finite elasticity is discussed by Adkins et al. [81] and the use of the procedures to the formulation of problems in second-order elasticity for both compressible and incompressible elastic solids is also briefly discussed. The work of Carlson and Shield [82] continued the application of the successive approximation technique to include third-order effects applicable to a general class of plane problems of hyperelastic solids and specific solutions to boundary value problems are given for materials with a strain energy function of the Mooney–Rivlin type. Extensive use of complex variable methods is employed to develop the higher-order solutions to a range of problems of engineering interest. Chan and Carlson [83] examined the second-order incompressible elastic torsion and used, through several changes of the dependent variables, a direct procedure that reduces the second-order torsion problem to the solution of a two-dimensional classical linear elasticity problem without a pseudo body force term. The procedure was used to examine the second-order torsion of a bar with a square cross-section. Hill [84] re-examined the approach proposed by Chan and Carlson [83] and proved that the technique is completely general if the strain energy function for the incompressible elastic material is a symmetric function of the remaining principal invariants.

The application of a displacement function approach for the analysis of the second-order problem for incompressible elastic materials was first identified by Spencer in 1968 (see Appendix 1). A formal development of the application of the displacement function technique to axisymmetric problems in second-order elasticity for incompressible elastic materials was presented by Selvadurai and Spencer [85] and a collection of applications to both axisymmetric and two-dimensional problems in incompressible elasticity were presented by Selvadurai [86]. The strain energy function that was adopted for the analysis of problems in second-order elasticity was of the Mooney–Rivlin form and considering the expansion of the response functions in the general constitutive relationship for isotropic incompressible elastic materials, it can be shown that the Mooney–Rivlin form of the strain energy function is sufficient for describing the second-order effects in incompressible elastic materials. The basic methodology proposed by Spencer was adopted to formulate plane strain problems in second-order elasticity theory [87,88], the torsion of annular regions [89], and the axisymmetric loading of a rigid spherical inclusion in an elastic solid of infinite extent [90]. The studies were extended [91] to include axisymmetric problems dealing with spherical cavities and rigid spherical inclusions embedded in incompressible elastic media.

Other approaches have been used for the formulation and solution of problems in second-order elasticity theory. For example, Shield [92] used an energy method with second-order effects and applied it to examine the response of bars composed of compressible hyperelastic materials and subject to the action of torsion under an initial tension. Choi and Shield [93] used an inverse deformation approach to examine the category of second-order elasticity problems involving both the indentation of half-space regions and problems related to an infinite elastic solid containing a spherical cavity. Carroll and Rooney [94] utilize procedures whereby the solution of the second-order elasticity problem can be further facilitated by adopting suitable representations of the second-order incompressibility condition. Carroll and Rooney [95] also used the Strain Potential approach proposed by Love [96] to develop the second-order solution to Lord Kelvin’s problem [97] (see also [98–105]), which was originally developed and formulated in terms of Spencer’s displacement function by Selvadurai [86]. The second-order problem of the torsional indentation of an incompressible elastic half-space by a bonded flat circular punch was considered in the paper by Lindsay [106]. Goodman and Naghdi [107] adopt displacement functions that are similar to the Neuber–Papkovich representations for the solution of problems in linear elasticity [98–105,108] to examine the use of displacement functions in second-order elasticity theory for compressible elastic materials, and present solutions to certain plane strain problems in second-order elasticity theory involving compressible elastic materials. Lindsay [109] also considered the problem of the torsion of an incompressible slab; the problem is formulated by recourse to Weber–Orr transforms that are numerically inverted to produce relevant numerical results. Guo and Kaloni [110] have presented exact closed form solutions to the second-order elasticity problem for a compressible elastic half-space, which is subjected to a non-uniform shear load. The approach used involves the direct application of integral transform techniques to solve the analogous first-order or classical elasticity problem. In addition to the application of the theory of second-order elasticity to the study of three-dimensional problems, the method of successive approximations was applied by Kydoniefs and Spencer [111] to examine the finite inflation of a toroidal membrane. Similar successive approximation techniques were used by Kydoniefs [112] to examine elastic membrane problems. The scope of second-order elasticity in terms of potential applications to rubber-like materials in a technological setting makes the development of analytical approaches a worthwhile exercise.

In this study, we illustrate the application of Spencer’s displacement function approach for the formulation and solution of the second-order problem related to Kelvin’s fundamental problem, which deals with the application of a concentrated force at the interior of an incompressible elastic medium of infinite extent. The strain energy function for the elastic medium is assumed to be of the Mooney–Rivlin form. In the case of an incompressible elastic medium, the similarity between Kelvin’s problem for the concentrated force acting at the interior of an infinite space region and Boussinesq’s problem for a normal force acting at the surface of an elastic half-space is well known [99–101,103,113,114]. The displacement function approach gives explicit solutions to Kelvin’s problem and the methodology is also applied to develop a formal second-order solution to Love’s doublet problem and Boussinesq’s problem for a concentrated normal force on the surface of a half-space.

2. Governing equations

We follow the developments documented in Green and Adkins [8], Green and Spratt [80], Selvadurai and Spencer [85], and Selvadurai [86] and adopt the presentation for axially symmetric problems in terms of the displacement function approach. We denote particles in the reference configuration by $(R,\mathrm{\Theta },Z)$ and particles in the deformed configuration by $(r,\theta ,z)$ . The matrix of deformation gradients in the directions of $R,\mathrm{\Theta },Z$ is given by:

$$\mathbf{F}=\left(\begin{array}{ccc}\frac{\partial r}{\partial R}\hfill & 0\hfill & \frac{\partial r}{\partial Z}\hfill \\ \multicolumn{1}{c}{0}& \frac{r}{R}\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\partial z}{\partial R}}& 0\hfill & \frac{\partial z}{\partial Z}\hfill \end{array}\right).$$

(1)

We consider incompressible elastic materials for which,

$$det\mathbf{F}=\frac{r}{R}(\frac{\partial r}{\partial R}\frac{\partial z}{\partial Z}-\frac{\partial r}{\partial Z}\frac{\partial z}{\partial R})=1.$$

(2)

The matrix of physical components of the left Cauchy–Green strain tensor, referred to the $(r,\theta ,z)$ coordinates, is:

$$\mathbf{B}=\mathbf{F}{\mathbf{F}}^T=\left(\begin{array}{ccc}{\left(\frac{\partial r}{\partial R}\right)}^2+{\left(\frac{\partial r}{\partial Z}\right)}^2\hfill & 0\hfill & \frac{\partial r}{\partial R}\frac{\partial z}{\partial R}+\frac{\partial r}{\partial Z}\frac{\partial z}{\partial Z}\hfill \\ \multicolumn{1}{c}{0}& \frac{r^2}{R^2}\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\partial r}{\partial R}\frac{\partial z}{\partial R}+\frac{\partial r}{\partial Z}\frac{\partial z}{\partial Z}}& 0\hfill & {\left(\frac{\partial z}{\partial R}\right)}^2+{\left(\frac{\partial z}{\partial Z}\right)}^2\hfill \end{array}\right).$$

(3)

and, using equation (2), the inverse of equation (3) can be written as:

$${\mathbf{B}}^{-1}=\left(\begin{array}{ccc}\frac{r^2}{R^2}[{\left(\frac{\partial z}{\partial R}\right)}^2+{\left(\frac{\partial z}{\partial Z}\right)}^2]\hfill & 0\hfill & -\frac{r^2}{R^2}[\frac{\partial r}{\partial R}\frac{\partial z}{\partial R}+\frac{\partial r}{\partial Z}\frac{\partial z}{\partial Z}]\hfill \\ \multicolumn{1}{c}{0}& {[\frac{\partial r}{\partial R}\frac{\partial z}{\partial Z}-\frac{\partial r}{\partial Z}\frac{\partial z}{\partial R}]}^2\hfill & 0\hfill \\ \multicolumn{1}{c}{-\frac{r^2}{R^2}[\frac{\partial r}{\partial R}\frac{\partial z}{\partial R}+\frac{\partial r}{\partial Z}\frac{\partial z}{\partial Z}]}& 0\hfill & \frac{r^2}{R^2}[{\left(\frac{\partial r}{\partial R}\right)}^2+{\left(\frac{\partial r}{\partial Z}\right)}^2]\hfill \end{array}\right).$$

(4)

The basis for the method of successive approximations is that the displacement field can be expanded in a power series in terms of a small parameter, $\epsilon$ , which would naturally evolve in the analysis of the first-order problem; we assume that for axisymmetric states of deformation, the coordinates $(r,\theta ,z)$ can be represented in the form:

$$\begin{array}{cc}\multicolumn{1}{c}{r=}& R+\epsilon u_1(R,Z)+{\epsilon }^2{u}_2(R,Z)+\mathrm{O}\left({\epsilon }^3\right),\hfill \\ \multicolumn{1}{c}{\theta =}& \mathrm{\Theta },\hfill \\ \multicolumn{1}{c}{z=}& Z+\epsilon w_1(R,Z)+{\epsilon }^2{w}_2(R,Z)+\mathrm{O}\left({\epsilon }^3\right),\hfill \end{array}$$

(5)

where the suffixes $1$ and $2$ refer to the first-order and second-order components, respectively. Using equation (5) in equation (2), we obtain the first- and second-order incompressibility conditions as follows:

$$\begin{array}{c}\frac{\partial u_1}{\partial R}+\frac{u_1}{R}+\frac{\partial w_1}{\partial Z}=0,\hfill \\ \frac{\partial u_2}{\partial R}+\frac{u_2}{R}+\frac{\partial w_2}{\partial Z}=\frac{u_1^2}{R^2}+\frac{\partial u_1}{\partial Z}\frac{\partial w_1}{\partial R}-\frac{\partial u_1}{\partial R}\frac{\partial w_1}{\partial Z}=G(R,Z).\hfill \end{array}$$

(6)

For an isotropic incompressible elastic solid, the general form of the constitutive equation for the symmetric contravariant stress tensor $\mathbf{T}$ [80] can be reduced to the form:

$$\mathbf{T}=-p\mathbf{I}+{\mathrm{\Phi }}_1\mathbf{B}+{\mathrm{\Phi }}_{-1}{\mathbf{B}}^{-1},$$

(7)

where ${\mathrm{\Phi }}_i(I_1,I_2)(i=1,-1)$ are the scalar functions of the principal invariants $I_1$ and $I_2$ of the strain tensor $\mathbf{B}$ , $p$ is a scalar pressure and, for incompressibility, $I_3=1$ . We assume that the isotropic stress and the stress tensor can also be expanded in power series in terms of the parameter $\epsilon$ in the forms:

$$p=\sum _{n=1}^{\infty }{\epsilon }^n{p}_n(R,Z);\mathbf{T}=\sum _{n=1}^{\infty }{\epsilon }^n{\mathbf{T}}_n(R,Z).$$

(8)

The first-order constitutive equation can be reduced to the forms:

$$\begin{array}{cc}\multicolumn{1}{c}{T_{\mathrm{rr}}^{(1)}=}& -p_1+2\mu \frac{\partial u_1}{\partial R},\hfill \\ \multicolumn{1}{c}{T_{\theta \theta }^{(1)}=}& -p_1+2\mu \frac{u_1}{R},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{zz}}^{(1)}=}& -p_1+2\mu \frac{\partial w_1}{\partial Z},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rz}}^{(1)}=}& \mu (\frac{\partial u_1}{\partial Z}+\frac{\partial w_1}{\partial R}),\hfill \end{array}$$

(9)

where $\mu (=2(C_1+C_2))$ is the linear elastic shear modulus, and the constants $C_1$ and $C_2$ can be identified with the material constants characterizing the Mooney–Rivlin form of the strain energy function:

$$W(I_1,I_2)=C_1(I_1-3)+C_2(I_2-3).$$

(10)

The constitutive equations for the second-order stress components can be written as:

$$\begin{array}{cc}\multicolumn{1}{c}{T_{\mathrm{rr}}^{(2)}=}& -p_2+2\mu \frac{\partial u_2}{\partial R}+t_{\mathrm{rr}},\hfill \\ \multicolumn{1}{c}{T_{\theta \theta }^{(2)}=}& -p_2+2\mu \frac{u_2}{R}+t_{\theta \theta },\hfill \\ \multicolumn{1}{c}{T_{\mathrm{zz}}^{(2)}=}& -p_2+2\mu \frac{\partial w_2}{\partial Z}+t_{\mathrm{zz}},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rz}}^{(2)}=}& \mu (\frac{\partial u_2}{\partial Z}+\frac{\partial w_2}{\partial R})+t_{\mathrm{rz}},\hfill \end{array}$$

(11)

and the components $t_{\mathrm{rr}},t_{\theta \theta },t_{\mathrm{zz}}$ , and $t_{\mathrm{rz}}$ are given by:

$$\begin{array}{cc}\multicolumn{1}{c}{t_{\mathrm{rr}}=}& \mu [{\left(\frac{\partial u_1}{\partial R}\right)}^2+{\left(\frac{\partial u_1}{\partial Z}\right)}^2]-2C_2[4{\left(\frac{\partial u_1}{\partial R}\right)}^2+{(\frac{\partial w_1}{\partial R}+\frac{\partial u_1}{\partial Z})}^2],\hfill \\ \multicolumn{1}{c}{t_{\theta \theta }=}& (\mu -8C_2){\left(\frac{u_1}{R}\right)}^2,\hfill \\ \multicolumn{1}{c}{t_{\mathrm{zz}}=}& \mu [{\left(\frac{\partial w_1}{\partial Z}\right)}^2+{\left(\frac{\partial w_1}{\partial R}\right)}^2]-2C_2[4{\left(\frac{\partial w_1}{\partial Z}\right)}^2+{(\frac{\partial w_1}{\partial R}+\frac{\partial u_1}{\partial Z})}^2],\hfill \\ \multicolumn{1}{c}{t_{\mathrm{rz}}=}& \mu [\frac{\partial u_1}{\partial R}\frac{\partial w_1}{\partial R}+\frac{\partial u_1}{\partial Z}\frac{\partial w_1}{\partial Z}]-4C_2[\frac{\partial u_1}{\partial R}+\frac{\partial w_1}{\partial Z}][\frac{\partial w_1}{\partial R}+\frac{\partial u_1}{\partial Z}].\hfill \end{array}$$

(12)

For axial symmetry, the equations of equilibrium for the symmetric contravariant stress $\mathbf{T}$ referred to the deformed coordinate system can be written as:

$$\begin{array}{c}\frac{\partial T_{\mathrm{rr}}}{\partial r}+\frac{\partial T_{\mathrm{rz}}}{\partial z}+\frac{T_{\mathrm{rr}}-T_{\theta \theta }}{r}=0,\hfill \\ \frac{\partial T_{\mathrm{rz}}}{\partial r}+\frac{\partial T_{\mathrm{zz}}}{\partial z}+\frac{T_{\mathrm{rz}}}{r}=0.\hfill \end{array}$$

(13)

Using the expansions in powers of $\epsilon$ and noting the derivatives:

$$\begin{array}{c}\frac{\partial }{\partial r}=(1-\epsilon \frac{\partial u_1}{\partial R})\frac{\partial }{\partial R}-\epsilon \frac{\partial w_1}{\partial R}\frac{\partial }{\partial Z},\hfill \\ \frac{\partial }{\partial z}=(1-\epsilon \frac{\partial w_1}{\partial Z})\frac{\partial }{\partial Z}-\epsilon \frac{\partial u_1}{\partial Z}\frac{\partial }{\partial R},\hfill \end{array}$$

(14)

the first-order equations of equilibrium take the forms:

$$\begin{array}{c}\frac{\partial T_{\mathrm{rr}}^{(1)}}{\partial R}+\frac{\partial T_{\mathrm{rz}}^{(1)}}{\partial Z}+\frac{T_{\mathrm{rr}}^{(1)}-T_{\theta \theta }^{(1)}}{R}=0,\hfill \\ \frac{\partial T_{\mathrm{rz}}^{(1)}}{\partial R}+\frac{\partial T_{\mathrm{zz}}^{(1)}}{\partial Z}+\frac{2T_{\mathrm{rz}}^{(1)}}{R}=0,\hfill \end{array}$$

(15)

and the second-order equations of equilibrium take the forms:

$$\begin{array}{c}\frac{\partial T_{\mathrm{rr}}^{(2)}}{\partial R}+\frac{\partial T_{\mathrm{rz}}^{(2)}}{\partial Z}+\frac{T_{\mathrm{rr}}^{(2)}-T_{\theta \theta }^{(2)}}{R}=H_R(R,Z),\hfill \\ \frac{\partial T_{\mathrm{rz}}^{(2)}}{\partial R}+\frac{\partial T_{\mathrm{zz}}^{(2)}}{\partial Z}+\frac{T_{\mathrm{rz}}^{(2)}}{R}=H_Z(R,Z),\hfill \end{array}$$

(16)

where:

$$\begin{array}{c}{H}_R(R,Z)=\frac{\partial u_1}{\partial R}\frac{\partial T_{\mathrm{rr}}^{(1)}}{\partial R}+\frac{\partial w_1}{\partial R}\frac{\partial T_{\mathrm{rr}}^{(1)}}{\partial Z}+\frac{\partial w_1}{\partial Z}\frac{\partial T_{\mathrm{zz}}^{(1)}}{\partial Z}+\frac{\partial u_1}{\partial Z}\frac{\partial T_{\mathrm{rz}}^{(1)}}{\partial R}+\frac{u_1}{R^2}(T_{\mathrm{rr}}^{(1)}-T_{\theta \theta }^{(1)}),\hfill \\ H_Z(R,Z)=\frac{\partial u_1}{\partial R}\frac{\partial T_{\mathrm{rz}}^{(1)}}{\partial R}+\frac{\partial w_1}{\partial R}\frac{\partial T_{\mathrm{rz}}^{(1)}}{\partial Z}+\frac{\partial w_1}{\partial Z}\frac{\partial T_{\mathrm{zz}}^{(1)}}{\partial Z}+\frac{\partial u_1}{\partial Z}\frac{\partial T_{\mathrm{zz}}^{(1)}}{\partial R}+\frac{u_1}{R^2}{T}_{\mathrm{rz}}^{(1)}.\hfill \end{array}$$

(17)

The second-order equations of equilibrium can be expressed in the forms:

$$\begin{array}{c}-\frac{\partial p_2}{\partial R}+\mu (2\frac{{\partial }^2{u}_2}{\partial R^2}+\frac{{\partial }^2{w}_2}{\partial R\partial Z}+\frac{{\partial }^2{u}_2}{\partial Z^2}+\frac{2}{R}\frac{\partial u_2}{\partial R}-\frac{2u_2}{R^2})+\frac{\partial t_{\mathrm{rr}}}{\partial R}+\frac{\partial t_{\mathrm{rz}}}{\partial Z}+\frac{t_{\mathrm{rr}}-t_{\theta \theta }}{R}=H_R(R,Z),\hfill \\ -\frac{\partial p_2}{\partial Z}+\mu (2\frac{{\partial }^2{w}_2}{\partial Z^2}+\frac{{\partial }^2{u}_2}{\partial R\partial Z}+\frac{{\partial }^2{w}_2}{\partial R^2}+\frac{1}{R}\frac{\partial u_2}{\partial Z}+\frac{1}{R}\frac{\partial w_2}{\partial R})+\frac{\partial t_{\mathrm{rz}}}{\partial R}+\frac{\partial t_{\mathrm{zz}}}{\partial Z}+\frac{t_{\mathrm{rz}}}{R}=H_Z(R,Z).\hfill \end{array}$$

(18)

3. Spencer’s displacement functions

In the formulation of the first-order problem in incompressible elasticity, Spencer introduced the approach involving Stokes’ stream function, which satisfies the incompressibility condition for any choice of the function. This approach for the first-order problem bears a similarity with the problem of slow viscous flow for Newtonian fluids in terms of Stokes’ stream function (Lamb [115], Happel and Brenner [116], Langlois and Deville [117], Constantinescu [118]). This aspect has been observed by many researchers including Lord Rayleigh [119], Goodier [120], Hill [121], Prager [122], Adkins [123], Collins [124], and Richards [125]. The extension of the approach to the formulation of problems in second-order elasticity for incompressible elastic materials is not straightforward and techniques have to be developed to formulate the expressions for the second-order displacement field in terms of the second-order displacement function in such way that the second-order incompressibility condition can be satisfied.

Considering the first-order problem, we introduce a displacement function ${\mathrm{\Psi }}_1(R,Z)$ such that the first-order displacement components $u_1(R,Z)$ and $w_1(R,Z)$ are given by:

$$u_1=\frac{1}{R}\frac{\partial {\mathrm{\Psi }}_1}{\partial Z};w_1=-\frac{1}{R}\frac{\partial {\mathrm{\Psi }}_1}{\partial R}.$$

(19)

Considering these representations and first-order stresses (9) and the first-order equations of equilibrium, we arrive at the following partial differential equations governing the first-order displacement function ${\mathrm{\Psi }}_1(R,Z)$ and the first-order isotropic stress $p_1(R,Z)$ , i.e.,

$$E^4{\mathrm{\Psi }}_1(R,Z)=0;{\nabla }^{2}p(R,Z)=0,$$

(20)

where $E^2$ and ${\nabla }^2$ are, respectively, Stokes’ operator and Laplace’s operator defined by:

$$\begin{array}{cc}\multicolumn{1}{c}{E^2=}& \frac{{\partial }^2}{\partial R^2}-\frac{1}{R}\frac{\partial }{\partial R}+\frac{{\partial }^2}{\partial Z^2};E^4=E^2{E}^2,\hfill \\ \multicolumn{1}{c}{{\nabla }^2=}& \frac{{\partial }^2}{\partial R^2}+\frac{1}{R}\frac{\partial }{\partial R}+\frac{{\partial }^2}{\partial Z^2}.\hfill \end{array}$$

(21)

The extension of the displacement approach to the formulation of the second-order problem requires a representation that can exactly satisfy the second-order incompressibility condition given by equation (6). Methods that involve other techniques in terms of displacement and stress functions were proposed by Green and Spratt [80] and Chan and Carlson [83], but Spencer’s formulation retains the basic displacement function approach adopted for the first-order problem; the representation of the second-order displacement components in terms of the second-order displacement function have the forms:

$$u_2=\frac{1}{R}\frac{\partial {\mathrm{\Psi }}_2}{\partial Z}+\frac{u_1}{2}(\frac{\partial u_1}{\partial R}-\frac{\partial w_1}{\partial Z});w_2=-\frac{1}{R}\frac{\partial {\mathrm{\Psi }}_1}{\partial R}+u_1\frac{\partial w_1}{\partial R}.$$

(22)

As was indicated by Selvadurai and Spencer [85], there are other representations of the type (22) that can satisfy the second-order incompressibility condition (6) and they are admissible without exception. The second-order equations of equilibrium (18) can now be reduced to forms in terms of the second-order displacement function ${\mathrm{\Psi }}_2(R,Z)$ and the second-order isotopic stress $p_2(R,Z)$ in the forms:

$$\begin{array}{c}-\frac{\partial p_2}{\partial R}+\mu \left(\frac{1}{R}\frac{\partial }{\partial Z}\left(E^2{\mathrm{\Psi }}_2\right)\right)=F_R(R,Z),\hfill \\ -\frac{\partial p_2}{\partial Z}-\mu \left(\frac{1}{R}\frac{\partial }{\partial R}\left(E^2{\mathrm{\Psi }}_2\right)\right)=F_Z(R,Z),\hfill \end{array}$$

(23)

where:

$$\begin{array}{cc}\multicolumn{1}{c}{F_R(R,Z)=}& H_R(R,Z)-(\frac{\partial {t\prime }_{\mathrm{rr}}}{\partial R}+\frac{\partial {t\prime }_{\mathrm{rz}}}{\partial Z}+\frac{{t\prime }_{\mathrm{rr}}-{t\prime }_{\theta \theta }}{R}),\hfill \\ \multicolumn{1}{c}{F_Z(R,Z)=}& H_Z(R,Z)-(\frac{\partial {t\prime }_{\mathrm{rz}}}{\partial R}+\frac{\partial {t\prime }_{\mathrm{zz}}}{\partial Z}+\frac{{t\prime }_{\mathrm{rz}}}{R}),\hfill \end{array}$$

(24)

and

$$\begin{array}{c}{t\prime }_{\mathrm{rr}}=\mu \{u_1(\frac{u_1}{R^2}-2\frac{{\partial }^2{w}_1}{\partial R\partial Z})+3{\left(\frac{\partial u_1}{\partial R}\right)}^2+{\left(\frac{\partial u_1}{\partial Z}\right)}^2\}-2C_2\{4{\left(\frac{\partial u_1}{\partial R}\right)}^2+{(\frac{\partial u_1}{\partial Z}+\frac{\partial w_1}{\partial R})}^2\},\hfill \\ {t\prime }_{\theta \theta }=-2\mu \left(\frac{u_1}{R}\right)\left(\frac{\partial w_1}{\partial Z}\right)-8C_2{\left(\frac{u_1}{R}\right)}^2,\hfill \\ {t\prime }_{\mathrm{zz}}=\mu \{2u_1\frac{{\partial }^2{w}_1}{\partial R\partial Z}+2\left(\frac{\partial u_1}{\partial Z}\right)\left(\frac{\partial w_1}{\partial R}\right)+{\left(\frac{\partial w_1}{\partial R}\right)}^2+{\left(\frac{\partial w_1}{\partial Z}\right)}^2\}-2C_2\{4{\left(\frac{\partial w_1}{\partial Z}\right)}^2+{(\frac{\partial u_1}{\partial Z}+\frac{\partial w_1}{\partial R})}^2\},\hfill \\ {t\prime }_{\mathrm{rz}}=\mu \{u_1(\frac{{\partial }^2{u}_1}{\partial R\partial Z}+\frac{{\partial }^2{w}_1}{\partial R^2})+2\left(\frac{\partial u_1}{\partial R}\right)\left(\frac{\partial w_1}{\partial R}\right)\}+4C_2\left\{\frac{u_1}{R}(\frac{\partial u_1}{\partial Z}+\frac{\partial w_1}{\partial R})\right\}.\hfill \end{array}$$

(25)

By successively eliminating $p_2(R,Z)$ and ${\mathrm{\Psi }}_2(R,Z)$ from equation (23), we obtain the following partial differential equations:

$$\begin{array}{cc}\multicolumn{1}{c}{E^4{\mathrm{\Psi }}_2(R,Z)=}& \frac{R}{\mu }(\frac{\partial F_R}{\partial Z}-\frac{\partial F_Z}{\partial R}),\hfill \\ \multicolumn{1}{c}{{\nabla }^2{p}_2(R,Z)=}& -(\frac{\partial F_R}{\partial R}+\frac{F_R}{R}+\frac{\partial F_Z}{\partial Z}),\hfill \end{array}$$

(26)

which govern the axisymmetric second-order problem for an incompressible elastic material. For the solution of specific boundary value problems, both traction and/or displacement boundary conditions need to be specified on the deformed boundaries.

3.1. Boundary conditions

The development of solutions to the first-order problem is straightforward and either traction or displacement boundary conditions could be specified on surfaces referred to the undeformed configuration. In the case of the second-order problem, the displacement and traction boundary conditions have to be specified in relation to surfaces that are prescribed on the deformed body. The components of the traction vector on a surface in the deformed body, resolved in the $r-$ and $z-\mathrm{directions}$ are:

$$F_r=n_r{T}_{\mathrm{rr}}+n_z{T}_{\mathrm{rz}};F_z=n_r{T}_{\mathrm{rz}}+n_z{T}_{\mathrm{zz}},$$

(27)

where:

$$n_r^2+n_z^2=1;\frac{n_r}{n_z}=\frac{\frac{\partial F}{\partial r}}{\frac{\partial F}{\partial z}},$$

(28)

and the traction boundary conditions are normally specified in relation to a surface in the deformed body defined by:

$$F(r,z)=\tilde{F}(R,Z)=0.$$

(29)

From equation (28), the normal to the deformed surface has components, $n_r^{(1)}+\epsilon n_r^{(2)}$ and $n_z^{(1)}+\epsilon n_z^{(2)}$ to order $\epsilon$ , where:

$$\begin{array}{c}{(n_r^{(1)}+\epsilon n_r^{(2)})}^2+{(n_z^{(1)}+\epsilon n_z^{(2)})}^2=1,\hfill \\ \frac{n_r^{(1)}+\epsilon n_r^{(2)}}{n_z^{(1)}+\epsilon n_z^{(2)}}=\frac{\frac{\partial \tilde{F}}{\partial R}-\epsilon (\frac{\partial u_1}{\partial R}\frac{\partial \tilde{F}}{\partial R}+\frac{\partial w_1}{\partial R}\frac{\partial \tilde{F}}{\partial Z})}{\frac{\partial \tilde{F}}{\partial Z}-\epsilon (\frac{\partial u_1}{\partial Z}\frac{\partial \tilde{F}}{\partial R}+\frac{\partial w_1}{\partial Z}\frac{\partial \tilde{F}}{\partial Z})}.\hfill \end{array}$$

(30)

Using these relationships, the expressions (27) can be expressed as:

$$\begin{array}{c}{F}_r=\epsilon \{n_r^{(1)}{T}_{\mathrm{rr}}^{(1)}+n_z^{(1)}{T}_{\mathrm{rz}}^{(1)}\}+{\epsilon }^2\{n_r^{(1)}{T}_{\mathrm{rr}}^{(2)}+n_z^{(1)}{T}_{\mathrm{rz}}^{(2)}+n_r^{(2)}{T}_{\mathrm{rr}}^{(1)}+n_z^{(2)}{T}_{\mathrm{rz}}^{(1)}\},\hfill \\ F_z=\epsilon \{n_r^{(1)}{T}_{\mathrm{rz}}^{(1)}+n_z^{(1)}{T}_{\mathrm{zz}}^{(1)}\}+{\epsilon }^2\{n_r^{(1)}{T}_{\mathrm{rr}}^{(2)}+n_z^{(1)}{T}_{\mathrm{zz}}^{(2)}+n_r^{(2)}{T}_{\mathrm{rz}}^{(1)}+n_z^{(2)}{T}_{\mathrm{zz}}^{(1)}\}.\hfill \end{array}$$

(31)

4. Kelvin’s problem

The state of stress in an elastic solid of infinite extent subjected to a concentrated force of magnitude $P_K$ applied at the origin of coordinates of a system of cylindrical polar coordinates was first solved by Lord Kelvin [97] (Figure 1). When the line of action of the concentrated force coincides with the $Z-\mathrm{axis}$ of the cylindrical polar coordinate system $(R,\mathrm{\Theta },Z)$ , the problem is axisymmetric and the solution to the first-order or classical elasticity problem can be found in a variety of ways, which are summarized in texts on classical elasticity and its applications [98–105]. For example, if Love’s strain potential approach is applied to the solution of Kelvin’s problem with axial symmetry, the analysis requires the use of the strain potential $\mathrm{\Omega }(R,Z)$ , which satisfies the bi-harmonic equation:

$${\nabla }^4\mathrm{\Omega }(R,Z)=0,$$

(32)

where the displacement and stress components are given by:

$$\begin{array}{c}2\mu u_1=-\frac{{\partial }^2\mathrm{\Omega }}{\partial R\partial Z};2\mu w_1=2(1-\nu ){\nabla }^2\mathrm{\Omega }-\frac{{\partial }^2\mathrm{\Omega }}{\partial Z^2},\hfill \\ T_{\mathrm{rr}}^{(1)}=\frac{\partial }{\partial Z}(\nu {\nabla }^2\mathrm{\Omega }-\frac{{\partial }^2\mathrm{\Omega }}{\partial R^2});T_{\theta \theta }^{(1)}=\frac{\partial }{\partial Z}(\nu {\nabla }^2\mathrm{\Omega }-\frac{1}{R}\frac{\partial \mathrm{\Omega }}{\partial R});T_{\mathrm{zz}}^{(1)}=\frac{\partial }{\partial Z}((2-\nu ){\nabla }^2\mathrm{\Omega }-\frac{{\partial }^2\mathrm{\Omega }}{\partial Z^2}),\hfill \\ T_{\mathrm{rz}}^{(1)}=\frac{\partial }{\partial R}((1-\nu ){\nabla }^2\mathrm{\Omega }-\frac{{\partial }^2\mathrm{\Omega }}{\partial Z^2}),\hfill \end{array}$$

(33)

and $\nu$ is Poisson’s ratio.

Figure 1.

Kelvin’s problem.

A function that furnishes a solution to Kelvin’s problem is Love’s potential:

$$\mathrm{\Omega }(R,Z)=B{(R^2+Z^2)}^{1/2},$$

(34)

where $B$ is an arbitrary constant. Restricting the analysis to the case of an incompressible elastic material, we can show that the first-order displacement and stress components take the forms:

$$\begin{array}{cc}\multicolumn{1}{c}{2\mu u_1=}& B\frac{\mathrm{RZ}}{{(R^2+Z^2)}^{3/2}};2\mu w_1=B\frac{R^2+2Z^2}{{(R^2+Z^2)}^{3/2}},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rr}}^{(1)}=}& B(-3\frac{R^{2}Z}{{(R^2+Z^2)}^{5/2}});T_{\theta \theta }^{(1)}=0,\hfill \\ \multicolumn{1}{c}{T_{\mathrm{zz}}^{(1)}=}& B(-3\frac{Z^3}{{(R^2+Z^2)}^{5/2}});T_{\mathrm{rz}}^{(1)}=B(-3\frac{R{Z}^2}{{(R^2+Z^2)}^{5/2}}).\hfill \end{array}$$

(35)

The constant $B$ can be determined by calculating the resultant force acting on a cylindrical or spherical surface enclosing the origin. This gives:

$$B=\frac{P_K}{4\pi }.$$

(36)

It is noted that the first-order displacement and stress fields reduce to zero as $(R,Z)\to \infty$ , and they are singular as $(R,Z)\to 0$ . This is a constraint of the first-order solution for Kelvin’s problem. The only means of maintaining the first-order solution bounded as $(R,Z)\to 0$ is to physically exclude the origin of coordinates through the provision of either a cavity or an inclusion of finite dimension. This approach was used by Chadwick and Trowbridge [126], Selvadurai [90,127,128] and Selvadurai and Dasgupta [129] by including either a rigid spherical or spheroidal inclusion at the origin and satisfying the appropriate displacement and/or traction boundary conditions on the surface of the inclusion. A second-order elasticity problem of the centrally loaded spherical rigid inclusion analogue of Kelvin’s problem was examined by Selvadurai [90]. In order to develop the second-order solution to the localized Kelvin force problem, we utilize the first-order solution given by equation (35). Avoiding details of lengthy algebraic manipulations, it can be shown that the partial differential equation governing the second-order displacement function takes the form:

$$E^4{\mathrm{\Psi }}_2=-\frac{3P_K^2}{4{\pi }^2{\mu }^2}(\frac{2C_2}{\mu }-1)\frac{Z{R}^2}{{(R^2+Z^2)}^4}.$$

(37)

The particular integral of equation (37) can be obtained by converting the equation into a spherical polar coordinate form, where:

$$\begin{array}{c}R=S\sin \mathrm{\Theta };Z=S\cos \mathrm{\Theta };S={(R^2+Z^2)}^{1/2}\hfill \\ E^2=(\frac{{\partial }^2}{\partial S^2}+\frac{1}{S^2}\frac{{\partial }^2}{\partial {\mathrm{\Theta }}^2}-\frac{\cot \mathrm{\Theta }}{S^2}\frac{\partial }{\partial \mathrm{\Theta }}).\hfill \end{array}$$

(38)

It can be shown that the relevant expressions for the displacement and stress components derived from the particular solution of equation (37) are:

$$\begin{array}{cc}\multicolumn{1}{c}{u_{2p}=}& \frac{K^2}{{(R^2+Z^2)}^3}\{\left(\frac{2C_2}{\mu }\right)(2R^3-6R{Z}^2)+\frac{1}{2}(-R^3+3R{Z}^2)\},\hfill \\ \multicolumn{1}{c}{w_{2p}=}& \frac{K^2}{{(R^2+Z^2)}^3}\{\left(\frac{2C_2}{\mu }\right)(-4Z^3+4R^{2}Z)-2R^{2}Z\},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rr}}^{(2)p}=}& K^2\mu \{\left(\frac{2C_2}{\mu }\right)\left(\frac{-14R^4+68R^2{Z}^2-2Z^4}{{(R^2+Z^2)}^4}\right)+\left(\frac{3R^4-18R^2{Z}^2+3Z^4}{{(R^2+Z^2)}^4}\right)+\left(\frac{R^6+Z^6}{{(R^2+Z^2)}^5}\right)\},\hfill \\ \multicolumn{1}{c}{T_{\theta \theta }^{(2)p}=}& K^2\mu \{\left(\frac{2C_2}{\mu }\right)\left(\frac{2R^2-2Z^2}{{(R^2+Z^2)}^3}\right)+\left(\frac{-R^2+4Z^2}{{(R^2+Z^2)}^3}\right)\},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{zz}}^{(2)p}=}& K^2\mu \{\left(\frac{2C_2}{\mu }\right)\left(\frac{6R^4-56R^2{Z}^2+22Z^4}{{(R^2+Z^2)}^4}\right)+\left(\frac{-4R^4+20R^2{Z}^2}{{(R^2+Z^2)}^4}\right)+\left(\frac{R^6+9R^4{Z}^2+12R^2{Z}^4+4Z^6}{{(R^2+Z^2)}^5}\right)\},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rz}}^{(2)p}=}& K^2\mu \{\left(\frac{2C_2}{\mu }\right)\left(\frac{44Z^{3}R-40Z{R}^3}{(R^2+Z^2){)}^4}\right)+\left(\frac{17Z{R}^3-10Z^{3}R}{{(R^2+Z^2)}^4}\right)\}.\hfill \end{array}$$

(39)

where:

$$K=\frac{P_K}{8\pi \mu }.$$

(40)

The displacement and stress components (39) satisfy the first- and second-order incompressibility conditions in equation (6) and the second-order equations of equilibrium (2), where:

$$\begin{array}{cc}\multicolumn{1}{c}{G(R,Z)=}& K^2\left(\frac{-R^2+3Z^2}{{(R^2+Z^2)}^3}\right),\hfill \\ \multicolumn{1}{c}{H_R(R,Z)=}& K^2\mu \left(\frac{-42Z^{2}R+6R^3}{{(R^2+Z^2)}^4}\right);H_Z(R,Z)=K^2\mu \left(\frac{-36Z^3+12Z{R}^2}{{(R^2+Z^2)}^4}\right),\hfill \end{array}$$

(41)

and

$$-p_2(R,Z)=K^2\mu \left\{\left(\frac{2C_2}{\mu }\right)\left(\frac{14Z^2-2R^2}{{(R^2+Z^2)}^3}\right)\right\}.$$

(42)

We note that the second-order solution is symmetric about $Z=0$ , with $w_{2p}(R,0)=0$ and $T_{\mathrm{rz}}^{(2)p}(R,0)=0$ . Therefore, on any closed surface enclosing the origin, the resultant traction in the Z-direction is zero, and there is no second-order contribution to $P_K$ .

It may be noted that the order of any singularities that may exist in the first-order solution [130–132] is usually increased in the second-order solution (England [133]). The second-order plane strain problem of an infinite plane containing an elliptic cavity was examined by Lianis [134] using the complex variable-based methodology presented by Adkins et al. [71]. The study was extended by De Hoff [135] to include the problem of an infinite plane containing an elliptical hole and subjected to an oblique uniaxial far field tension. The article by Varley and Cumberbatch [136] examines the elliptical cavity problem for a Fritz John-type hyperelastic material from which the stress magnification for a crack can be recovered. Both analyses indicate that in the limit of the elliptical cavity degenerating to a crack, the order of the singularity is increased for the second-order solution. Similar conclusions are presented by Knesl and Semela [137] for the case of the second-order results for a Westergaard-type crack. The finite deformation of a crack into the shape of a wedge is also treated by Blatz [138] and the stress states corresponding to the infinitesimal and finite strain solutions are compared. A related problem of inhomogeneous deformations of hyperelastic wedges is discussed by Rajagopal and Carroll [139]. A useful discussion of the limits of successive approximate techniques similar to that encountered in the Poincare–Lighthill–Kuo method in aerofoil theory was presented by Tsien [140]. The mathematical modelling of cracks in hyperelastic materials was examined in a series of elegant articles by Knowles and Sternberg [141–145] and Knowles [146]. Complimentary investigations of the plane crack problem were also discussed by Amazigo [147], Lo [148], Le and Stumpf [149], Geubelle and Knauss [150] and Tarantino [151]. Finite deformation analyses at the crack tip located at a bimaterial elastic interface consisting of hyperelastic materials were also developed by Knowles and Sternberg [152]. Ravichandran and Knauss [153], Herrmann [154, 155], Geubelle and Knauss [156, 157], Gao and Shi [158] and Gao [159] to investigate the oscillatory nature of the singularity at the tip of a crack located in a bimaterial region under far-field stresses. The vanishing of the oscillatory stress singularity can result from either the incompressibility constraint usually associated with hyperelastic materials and/or the finite curvature that the crack tip can acquire from the finite deformation.

4.1. Solutions of the homogeneous equation

Considering the homogeneous equation (37), expressed in spherical polar coordinates, we can write:

$$E^4{\mathrm{\Psi }}_{2h}=(\frac{{\partial }^2}{\partial S^2}+\frac{1}{S^2}\frac{{\partial }^2}{\partial {\mathrm{\Theta }}^2}-\frac{\cot \mathrm{\Theta }}{S^2}\frac{\partial }{\partial \mathrm{\Theta }})(\frac{{\partial }^2}{\partial S^2}+\frac{1}{S^2}\frac{{\partial }^2}{\partial {\mathrm{\Theta }}^2}-\frac{\cot \mathrm{\Theta }}{S^2}\frac{\partial }{\partial \mathrm{\Theta }}){\mathrm{\Psi }}_{2h}=0.$$

(43)

Considering the substitution $\rho =\cos \mathrm{\Theta }$ , equation (40) can be written as:

$$(\frac{{\partial }^2}{\partial S^2}+\frac{(1-{\rho }^2)}{S^2}\frac{{\partial }^2}{\partial {\rho }^2})(\frac{{\partial }^2{\mathrm{\Psi }}_{2p}}{\partial S^2}+\frac{(1-{\rho }^2)}{S^2}\frac{{\partial }^2{\mathrm{\Psi }}_{2p}}{\partial {\rho }^2})=0.$$

(44)

We seek solutions of equation (44) of the form:

$${\mathrm{\Psi }}_{2p}=\frac{f\left(\rho \right)}{S},$$

(45)

which can be used to reduce equation (44) to:

$$((1-{\rho }^2)\frac{d^2}{d{\rho }^2}+12)((1-{\rho }^2)\frac{d^2}{d{\rho }^2}+2)f\left(\rho \right)=0.$$

(46)

We note that the operators in equation (46) commute, i.e.,

$$((1-{\rho }^2)\frac{d^2}{d{\rho }^2}+2)((1-{\rho }^2)\frac{d^2}{d{\rho }^2}+12)f\left(\rho \right)=0.$$

(47)

Solving the resulting ordinary differential equations (ODEs), we obtain the following homogeneous solutions:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Psi }}_{2h}=A\left(\frac{R^2}{{(R^2+Z^2)}^{3/2}}\right)+B(\frac{2Z}{(R^2+Z^2)}+\frac{R^2}{{(R^2+Z^2)}^{3/2}}{\log }_e\left(\frac{\sqrt{R^2+Z^2}+Z}{\sqrt{R^2+Z^2}-Z}\right))+C\left(\frac{R^2(R^2-4Z^2)}{{(R^2+Z^2)}^{5/2}}\right)}\\ +D(\frac{26Z{R}^2-4Z^3}{{(R^2+Z^2)}^2}+\frac{3R^2(R^2-4Z^2)}{{(R^2+Z^2)}^{5/2}}{\log }_e\left(\frac{\sqrt{R^2+Z^2}+Z}{\sqrt{R^2+Z^2}-Z}\right)),\hfill \end{array}$$

(48)

where $A,B,C$ , and $D$ are the arbitrary constants. Evaluating the displacement and stress components from equation (48) we note that, in order to maintain the second-order displacement component $u_{2h}$ and the second-order shear stress component $T_{\mathrm{rz}}^{(2)h}$ to be bounded as $R\to 0$ , we require $B=D=0$ . The remaining solutions give the following second-order displacement and stress components:

$$\begin{array}{c}{u}_{2h}=A(-\frac{3\mathrm{RZ}}{{(R^2+Z^2)}^{5/2}})+C\left(\frac{2R^{3}Z-3R{Z}^3}{{(R^2+Z^2)}^{7/2}}\right),\hfill \\ w_{2h}=A\left(\frac{R^2-2Z^2}{{(R^2+Z^2)}^{5/2}}\right)+C\left(\frac{-2Z^4+3Z^2{R}^2}{{(R^2+Z^2)}^{7/2}}\right),\hfill \end{array}$$

(49)

and

$$\begin{array}{c}{T}_{\mathrm{rr}}^{(2)h}=A\left(\frac{24R^{2}Z-6Z^3}{{(R^2+Z^2)}^{7/2}}\right)+C\left(\frac{-22R^{4}Z+46Z^3{R}^2-2Z^5}{{(R^2+Z^2)}^{9/2}}\right),\hfill \\ T_{\theta \theta }^{(2)h}=A\left(\frac{6R^{2}Z-6Z^3}{{(R^2+Z^2)}^{7/2}}\right)+C\left(\frac{-2R^{4}Z-4Z^3{R}^2-2Z^5}{{(R^2+Z^2)}^{9/2}}\right),\hfill \\ T_{\mathrm{zz}}^{(2)h}=A\left(\frac{-18R^{2}Z+12Z^3}{{(R^2+Z^2)}^{7/2}}\right)+C\left(\frac{6R^{4}Z-48Z^3{R}^2+16Z^5}{{(R^2+Z^2)}^{9/2}}\right),\hfill \\ T_{\mathrm{rz}}^{(2)h}=A\left(\frac{-6R^3+24Z^{2}R}{{(R^2+Z^2)}^{7/2}}\right)+C\left(\frac{2R^5-36Z^2{R}^3+32Z^{4}R}{{(R^2+Z^2)}^{9/2}}\right).\hfill \end{array}$$

(50)

We observe that since $u_{2h}(R,0)=0$ and $T_{\mathrm{zz}}^{(2)h}(R,0)=0$ , the homogeneous solution will involve only states of deformation and stress that are asymmetric about $Z=0$ , which can be caused only by resultant forces similar in character to the Kelvin force and this includes higher-order singularities in both the displacement and stress fields. The first-order displacement function for a concentrated force acting at the interior of an elastic infinite space of infinite extent is given by:

$${\left({\mathrm{\Psi }}_1\right)}_{\mathrm{Kelvin}\mathrm{force}}=C_1^{*}\left(\frac{R^2}{{(R^2+Z^2)}^{1/2}}\right),$$

(51)

where $C_1^{*}$ is a constant. The first-order displacement potential for a doublet (see section 5), which consist of two collinear Kelvin forces separated by a small distance $d$ (Figure 2), is given by the partial derivative of equation (51) with respect to $Z$ , i.e.,

$${\left({\mathrm{\Psi }}_1\right)}_{\mathrm{Doublet}}=C_2^{*}\left(\frac{R^{2}Z}{{(R^2+Z^2)}^{3/2}}\right),$$

(52)

where $C_2^{*}$ is a constant. The first-order displacement potential for a combination of doublets (i.e., a doublet in tension in close proximity to a doublet in compression) can be obtained by taking the partial derivative of equation (52) with respect to $Z$ , i.e.,

$${\left({\mathrm{\Psi }}_1\right)}_{\mathrm{Doublet}}=C_3^{*}(\frac{R^2}{{(R^2+Z^2)}^{3/2}}-\frac{3R^2{Z}^2}{{(R^2+Z^2)}^{5/2}}),$$

(53)

where $C_3^{*}$ is a constant. The result (53) has a form similar to a homogeneous solution that can be obtained by suitably adjusting the coefficients $A$ and $C$ in the reduced solution for equation (48). Therefore, the second-order solutions derived from the homogeneous solution (45) become relevant only if additional second-order local forces are applied in the vicinity of the Kelvin force. It can be concluded that the second-order displacement and stress fields corresponding to the Kelvin force problem are obtained from the particular solution of equation (37). These solutions agree with the results obtained by Carroll and Rooney [94] for the second-order Kelvin’s problem, using a Love Strain potential approach along with Signorini’s compatibility condition.

Figure 2.

The doublet problem.

5. The doublet problem

Kelvin’s classical elasticity solution for the single point force acting at the interior of an elastic solid of infinite extent can be successfully used to find novel solutions of the equations of classical elasticity corresponding to other singularities. One such problem, described by Love [96] as a “double force without moment,” is obtained by superposing on the Kelvin force problem (the force denoted by $P_D$ ) a point force of equal magnitude, acting in the opposite direction but located at a small distance $d$ from the origin, in the negative Z-direction (Figure 2). Since the second force is acting in the opposite direction and considering $d$ to be infinitesimally small, the first-order solutions to the doublet problem is obtained by replacing any term in the solution to Kelvin’s problem, e.g., $f(R,Z)$ by $f+(\partial f/\partial Z)d$ (see, e.g., [99]).

For an incompressible elastic material, the first-order displacement and stress components for the doublet problem can be written as:

$$\begin{array}{cc}\multicolumn{1}{c}{u_1=}& \tilde{C}\frac{(-R^3+2R{Z}^2)}{{(R^2+Z^2)}^{5/2}};w_1=\tilde{C}\frac{(2Z^3-Z{R}^2)}{{(R^2+Z^2)}^{5/2}},\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rr}}^{(1)}=}& 2\tilde{C}\mu \left(\frac{3R^4-12R^2{Z}^2}{{(R^2+Z^2)}^{7/2}}\right);T_{\theta \theta }^{(1)}=0,\hfill \\ \multicolumn{1}{c}{T_{\mathrm{zz}}^{(1)}=}& 2\tilde{C}\mu \left(\frac{9R^2{Z}^2-6Z^4}{{(R^2+Z^2)}^{7/2}}\right);T_{\mathrm{rz}}^{(1)}=2\tilde{C}\mu \left(\frac{6R^{3}Z-9R{Z}^3}{{(R^2+Z^2)}^{5/2}}\right),\hfill \end{array}$$

(54)

where:

$$\tilde{C}=\frac{P_{D}d}{8\pi \mu }.$$

(55)

Omitting details, the partial differential equation governing the second-order displacement function ${\mathrm{\Psi }}_2(R,Z)$ takes the form:

$$E^4{\mathrm{\Psi }}_2=\frac{{\tilde{C}}^2}{{(R^2+Z^2)}^6}[144\left(\frac{2C_2}{\mu }\right)(R^{4}Z-4R^2{Z}^3)-27(3R^{4}Z-2R^2{Z}^2+5Z^5)].$$

(56)

The second-order expressions for the displacement and stress components obtained from the particular solution of equation (56) can be written in the forms:

$$\begin{array}{cc}\multicolumn{1}{c}{u_{2p}=}& \frac{{\tilde{C}}^2}{{(R^2+Z^2)}^5}[\left(\frac{2C_2}{\mu }\right)(-5R^5+47R^3{Z}^2-20R{Z}^4)-\left(\frac{3}{8}(-4R^5+52R^3{Z}^2-16R{Z}^4)\right)]\hfill \\ \multicolumn{1}{c}{w_{2p}=}& \frac{{\tilde{C}}^2}{{(R^2+Z^2)}^5}[\left(\frac{2C_2}{\mu }\right)(-20R^{4}Z+44R^2{Z}^3-8Z^5)-\left(\frac{3}{8}(-24R^{4}Z+48R^2{Z}^3)\right)]\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rr}}^{(2)p}=}& \frac{{\tilde{C}}^2}{{(R^2+Z^2)}^6}[2C_2(48R^6-684R^4{Z}^2+528R^2{Z}^4)+\mu (-12R^6+324R^4{Z}^2-237R^2{Z}^4+12Z^6)]\hfill \\ \multicolumn{1}{c}{T_{\theta \theta }^{(2)p}=}& \frac{{\tilde{C}}^2}{{(R^2+Z^2)}^6}[2C_2(72R^4{Z}^2+72R^2{Z}^4)+\mu (3R^6-36R^4{Z}^2-27R^2{Z}^4+12Z^6)]\hfill \\ \multicolumn{1}{c}{T_{\mathrm{zz}}^{(2)p}=}& \frac{{\tilde{C}}^2}{{(R^2+Z^2)}^6}[2C_2(-30R^6+540R^4{Z}^2-618R^2{Z}^4+72Z^6)+\mu (18R^6-279R^4{Z}^2+300R^2{Z}^4+12Z^6)]\hfill \\ \multicolumn{1}{c}{T_{\mathrm{rz}}^{(2)p}=}& \frac{{\tilde{C}}^2}{{(R^2+Z^2)}^6}[2C_2(240R^{5}Z-80R^3{Z}^3+216R{Z}^5)+\mu (-111R^{5}Z+402R^3{Z}^3-72R{Z}^5)].\hfill \end{array}$$

(57)

The second-order displacement and stress components (57) satisfy the second-order incompressibility condition (6) and the second-order equations of equilibrium (23), where:

$$\begin{array}{cc}\multicolumn{1}{c}{G(R,Z)=}& {\tilde{C}}^2\left(\frac{3R^4-12R^2{Z}^2+12Z^4}{{(R^2+Z^2)}^5}\right),\hfill \\ \multicolumn{1}{c}{H_R(R,Z)=}& {\tilde{C}}^2\left(\frac{-54R^5+162R^3{Z}^2-324R{Z}^4}{{(R^2+Z^2)}^6}\right),\hfill \\ \multicolumn{1}{c}{H_Z(R,Z)=}& {\tilde{C}}^2\left(\frac{-108R^{4}Z+216R^2{Z}^3-216Z^5}{{(R^2+Z^2)}^4}\right).\hfill \end{array}$$

(58)

6. Boussinesq’s problem

The solution to the problem of an isotropic elastic half-space that is subjected to a concentrated normal force at the origin of coordinates is identified as Boussinesq’s problem (Boussinesq [160]) (Figure 3). The solution to the linear elasticity problem can be obtained in diverse ways, including through the application of integral transform techniques [98–105, 161–163] and through the successive applications of superposition techniques involving Kelvin’s solution and Love’s strain potential technique [96]. The recent study [113] develops a solution to Boussinesq’s problem using Kelvin’s solution and dimensional considerations in choosing additional solutions required to render the surface of the half-space region traction free. While methods have been developed to apply Boussinesq’s solution to examine the decay in the distribution of stress in earth masses subjected to surface loads, the modifications proposed fail to satisfy the requirement of classical elasticity, particularly in terms of satisfying the Beltrami–Michell compatibility conditions [114].

Figure 3.

Boussinesq’s problem.

The application of the theory of second-order elasticity to Boussinesq’s problem is a natural extension of the study dealing with Kelvin’s problem. The first-order solution for Boussinesq’s problem is identical in form to Kelvin’s solution given by equation (35) and the arbitrary constant corresponding to equation (36) is given by:

$$B^{*}=\frac{P_B}{4\pi }.$$

(59)

The stress field given in equation (39) must satisfy the traction-free boundary conditions applicable to the deformed surface of the half-space region. The equation of the deformed surface of the half-space valid to order $\epsilon$ is given by:

$$F(r,z)=z-\epsilon w_1(R,Z)=\mathrm{const}.$$

(60)

The components of the unit normal in the r- and z-directions are given by:

$$\mathbf{n}=\frac{(-\epsilon \frac{\partial w_1}{\partial r},0,\{1-\epsilon \frac{\partial w_1}{\partial z}\})}{\sqrt{{\epsilon }^2{\left(\frac{\partial w_1}{\partial r}\right)}^2+{\epsilon }^2{(1-\epsilon \frac{\partial w_1}{\partial z})}^2}},$$

(61)

which reduces to:

$$\mathbf{n}=(-\epsilon \frac{\partial w_1}{\partial R},0,1).$$

(62)

The traction-free boundary conditions for the second-order stress components on the deformed surface of the half-space reduce to:

$$\begin{array}{c}{T}_{\mathrm{rz}}^{(2)}-\frac{\partial w_1}{\partial R}{T}_{\mathrm{rr}}^{(1)}=0\mathrm{on}Z=0,\hfill \\ T_{\mathrm{zz}}^{(2)}-\frac{\partial w_1}{\partial R}{T}_{\mathrm{rz}}^{(1)}=0\mathrm{on}Z=0.\hfill \end{array}$$

(63)

Since $T_{\mathrm{rr}}^{(1)}=T_{\mathrm{rz}}^{(1)}=0\mathrm{on}Z=0$ , the second-order stress field should satisfy the boundary conditions $T_{\mathrm{rz}}^{(2)}=T_{\mathrm{zz}}^{(2)}=0\mathrm{on}Z=0$ . The results from the particular solution (39) give:

$$\begin{array}{c}{T}_{\mathrm{rz}}^{(2)}=0\mathrm{on}Z=0,\hfill \\ T_{\mathrm{zz}}^{(2)}=\frac{3K^2}{R^4}(4C_2-\mu )\mathrm{on}Z=0.\hfill \end{array}$$

(64)

The first boundary condition of equation (64) is identically satisfied by the expression for $T_{\mathrm{rz}}^{(2)p}$ given by equation (39). The second boundary condition needs to be satisfied by selecting a suitable solution of the homogeneous equation of equation (37) that will have the correct distribution of $T_{\mathrm{zz}}^{(2)h}(R,0)$ . The reduced form of equation (48) will not provide the necessary solution. The alternative is to use Love’s stress function approach and obtain the stress field that can satisfy the second boundary conditions of equation (64). Selvadurai [86] showed that a solution based on the Hankel integral transforms (Lamb [161], Sneddon [162,163], Terezawa [164]) does not provide an approach to satisfy the required traction-free boundary condition. A direct integration of Boussinesq’s solution applicable to a half-space is a possible approach. For example, the expression for the required homogeneous solution for the displacement $u_{2h}(R,Z)$ will give results of the type

$$u_{2h}(R,Z)=\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\infty }{\int }}\frac{\mathrm{CZ}}{8\pi \mu {(R\prime )}^4}\frac{{\{R^2+{(R\prime )}^2-2\mathrm{RR}\cos \theta \prime \}}^{1/2}}{{\{Z^2+R^2+{(R\prime )}^2-2\mathrm{RR}\prime \cos \theta \prime \}}^{3/2}}R\prime \mathrm{dR}\prime d\theta \prime .$$

(65)

While this approach is feasible, the divergent integrals can only be evaluated numerically, which is a limitation. Also, consideration of a hemi-spherical inclusion at the surface of the half-space results in the first-order solution that is a series in the Legendre polynomials, which cannot be manipulated easily in constructing the second-order formulation. Subsequently, the second-order solution to Boussinesq’s problem for a concentrated normal force was provided in a compact and closed form by Carroll and Rooney [95], using the strain potential approach of Love [96].

7. Conclusion

The theory of second-order elasticity is a mathematically consistent theory for modeling hyperelastic materials undergoing moderately large strains, as opposed to either media exhibiting large deflections and rotations but with small strains or small deformations superposed on large. The expansion of the dependent variables in terms of a small non-dimensional parameter forms the basis of the approach and if a particular problem is formulated in a consistent manner, the small parameter will naturally evolve in the formulation. The solution of problems in second-order elasticity theory also becomes meaningful in terms of technological applications involving rubber-like elastic materials used as load bearing components or mountings. Also, when the incompressibility constraint is introduced, the strain energy function of the Mooney–Rivlin form can completely accommodate the second-order moderately large strain phenomena. The formulation and solution of problems in second-order elasticity can be facilitated through the use of complex potentials and complex variable theory, stress functions based on the adaptation of the functions proposed by Love, Neuber–Papkovich, and integral transform techniques. Spencer’s displacement function approach is also a convenient method for the treatment of the second-order incompressible elastic problem and reduces the governing partial differential equation for the function to a canonical linear form. A similar result is obtained for the second-order component of the isotropic stress. The displacement function approach is applied to develop solutions to certain axisymmetric localized loading problems of an infinite space and a half-space. The first-order problem corresponds to the classical elasticity solution and the order of the stress singularities present in the first-order problems are generally increased in the second-order problem. The issue of singularities in the first-order solution can be alleviated by selecting suitable boundary conditions where the loads are applied over finite regions. The ensuing formulations can involve algebraic complexity. In such cases the solution of the second-order problem entails a great deal of routine mathematical operations and computer-aided symbolic mathematical manipulation techniques can be used to solve the second-order problem with speed and accuracy.