On well-posedness of two-phase nonlocal integral models for higher-order refined shear deformation beams

Due to the conflict between equilibrium and constitutive requirements, Eringen’s strain-driven nonlocal integral model is not applicable to nanostructures of engineering interest. As an alternative, the stress-driven model has been recently developed. In this paper, for higher-order shear deformation beams, the ill-posed issue (i.e., excessive mandatory boundary conditions (BCs) cannot be met simultaneously) exists not only in strain-driven nonlocal models but also in stress-driven ones. The well-posedness of both the strain- and stress-driven two-phase nonlocal (TPN-StrainD and TPN-StressD) models is pertinently evidenced by formulating the static bending of curved beams made of functionally graded (FG) materials. The two-phase nonlocal integral constitutive relation is equivalent to a differential law equipped with two restriction conditions. By using the generalized differential quadrature method (GDQM), the coupling governing equations are solved numerically. The results show that the two-phase models can predict consistent scale-effects under different supported and loading conditions.


Introduction
Beam-like micro-/nano-structures are essential parts and are often utilized in the engineering design of micro-/nano-electromechanical systems (MEMSs/NEMSs). In such applications, both experiments and atomistic simulations show that the size-effect must be comprehensively and rigorously considered. Owing to the lack of any inherent length characteristic parameter, the classical elasticity theory fails to address such size-dependent problems. Molecular dynamics (MD) simulations require high computational costs, and experiments are often difficult to implement at a micro-/nano-scale. Therefore, several non-classical continuum theories are suggested to account for the size-effects, thereby overcoming the breakdown of classical elasticity. Among the higher-order continuum theories, the nonlocal elasticity model [1][2][3] is a popular tool.
Unlike the classical local elasticity theory, the nonlocal stress is expressed as a convolution integral of a decaying kernel and the strain fields at each point, which is the so-called averaging kernel. At the same time, a length-scale parameter is introduced to evaluate the size-effect. In practical applications [4][5][6] , integral constitutive relations are usually converted into a differential form that is easier to solve. However, some paradoxes appear when applying the differential form of Eringen's strain-driven model to bounded beam structures. In particular, Refs. [7] and [8] have shown that the results based on Eringen's model coincided with the local one in the case of a cantilever beam subjected to tip forces, which are inconsistent with the results under other supported conditions. As for a reason, Romano et al. [9] systematically studied Eringen's nonlocal differential model over the bounded continuous domains and recently pointed out that such inconsistencies are induced by the neglect of the two constraint conditions related to the differentiation process of the originally integral equation, and these two constraint conditions happen to conflict with the equilibrium requirement in the case of cantilever beams under point-loading conditions. Alternatively, the two-phase formulation of the nonlocal model suggested by Eringen is used to bypass such ill-posed issues. Wang et al. [10][11] employed such a model to study the Euler-Bernoulli and Timoshenko beams, in which the integral constitutive equation is utterly equivalent to a differential equation with two higher-order constitutive boundary conditions (BCs). Zhang et al. [12] investigated the scale-dependent bending response of circular beams. The two-phase nonlocal integral equations are solved directly via the Laplace transform method.
Another option for removing the defects mentioned above is the so-called stress-driven nonlocal integral elasticity theory [13][14] . The philosophy of such a novel theory is similar to Eringen's strain-driven nonlocal model, but by exchanging the positions of stress and strain terms, the integral constitutive relation is expressed as an integral convolution of the stress field instead of the strain field. Therefore, the constitutive BCs are converted from force-described ones to displacement-described ones, thereby eliminating the conflicts in BCs. So far, the stress-driven nonlocal model has been successfully used to investigate static and dynamic mechanical behaviors of beam structures, in which a consistently stiffening effect is captured [15][16] . However, the research objects of these studies focus on the classical or first-order shear beam models (i.e., the Euler-Bernoulli beam theory and the Timoshenko beam theory, respectively).
In fact, the Euler-Bernoulli beam theory is the most straightforward beam theory and is inaccurate for estimating the deformation of non-slender beams owing to the neglect of the shear-deformable effect. Although the Timoshenko beam theory can overcome the Euler-Bernoulli beam theory's limitations by taking into account the shear-deformable effect, it does not satisfy free shear stress on the upper and lower surfaces of the beam. As a result, a shear correction factor needs to be introduced to compensate for the difference between the actual stress state and the assumed one. To remove these defects and better predict beam behaviors, various types of higher-order shear deformation beam theory are developed, such as the third-order one, the sinusoidal one, the hyperbolic one, the exponential one, and the exponential one [17] . The shear strain distribution across the beam section is governed by a shear shape function, which meets the free traction conditions on surfaces, and thus no shear correction factor is required. Furthermore, Vo and Thai [18] presented a refined shear deformation theory, which divided the components of transverse displacement into bending components and shear components. In this model, the bending and shear forces are only determined by the bending component and the shear component of the beams, respectively. Also, the bending component is similar to the expressions in the Euler-Bernoulli beam theory, and the shear component is expressed as a specific high-order variation of shear strain [18] . Therefore, the shear stress will also vanish at the top and bottom surfaces. After that, various higher-order refined shear deformation beam theories (HORSDBTs) have been proposed [19][20][21] . Up to now, although several authors have studied the mechanical behavior of nanostructures based on the higher-order shear deformation theories [22][23][24] , most of the existing literature is still based on Eringen's ill-posed constitutive relation, that is, the inherent constitutive BCs are rarely considered.
Motivated by these reasons mentioned above, this work aims to extend those well-posedness nonlocal models to the application of higher-order refined shear deformation beams. The main novelty of this paper is a discovery that the stress-driven model is also ill-posed in higher-order refined shear deformation beams. As a remedy, the well-posedness and the consistency of both the strain-and stress-driven two-phase local/nonlocal mixed models are pertinently evidenced by studying the size-dependent bending of various high-order refine shear deformation beams. The content of this article is arranged as follows. In Section 2, the mathematical formulation of the bending problem of a functionally graded (FG) curved beam is established by using different versions of nonlocal elasticity. Such a beam model can degenerate into a homogeneous beam or a straight beam by adjusting related parameters. In Section 3, the ill-posedness of pure nonlocal elasticity of both the stress-and strain-driven types is pointed out, and the necessity to use the two-phase nonlocal models is illustrated. In Section 4, an effective method named as the generalized differential quadrature method (GDQM) is introduced to solve the coupling governing equations of curved beams. In Section 5, several numerical examples are considered to validate the consistency of the predicted scale-effects. Finally, the conclusions are given in Section 6.
2 Problem formulation 2.1 Higher-order refined shear deformation curved beam made of FG materials A slightly curved FG beam having the length L, the width b, the thickness h, and the constant curvature radius R together with the coordinate system is shown in Fig. 1. The properties of the beam are assumed to vary smoothly through the thickness of the beam by a power-law. Therefore, the useful material properties of the beam, including Young's modulus E, Poisson's ratio ν, and the shear modulus G are given by in which subscripts 1 and 2 denote the properties of materials 1 and 2, respectively. α is a non-negative power-law index, which governs the material distribution along the beam section. It is easy to know that α = 0 indicates homogenous cases (i.e., the full material 2), while α → ∞ indicates the full material 1.
According to the assumption of the HORSDBTs, the general forms of the circumferential displacement field u x and the radial displacement field u z of a curved beam can be expressed as follows [25] : in which u(x) is the tangential displacement of any point on the central axis, and w b (x) and w s (x) stand for the bending and shear components of radial displacement, respectively. Also, a shape function f (z) is taken into account to determine the shear strain and stress distributions across the beam's thickness. In this paper, the following various shape functions [26][27] are adopted.
(I) The third-order polynomial type of the HORSDBT (HORSDBT-T) (II) The sinusoidal type of the HORSDBT (HORSDBT-S) (III) The hyperbolic type of the HORSDBT (HORSDBT-H) (IV) The exponential type of the HORSDBT (HORSDBT-E) Moreover, a limit case f (z) = z corresponds to the Euler-Bernoulli beam theory. The non-zero strains of the curved beam can be given by [28]      where the prime symbol represents the differential with respect to the coordinate x.
It can be seen that the shape function form depends on the satisfaction of stress-free BCs on the top and bottom surfaces of the beam without utilizing a shear correction coefficient.
The variational principle of virtual work is used to derive the governing equations and BCs, which states where the variation of strain energy is On well-posedness of two-phase nonlocal integral models 935 and the virtual work done by the distributed external force q 0 and the concentrated force at beam ends is given by By inserting Eqs. (9) and (10) into Eq. (8) and setting the coefficients of δu, δw b , and δw s to zero, one can obtain the following governing equations: and the corresponding BCs at Those resultants in the above equations are defined as where A stands for the cross-sectional area of the beam.

Different versions of nonlocal constitutive relation 2.2.1 Strain-driven type
Eringen's strain-driven nonlocal model assumes that the nonlocal stress σ ij at any point is assumed to be a convolution integral of the strain ε ij at each point in the domain and a decaying kernel ψ. In this framework, the constitutive relation for a one-dimensional case is given by in which the so-called nonlocal parameter χ is used to depict long-range interactions in the region [x 1 , x 2 ]. However, some ill-posed issues appear when studying bounded structures. As a remedy, the strain-driven two-phase nonlocal (TPN-StrainD) model is proposed initially by Eringen, in which the stress-strain relation is expressed as a combination of classical local elasticity and nonlocal elasticity through a volume fraction ξ, In this paper, the classical Helmholtz kernel is adopted for all the nonlocal models mentioned, in which κ = χl e = e 0 a is a nonlocal length-scale parameter, and a and l e denote the internal and external characteristic lengths, respectively. e 0 is the non-negative nonlocal material constant. A combination of Eqs. (7), (13), (15), and (16) gives the strain-driven nonlocal constitutive equations of the FG curved beams, where the sectional rigidities are determined by Such a formula can revert to the strain-driven pure nonlocal and classical local theories by setting ξ = 1 and ξ = 0, respectively.

Stress-driven type
As another innovative strategy, a stress-driven nonlocal integral model was presented in Refs. [9], [13], and [14], which expresses the strain field at a certain point as a convolution integral of the stress field and the kernel at all points in the region, Following the successful application of the stress-driven nonlocal integral model to various size-dependent static and dynamic problems of simple beam models, Barretta et al. [29] also developed the stress-driven type of two-phase local/nonlocal mixed formulation, Such a two-phase strategy seems more applicable than the pure nonlocal model, because two control parameters can be used to reflect the nonlocal effect of different materials and structures.
Using the TPN-StressD strategy with the Helmholtz kernel (see Eq. (16)), the constitutive relations of the one-dimensional FG curved beams become On well-posedness of two-phase nonlocal integral models 937 This formulation corresponds to the stress-driven pure nonlocal theory and the classical local theory as ξ = 1 and ξ = 0, respectively.

Non-dimensional problem formulation
By introducing the following non-dimensional quantities: one can obtain the equilibrium equations in non-dimensional forms as (hereafter, dropping the asterisks for convenience) and the corresponding BCs at x = 0, 1 are Besides, the non-dimensional constitutive equations under the two different types of nonlocal models are, respectively, as follows.
(i) TPN-StrainD model According to Ref. [29], one can derive the following equivalent differential TPN-StrainD constitutive equations: with the following constitutive BCs:

938
Pei ZHANG and Hai QING (ii) TPN-StressD model According to Ref. [29], the TPN-StressD constitutive equations can be equivalent to with the following constitutive BCs: 3 Ill-posedness of pure nonlocal models

Strain-driven type
In the frame of Eringen's strain-driven pure nonlocal model, the static bending of higher-order refined shear deformation curved beams is governed by the equilibrium equation (23) and the constitutive law (25) with ξ = 1. At the same time, the natural BCs (24) and the constitutive BCs (26)-(29) with ξ = 1 are mandatory. Now, assuming ξ = 1 and substituting first three equations of Eq. (25) into Eq. (23), one can easily find that the paradox existing in classical simple beam theories is still unavoidable; that is, the problem has more BCs than those required (i.e., 18 mandatory BCs and 12 required BCs). Since all BCs cannot be satisfied simultaneously, the problem may be over-constrained and usually cannot be solved. Besides, the influence of nonlocal parameter on the deformation of the cantilever beams under end-point load is inconsistent with that of other supported boundaries, which is mainly due to the conflict between the natural BCs (see the second equation to the fifth one of Eq. (24)) and the constitutive BCs (see Eqs. (27)-(28)).
In the frame of Euler-Bernoulli and Timoshenko beam theories, such paradoxes of Eringen's strain-driven pure nonlocal model can be entirely overcome by using the stress-driven model. By exchanging stress and strain terms in the integral constitutive relation, the inherent constitutive BCs are transformed from stress-described ones into strain-described ones, thereby eliminating the conflict between the natural BCs and the constitutive BCs. Also, swapping the positions of the strain field and the stress field in the constitutive equation increases the number of required BCs, and thus the problem becomes solvable [16,30] . However, for higher-order shear deformation beam models, this is not the case.

Stress-driven type
In the frame of the stress-driven pure nonlocal model, the equilibrium equation and constitutive laws for the problem can be obtained from Eqs. (23) and (30)  All in all, the two commonly used pure nonlocal elasticity models are both ill-posed when analyzing nonlocal higher-order refined beam models. Therefore, it is necessary and mandatory to try another strategy, i.e., two-phase nonlocal formulation.
4 Solution procedure for two-phase nonlocal models 4

.1 Equations in terms of displacements
Since the derived equations with two-phase nonlocal models are lengthy and complex, it is hard to obtain their analytical solutions accurately. As an alternative, a numerical solution method, named as the GDQM, is adopted in this paper. Before this process, the governing equations and the corresponding BCs should be expressed in terms of displacements u, w b , w s and variable Q.

With TPN-StrainD model
Using the TPN-StrainD model, the governing equations can be re-expressed in terms of the displacement components, g 11 u + g 12 u + g 13 w b + g 14 w b + g 15 w b + g 16 w s + g 17 w s + g 18 w s = 0, (35) g 21 u + g 22 u + g 23 w b + g 24 w b + g 25 w b + g 26 w b + g 27 w s + g 28 w s + g 29 w s + g 210 w s − R 2 q 0 /(R 2 + κ 2 ) = 0, (36) g 31 u + g 32 u + g 33 w b + g 34 w b + g 35 w b + g 36 w b + g 37 w s + g 38 w s and the constitutive BCs become The natural BCs for different supported types can also be expressed in terms of displacements, which are not listed here for simplicity. Moreover, the explicit expressions for the coefficients n ij , m ij , g ij , and c ij are not listed here owing to the limited space of the article.

With TPN-StressD model
For the TPN-StressD model, the governing equations can also be expressed in terms of displacements, namely, and the constitutive BCs become On well-posedness of two-phase nonlocal integral models Similarly, the natural BCs in displacement forms are not listed here for brevity. Moreover, the explicit expressions for the coefficients N ij , M ij , G ij , and C ij are not listed here owing to the limited space of the article.
Remark 1 According to the above formulation, it can be found that, when using the TPN-StrainD or TPN-StressD model, the orders of the unknown variables, i.e., u, w b , w s , and Q, are 4, 6, 6, and 2, respectively, which means that 18 BCs are required. This number is exactly equal to the number of mandatory BCs for the problem (i.e., 10 natural BCs and 8 constitutive BCs). That is to say, the ill-posedness in the pure nonlocal models is utterly avoided by using two-phase nonlocal strategies. Therefore, when solving the size-dependent problems of the higher-order refined beam models, a two-phase nonlocal strategy must be used to replace the pure nonlocal model, regardless of the strain-and stress-driven nonlocal types.

Solution procedure
By using the GDQM [31] , the coupled differential equations with arbitrary BCs can be solved. First, for faster convergence and better calculation accuracy, the grid points are considered as Thereafter, the following vectors are defined: , · · · , w bn , w b1 , w bn , w b1 , w bn ) , d ws = (w s1 , w s2 , · · · , w sn , w s1 , w sn , w s1 , w sn ) , while the approximation forms of the unknown variables can be expressed by using different interpolation functions, where n is the number of the grid points, ϕ u/w (x), ψ u/w (x), and ϑ w (x) are the Hermite interpolation functions [31][32] , and l j (x) denotes the Lagrange interpolation function. Based on these assumptions, the governing equations and BCs of the problem can be expressed in matrix forms as below, through which the unknown variables can be carried out.

Validation
Since the ill-posedness of the pure nonlocal models exists, and there is no relevant research on higher-order refined beams related to the two-phase nonlocal models, the efficiency and accuracy of the present solution procedure are validated by comparing the results of a degradation case (i.e., ξ = 0.000 1, κ = 0.001) with those of the local elasticity theory. In calculations, the number of grid points is set as n = 21. For comparison, the FG parameter is set as α = 0, and the material properties of the beam are assumed to be the same as in Ref. [26], while the following non-dimensional deflections are defined [26] : From Table 1, one can see that, for different values of the length to height ratio L/h, as the curvature radius R = 10 000 (i.e., approximately regarded as infinite), the difference between the present solutions and the results of those straight beams in Ref. [26] can be negligible, regardless of TPN-StrainD and TPN-StressD models. Table 1 Comparison of non-dimensional midpoint deflections w b (0.5) + ws(0.5) of simply-supported (SS) curved beams (R = 10 000) subjected to uniformly distributed loads, predicted by TPN models (with ξ = 0.000 1, κ = 0.001) and local theory Beam theory TPN-StrainD TPN-StressD Local theory [26] L/h = 5 L/h = 20 Moreover, a comparison between a limit case (f (z) = z, g(z) = 0) of present models and the two-phase nonlocal Euler-Bernoulli beam model [29] is made in Table 2. Also, good agreement is obtained for both the TPN-StrainD and TPN-StressD models. These comparative results show that the present solution procedure is effective and accurate.

Consistency of nonlocal effects
In this subsection, the numerical examples are considered to study the consistency of current two-phase nonlocal integral models for addressing the size-effect in bending behavior of various higher-order refined shear deformation curved beams (i.e., based on all of the shear functions HORSDBT-T, HORSDBT-S, HORSDBT-H, and HORSDBT-E) made of FG materials. Different boundary edges and loading conditions are considered, including SS, clamped-clamped (CC), and CF beams subjected to a uniformly distributed load q 0 , as well as CF beams under an end-point load Q. Moreover, it is worth mentioning that, unless otherwise stated, the definitions (see Eqs. (55) and (56)) in Subsection 5.1 will continue to be used in subsequent sections. In Figs. 2-5, the subfigure (a) depicts the variation of the dimensionless midpoint or tip deflections versus the volume fraction parameter ξ with a fixed dimensionless nonlocal parameter κ = 0.15, and the subfigure (b) presents the dimensionless midpoint or tip deflections versus the dimensionless nonlocal parameter κ with a fixed volume fraction parameter ξ = 0.3. In calculations, the FG power-law index, the dimensionless curvature radius, and the length to height ratio of the beam are assumed to be α = 0, R = 10, and L/h = 10. It can be seen from these pictures that, for all boundary types, the TPN-StainD nonlocal integral model shows a consistent softening impact on beam deformations. Increasing the value of the dimensionless length-scale nonlocal parameter κ or the volume fraction of nonlocal part ξ leads to an increase in dimensionless deflections of the beam. On the contrary, there is a stiffening effect on the beams for the stress-driven strategy as the nonlocal length-scale parameter increases. Moreover, the larger the volume fraction of the nonlocal stress part, the more significant the stiffening effect of the nonlocal length-scale parameter on beams. That is to say, both types of two-phase nonlocal models can predict consistent size-dependent responses. Therefore, the two-phase nonlocal integral strategies are suitable for addressing the is considered in Fig. 6. One can find that the dimensionless deflections of SS and CC beams decrease significantly with an increased length to height ratio. When the slenderness ratio is greater than 10, the deflection change of the CF beams is relatively tiny. Figure 7 illustrates the variation of the non-dimensional midpoint or tip deflections versus the curvature radius of different supported curved beams based on both TPN-StrainD and TPN-StressD models with different volume fractions ξ = 0.25, 0.5, and 0.75. The HORSDBT-S is utilized. In contrast to the stress-driven pure nonlocal model [30] , increasing the curvature radius leads to the increased beam deflections for all boundary types. Moreover, with the increase in the curvature radius, the nonlocal effects on SS and CC curved beams are enhanced by increasing the volume fraction ξ. However, for cantilever beams, the effect of nonlocal volume fraction seems to have nothing to do with the radius of curvature of the curved beams, because the shapes of the deflection curves are almost the same at different values of the radius of curvature.  volume fractions ξ = 0.25, 0.5, and 0.75. It can be seen that the non-dimensional deflections of all beams consistently increase with the increase in the FG parameter. As the FG index α increases, the proportion of the material 1 increases, making the sectional stiffness of the beams decrease. As the index α continues to increase, the FG influence gradually stabilizes.

Conclusions
In this paper, the ill-posedness of both the strain-and stress-driven models for higher-order refined beam models is uncovered. As a remedy, the well-posedness of the two-phase strategy of nonlocal integral models is illustrated by analyzing the scale-effected bending of a curved beam made of FG materials. The governing equations and the corresponding BCs are established by invoking the variational principle of virtual work. The two-phase nonlocal integral constitutive law is converted to an equivalent differential one with two constraint BCs. With the GDQM, the coupling governing equations are solved numerically. Then, several numerical examples are given for investigating the consistency of the two-phase nonlocal models. From the above results, the conclusions are given below.
(i) Unlike the case based on the classical simple beam theories, the strain-and stress-driven pure nonlocal models are both ill-posed for higher-order refined shear deformation beam models. There are too many BCs in such problems to be met simultaneously, thereby leading to ill-posed issues.
(ii) With the utilization of the two-phase nonlcoal strategy, the number of mandatory BCs for the problem happens to equal the number of BCs required, and the consistent size-effects can be predicted for different boundary and loading conditions. Therefore, it is necessary and mandatory to adopt the two-phase local/nonlocal formulation when analyzing the higher-order shear deformation beam models, regardless of the strain-and stress-driven types.
(iii) Furthermore, the effects of the length to height ratio, the curvature radius, and FG materials are investigated extensively. It is found that the dimensionless deflections decrease with the increase in the length to height ratio. Increasing the curvature radius leads to an increase in beam deflections. With the increase in the curvature radius, the nonlocal effect of the volume fraction ξ on SS and CC curved beams is significantly enhanced. The non-dimensional deflections of curved beams under all supported conditions consistently increase with the increase in the FG parameter, and as the index α continues to increase, the FG effect on the beam gradually stabilizes.
This study reveals the ill-posedness of both the strain-and stress-driven pure nonlocal models for various higher-order refined beam models. The conclusions and the present solution methods are also applicable to other structural problems related to high-order shear deformation assumptions.