In this paper, we consider the problem of an infinite isotropic elastic matrix that contains a material surface described by the Steigmann-Ogden theory. The theory treats the surface as a shell of vanishing thickness that possesses a boundary and is characterized by membrane and bending stiffness as well as residual surface tension. The Steigmann-Ogden theory could be useful for modeling emerging class of nano-and bio-materials that use ultrathin nanoplatelets or membranes (e.g., graphene or graphene-oxide sheets) as reinforcements, e.g., [6,11,28,29,32]. While mechanical properties of the nanoplatelets/sheets and the overall elastic properties of this class of composites have been investigated by several research groups, until recently only a few attempts have been made to study their local mechanical behavior, see [36]. In the last few years, it was suggested to use the Gurtin-Murdoch theory [12,13] for modeling materials with flexible and extensible reinforcements; the corresponding models have been investigated in [2] for the anti-plane setting and in [26] for the plane strain setting. The choice of the Gurtin-Murdoch model was motivated by the following facts: i) the graphene lattice has a hexagonal structure and, therefore, the material surface behaves as an isotropic elastic material under small deformation, see [6,34,35], ii) in many practical applications, the assumption of isotropic linear elasticity also holds for the matrix materials.
As the material surface in the Gurtin-Murdoch theory does not include bending resistance, the theory might not be appropriate for modeling of some reinforcements, as theoretical studies, see e.g. [16,18,21], have suggested that bending stiffness of graphene might, in some cases, be critical. In order to enable such modeling, one can use the extension of the Gurtin-Murdoch theory that was proposed by Steigmann and Ogden in [30,31]. Instead of a membrane, the Steigmann-Ogden theory treats the material surface as a shell of vanishing thickness.
In this paper, we present, for the first time, numerical solution of the two-dimensional plane-strain Steigmann-Ogden model and discuss its applications for simulating in-plane behavior of the composites with ultra-thin and stiff reinforcements. The model includes an isotropic elastic matrix that contains a material surface described by the Steigmann-Ogden theory, where the material surface possesses a boundary. We assume that the volume fraction of reinforcements is small (which is typically the case in graphene-based composites) and model the matrix material as an infinite isotropic elastic medium.
First, we use the theories of elastic layer potentials and integral equations and provide the governing equations for the model. Exact expressions for the elastic fields everywhere in the material are provided in terms of integral representations written in terms of complex variables. We then particularize those expressions for the case of a surface along a straight segment. Using the approximations for the boundary data that involve the series of Chebychev's polynomials and allow for automatic incorporation of tip conditions, we obtain the infinite system of linear algebraic equations for the unknown series coefficients. The elastic fields in the matrix are then found using appropriate complex integral representations.
We then use the obtained solution to illustrate the influence of governing dimensionless parameters and to study the local mechanical behavior of composite materials reinforced with graphene-based nanoplatelets.
The paper is structured as follows. In Section 2, we review the governing equations, jump, and tip conditions of the Steigmann-Ogden model for general three-dimensional case. We then adopt the plane strain assumption, particularize the general equations to meet that assumption, formulate the problem under study, and describe the connection of the problem with that of inextensible shell type elastic inhomogeneity/layer. In Section 3, we present integral representations for the fields everywhere in the material system, introduce complex variables combinations for the problem geometry and those for the elastic fields. The case of a material surface along a straight segment is considered in more details in Section 4. There, we introduce the dimensionless variables that govern the problem, reformulate the system of boundary integral equation in the dimensionless setting, discuss the uniqueness and existence of its solution, and describe the numerical solution technique. In Section 5, we present several examples of numerical simulations. Concluding remarks are presented in Section 6.
The Steigmann-Ogden theory describes an infinite isotropic elastic medium that contains a material surface treated as a shell of vanishing thickness and characterized by membrane and bending stiffness as well as residual surface tension. The governing equations of the theory obtained in [30,31], and also reported in e.g., [8,9,25,37], are briefly reviewed below.
The governing equations for the matrix material are those of standard linear elasticity. In particular, the displacement vector u satisfies the Navier equation that, for the case of vanishing body forces, has the following form:
where λ and µ are the Lamé parameters, the symbol "•" identifies the dot product of two vectors, ∇ is the gradient operator, and ∇ 2 is the Laplacian. Eq. ( 1) is supplemented by the following conditions describing the behavior of the elastic fields at the material surface S, across S, and at its boundary δS: a) Continuity of the displacements across S
where the superscripts "+", "-" here and below describe the limit values of the fields when S is approached from the direction of that of the normal vector or from the opposite direction, respectively. b) Traction jump conditions across
where n is a unit vector normal to S, σ is the Cauchy stress tensor of the matrix, ∆σ = σ + -σ -is the jump in stresses across S, T represents the (first) Piola-Kirchhoff surface stress tensor, M is the surface couple stress tensor, ∇ S = ∇ -n (n • ∇) is the surface gradient operator. c) Constitutive equations at S
in which I S = In ⊗ n is the unit tangent tensor, I is the three-dimensional unit tensor, µ S and λ S are the surface Lamé parameters, σ 0 is the residual surface tension, ζ S and χ S are the surface bending stiffness parameters, and "tr" is the trace operator. The surface tensors ε S and κ S are given by the formulas:
where the subscript "T " identifies transposition and
d) Conditions at the boundary ∂S
where ν, τ are respectively the tangent normal and tangent unit vectors of the Darboux frame for the curve ∂S and s is the arc length of ∂S.
Consider a cylindrical surface in an infinite isotropic elastic matrix whose sufficiently smooth cross-section represents an open arc L that is located inside Ox 1 x 2 plane of the global Cartesian coordinate system with the unit basis vectors e 1 , e 2 , e 3 . Introduce the local coordinates with the mutually orthogonal unit vectors n, ℓ, e 3 shown on Fig. 1 and treat s as the arc length parameter for L.
Under plane strain assumptions, the only nonzero components of the displacement vector u are u 1 (x 1 , x 2 ) , u 2 (x 1 , x 2 ) and properly modified Eq. (1) remains valid. The displacement vector can be represented in the local coordinate system as
where u n , u ℓ are its normal and tangential components. Thus, the expressions for the surface operators I S and ∇ S become
Using Eqs. ( 10)-( 12), one can rewrite constitutive equations ( 4)-( 7) as, see [25,26] for the details,
with
)
and R = R(s) being the local radius of curvature on L. Using Eqs. ( 10)-( 17), one can rewrite Eqs. ( 2)-( 3) as, see [25,26] for the details,
where σ n and σ ℓ are the normal and shear tractions of the bulk material.
Taking into account that the vectors ν, τ of Eqs. ( 8), ( 9) are given at the end-points of L by ν = ±ℓ (pointing out of L) and τ = e 3 and using Eqs. ( 18), (19), we obtain the following conditions at the tips ξ= a and ξ= b of L :
Thus, the boundary value problem of plane strain elasticity for the model considered is governed by classical equation for the matrix material, comple-mented by Eqs. ( 18)-( 21).
It is known (see [3,14,25,37]) that, in the case of zero surface tension, the jump conditions in the Steigmann-Ogden model represent specific case of the conditions that appear in elasticity problems containing the so-called inextensible shell type inhomogeneity/layer, see [3,4,5]. In such case, the thickness h of the inhomogeneity/layer and its properties must behave as
where µ I (λ I ) is the shear modulus (Lamé parameter) of the inhomogeneity/layer, H is some characteristic length-scale of the problem. The big-Θ symbol means an asymptotically tight bound. The elastic parameters of the Steigmann-Ogden model are related to those for the inhomogeneity/layer as, see e.g. [25,37],
Boundary conditions Eqs. ( 18)-( 20) are naturally satisfied, if the displacements in the matrix are sought in the form of a single layer elastic potential.
Using the theory of elastic potentials, introduced in [17], it was shown in [23,24] that the displacements everywhere in the domain can be expressed via the traction jumps as
where the repeated index implies summation, u ∞ k (x) are the components of the displacements in the homogeneous plane (without the material surface) due to the far-field load, t 1 , t 2 are the components of the two-dimensional traction vector t in the Cartesian coordinate system, and G kj (x, ξ) are the following displacements in the k -th direction at point x due to a unit concentrated force in the j -th direction at point ξ:
in which δ kj is Kronecker's symbol,
Complex variables form of the above representation, obtained in [19,20,24], can be written in terms of complex displacements u = u 1 + iu 2 and complex tractions that, in the Cartesian coordinates, are introduced as (here and in some expressions below we omitted the arguments, for brevity)
and related to the complex tractions in the local coordinates σ = σ n + iσ ℓ as
where α = α(s) is the angle between the axis Ox 1 and the tangent at the point ξ ∈ L.
Taking into account that
where τ = ξ 1 + iξ 2 , we obtain
and performing some algebraic manipulations, see [19,20,24], one can arrive at the following integral representation for the complex displacements outside of L:
in which z = x 1 +ix 2 is the complex combination of the Cartesian coordinates of the point z / ∈ L,
and, assuming the case of uniform far-field stresses σ
The expressions for the complex tractions on some lines outside of L can be obtained as, see [20,24,26]
where
In the following, we will need the equation for u ′ (z), which is
where
The boundary integral equations resulting from the limiting procedure in which the field point is allowed to reach the boundary point τ 0 = τ 0 1 + iτ 0 2 have the following forms, see [26]:
Using Eqs. ( 19), (20), one can express complex density ∆σ(z) via a com-bination of various derivatives of complex displacement u(z) as
Thus, Eq. ( 39), in combination with Eq. (40) provides the governing complex boundary integro-differential equation in terms of boundary displacements. It must be supplemented by the following conditions at the tips ξ= a and ξ= b, see Eqs. (21):
In the following we will require that the real and imaginary parts of u ′ (z) are the functions of the space C 1,α (L) with 0 < α < 1 at any point of L, except for its tips. This requirement, in combination with Eq. (41), guarantees that u ′ ∈ h 2 class on L, see [27] for the definition. Thus, the solution of Eq. (39) should be sought in that class of functions.
Assume that L is straight segment [-a, a] located on the axis Ox 1 , Fig. 2. Then, omitting the arguments of the functions for brevity, we arrive at and Eqs. (39), (40) are reduced to
Instead of directly solving the system of Eqs. ( 46), (47), let us first construct and analyze the integral representations for σ S and ω S . We note that from Eqs. (21) it follows that ε S , σ S ∈ h 2 class on L.
Using Eqs. ( 46), (47), and the fact that τ, τ 0 are real variables, one gets
where
As σ S vanishes at the tips, singular integral equation ( 48) is equivalent to the following hypersingular integral equation:
where symbol = denotes the Hadamard integral (the finite part of hypersingular integral).
Eqs. ( 48)-( 50) are identical to the corresponding equations for the Gurtin-Murdoch material surface, see Eqs. ( 62)-(64) in [26]. However, this is not the case for the corresponding expressions related to ω S . For example, the counterpart of Eq. ( 48) can be obtained using Eqs. ( 46), (47) as
As the expression σ 0 ω S -(2χ S +ζ S ) d 2 ω S dτ 2 also vanishes at the tips, singular integral equation ( 51) is equivalent to the following hypersingular integral equation, which is the counterpart of Eq. ( 50):
With the use of the following dimensionless variables (k, j = 1, 2):
one can rewrite Eqs. ( 50) and (53) as (-1 ≤ τ 0 ≤ 1)
One can see that Eqs. (55) and (56) are hypersingular equations of the same type as those considered in e.g. [2,26] and they are well-studied in the literature, see the references in the latter papers. Their singular integro-differential counterparts, Eqs. ( 48) and (51) are equations of the type used in the theory of aircraft wings of finite span. Thus, it can be proved, in the same manner as it was done in Section 5.1 of [2], that Eqs. ( 55) and ( 56) have unique solutions in class h 2 class on [-1, 1]. In the case of Σ 2 = 0, Eq. ( 56) is homogeneous equation that only has trivial solution ω S = 0. This can be proved in the same manner as it was done in Section 5.1 of [2] for the homogeneous equation of similar type.
Noteworthy, Eqs. ( 55) and ( 56) are uncoupled hypersingular integral equations for σ S and ω S , respectively. Those must be complemented by tip conditions (21) rewritten as:
After Eqs. ( 55) and ( 56) are solved, the values of ∂u ℓ ∂s = ε S and ∂un ∂s = ω S can be found from ( 15)-( 17) which, after an elementary integration will lead to the values (up to a rigid body translation) of u ℓ and u n .
When γ 1 = 0, Eq. ( 55) has the following analytical solution for the dimensionless surface stress:
As it was noted in [26], in the case when γ 1 = 0 and σ 0 = 0, Eq. ( 55) describes the problem of a rigid stiffener (anticrack) in an elastic medium, see e.g. [1,7,10,22,33] and the references therein. In the case when γ 1 = ∞, it is easy to show that Eq. ( 55) has the simple solution σ S = σ 0 . Under the condition γ 2 = 0, Eq. ( 56) has the solution ω S = -Σ 2 /µ. If γ 2 = 0 and σ0 ̸ = 0, then the equation (56) reduces to the corresponding equation for the Gurtin-Murdoch theory, see [26].
As Eq. ( 55) is identical to the one obtained in [26] for the case of the Gurtin-Murdoch membrane, the unknown function σ S ( τ ) can be approximated as
where U m ( τ ) is the Chebyshev polynomial of the second kind of degree m, A m are the unknown coefficients, and N is the maximal degree of the approximation.
In the case of γ 2 = 0, Eq. ( 56) reduces to the one obtained in [26] and its solution can be obtained by using the following approximation:
where C m are again the unknown coefficients.
Here we will concentrate on the solution of Eq. ( 56) for the case of γ 2 ̸ = 0, σ0 ̸ = 0 that describes the complete Steigmann-Ogden model. In such case, we will use the following approximation:
where k = σ0 /γ 2 . Notice, that the chosen approximations automatically satisfy two out of the three conditions at the tips of the surface, namely
Treating Eq. ( 63) as a differential equation with respect to the unknown function ω S , we solve it as:
where C is a real integration constant that can be found from the remaining third condition at the tips of the surface
Substituting the approximations of Eqs. ( 61), (63), and (66) into Eqs. ( 55) and ( 56), multiplying them by U j ( τ ), j = 0, ...N , and using orthogonality properties of the Chebyshev polynomials, we obtain the following system of linear algebraic equations with respect to the unknown coefficients A j , B j :
where δ j,m is the Kronecker symbol, and the constants I 1 jm , I 2 j are given by the formulas:
The inner integral in Eq. ( 70) can be re-scaled onto the segment [0, 1], and both integrals of Eqs. ( 70), (71) can then be computed using Chebyshev quadratures:
where ti and τl are the Chebyshev nodes of the first and second kind:
and M is taken to be even. The additional tip condition for dω S /dτ becomes:
Once the system of Eqs. ( 68), (69) , and ( 74) is solved and the numerical values of σ S , ω S are obtained, the real and imaginary components of the derivative of the displacements on the boundary can be found from Eqs. (43), (45). The real and imaginary components of ∆ σ can be found from Eq. (47). The stresses at any line outside of the surface can be then evaluated using properly modified Eq. (33).
This section presents the results of several numerical studies. As Eq. (55) has been studied in detail in [26], the main difference between the case of the Steigmann-Ogden model and the Gurtin-Murdoch model is in Eq. (56). Thus, here we will concentrate on the study of the latter equation. We also note that, due to the fact that approximation (63) is not appropriate for the Gurtin-Murdoch model, it is impossible to directly reduce the solution for the Steigmann-Ogden model to its counterpart for the Gurtin-Murdoch model, by taking the limit (k = ∞).
The parameters used in the current studies were chosen within the intervals reported in the literature for the nanocomposites with either epoxy or aluminum matrix. We took the elastic properties of the matrix material to be µ = 2 GPa, ν = 0.35 for epoxy and µ = 28 GPa, ν = 0.33 for aluminum. The value of a was assumed to vary within the interval [5 nm -50 nm].
For the reinforcements, the reported values of the two-dimensional elastic modulus, Poisson ratio, and pretension (that we interpret as surface tension) are [6,28,29,32]. We adopted a wider interval of variation for the prestress by assuming the following range of its dimensionless values: σ 0 ∈ [10 -4 , 10 -1 ]. The values of the stresses at infinity are taken to be σ ∞ 11 = 0, σ ∞ 12 = 1GPa, and σ ∞ 22 = 1GPa. We start with the convergence study of the algorithm for several values of γ 2 ̸ = 0 and σ0 ̸ = 0 and investigate influence of these parameters on the elastic fields everywhere in the elastic plane. Then, we perform similar studies by varying the parameters γ 2 and σ 0 . Finally, we present the numerical results for the local stresses in the composite material composed of epoxy matrix reinforced by graphene-oxide nanoplatelets.
5.1. Stability of the algorithm for several values of the parameters σ 0 , γ 2 and the number of Chebyshev polynomials N . Fig. 3 illustrates the convergence of the numerical results for ω S for various values of the parameters σ0 and γ 2 and different numbers N of Chebyshev polynomials used in approximations. Aluminum (µ = 28GPa, ν = 0.33) is assumed for the material matrix and the number of points M used in Chebyshev quadratures to compute the integrals of Eqs. (72), (73) is calculated as M = 4N . It can be seen that, for reasonable values of the parameters σ0 , γ 2 , the results converge very fast even for relatively small values of N . The variation of the parameter γ 2 here and below is achieved by varying the parameter k = σ0 /γ 2 that measures the relationship between the dimensionless values of the prestress and bending rigidity.
Analysis of Fig. 3 leads to the conclusion that the smaller the value of k is, the faster the convergence is achieved. Observe that for very large values of k, the matrix of the system of Eqs. ( 68), (69) can become ill-conditioned, which reflects the above mentioned fact that such case should be treated with the Gurtin-Murdoch model [26] rather than that of the Steigmann-Ogden.
In the numerical examples presented in the following subsections, we adopted the values N = 40, M = 160 in order to ensure accuracy of the results.
In this section we assume that the matrix material is made from epoxy with µ = 2GPa, ν = 0.35 and investigate the influence of the parameters γ 2 , σ 0 on the values of ω S and Re ∆σ.
Fig. 4 illustrates the distributions of ω S and Re ∆σ along the surface for several values of the normalized surface prestress σ 0 , when the parameter k = 10 is held constant.
The distributions of ω S and Re ∆σ along the surface for several values of k are illustrated on Fig. 5 for the case of σ 0 = 0.01. It is obvious from Figs. 4, 5 that both parameters σ 0 and k have strong influence on the distributions of stresses on the material surface. The plots tend to become more sharp for small values of the normalized surface prestress σ 0 and large values of k.
nanoplatelets. For illustration purposes, we again consider epoxy matrix with µ = 2 GPa, ν = 0.35, assume that the far-field loading is the same as in previous example, and take a = 5 nm. The 2D elastic properties of the reinforcement have been chosen in such a way that k = 0.1, σ 0 = 0.01.
The contours of normalized stresses σ kj in the domain -1.5 ≤ x/a ≤ 1.5, -0.5 ≤ y/a ≤ 0.5 are plotted on Fig. 6. It can be seen from Fig. 6a,b that the distributions of σ 11 and σ 22 are antisymmetric with respect to the two coordinate axes, while the distribution of σ 12 is symmetric (Fig. 6c). All stress components exhibit singular behavior near the tips with the concentration of either compressive stress σ 11 (Fig. 6a) or tensile stress σ 22 (Fig. 6b); the stress σ 12 undergoes jump across the material surface (Fig. 6c). It seems that the Poisson effect in the bulk (due to the far-field tension) is larger than the effect of positive surface tension σ0 on the surface, which leads to the compressive stress σ11 ahead of the surface tips.
In this paper we derived the solution of the plane strain problem involving an infinite matrix, subjected to uniform stresses at infinity, and a Steigmann-Ogden material surface with a boundary. The solution might have possible applications in modeling mechanical processes in nanocomposite materials reinforced by two-dimensional nanoplatelets. We showed that the governing equations for the problem can be reduced to the complex variables integro-differential equation in terms of complex surface displacements and the boundary (tip) conditions can be expressed via various complex derivatives of those displacements. For the case of a surface along a straight segment, we reduced the governing equations to the system of two uncoupled real variables hypersingular equations in terms of the components of surface stress tensor and introduced the dimensionless parameters that govern the problem. The system of equations was then reformulated in dimensionless setting and a numerical solution technique, based on the approximations of the boundary unknowns by the series of Chebyshev polynomials was presented. Finally, we presented the results of several numerical studies that were designed to test the algorithm, to investigate influence of governing parameters, and to demonstrate potential for applications.
The present work can be extended in several directions. First, the numerical algorithm could be modified to allow for the solutions of the problems with material surfaces of more general shapes. Second, the approach could be extended to solve the problems with multiple surfaces, which will allow for estimations of the effective properties (the only measurable properties so far) of composites reinforced by nanoplatelets. Third, the approach could be extended, in straightforward manner, to the problems involving surfaces inside the half-planes and bonded dissimilar half-planes using the complex integral equations reported in [19,20,24]. We also plan to extend the approach to three dimensions. Finally, we are working on leveraging the advantages of higher-order continuous spline approximations that have recently seen increasing popularity in the context of isogeometric analysis. Our goal is to enable numerical schemes to accommodate the higher-order derivatives that appear in the Gurtin-Murdoch and Steigmann-Ogden models. On the one side, we focus on isogeometric finite element discretizations, where we are working on formulating robust variational interface conditions that can take into account the specific Gurtin-Murdoch and Steigmann-Ogden specific continuity and jump conditions, see [15]. On the other side, we focus on the isogeometric boundary element method, where we are focusing on accounting for the multiple boundary conditions of the material surfaces via new collocation strategies.