Superposed Incremental Deformations of an Elastic Solid Reinforced with Fibers Resistant to Extension and Flexure

A comprehensive linear model for an elastic solid reinforced with fibers resistant to extension and flexure is presented. ,is includes the analysis of both unidirectional and bidirectional fiber-reinforced composites subjected to in-plane deformations. Within the prescription of the superposed incremental deformations, the fiber kinematics are approximated and used to determine the Euler equilibrium equations. ,e constraints of bulk incompressibility and admissible boundary conditions are also discussed for completeness. In particular, the complete systems of differential equations are obtained for the cases of NeoHookean and Mooney–Rivlin types of materials from which analytical solutions can be obtained.


Introduction
e mechanics of microstructured solids have consistently been the subject of intense research [1][2][3][4][5] for their practical importance in materials science and engineering.In particular, a considerable amount of attention is committed to the development of continuum models and analyses in an effort to predict the mechanical responses of fiber-matrix systems subjected to external forces and/or induced deformations (see, for example, [6,7] and references therein).Continuum-based approaches postulate continuous distribution of fibers within the matrix materials so as to establish the idealized description of homogenized fiber-matrix composites. is is framed in the setting of anisotropic elasticity where the response function depends on the first gradient of deformations, typically augmented by the constraints of bulk incompressibility and/or fiber inextensibility.e latter condition often results highly constrained prediction models so that the corresponding deformation fields are essentially kinematically determinate, especially that arises in fibers [6,7].e approach has clear advantages in the prediction of deformation profiles of the system via the deformation mapping of an individual fiber, yet rather insufficient in the estimation of overall material properties of the fiber-matrix system.Nonetheless, the aforementioned models have been widely adopted in the analysis of composite materials for their merit in the continuum description and the associated mathematical framework [5][6][7].For the estimation of resultant properties, one may also consider multiscale modeling method which integrates the material properties of continua assigned in different length scales by means of the Cauchy-Born rule.Examples of such practice can be found in the works of Shahabodini et al. [8,9].
A considerable advance in the continuum theory of fiber-reinforced solids was made in recent years. is includes the incorporation of the bending resistance of fibers into the models of deformations where elastic resistance is assigned to changes in curvature of the fibers [10][11][12].More precisely, the fibers are regarded as convected curves so that the bending deformation of fibers can be formulated via the second gradient of deformations explicitly [13].e concept has been successfully adopted in a wide range of problems arising in materials science [14][15][16][17], and the mathematical perspective of the subject is discussed in [18][19][20].e authors in [21] proposed a general theory for the mechanics of an elastic solids with fibers resistant to flexure, stretch, and twist within the simplified setting of the Cosserat theory of nonlinear elasticity [10,22].
Further, the second gradient theory of elasticity for the mechanics of meshed structures is presented in [23][24][25].To this end, the authors in [26][27][28][29] developed continuumbased models in the analysis of fiber-reinforced composites, where the extension and bending resistance of fibers are incorporated via the computations of the first and second gradient of deformations.e majority of the aforementioned studies address nonlinear continuum theory of fiber composites so that little has been devoted to the development of compatible linear models describing the mechanics of an elastic medium reinforced with fibers.
In the present work, we present comprehensive linear theory of the strain gradient elasticity for the mechanics of fiber-reinforced composites with fibers resistant to bending and extension.e kinematics of fibers is approximated with the prescription of the superposed incremental deformations.We formulate, in cases of Neo-Hookean types of materials, complete expressions of the linearized Piola stresses from which the Euler equilibrium equations and the admissible boundary conditions are obtained.Bidirectional fiber composites are accommodated via the decomposition of the deformation gradient tensor along the directions of fibers.In addition, the reduction of simultaneous differential equations into a single differential equation is demonstrated by utilizing the compatibility conditions of the response function.More importantly, we show that, even with the introduction of the second gradient of deformations, two different bases (referential and current coordinate) do indeed merge for "small" deformations superposed in large.Lastly, Mooney-Rivlin types of materials are elaborated where we show that the resulting Euler equilibrium equations are of the same form as those obtained from the Neo-Hookean model.e corresponding Piola stresses, on the other hand, are distinguished in that they demonstrate clear dependency on both the first and second invariant of the deformations gradient tensor.e obtained model can be easily adopted in the study of composite structures subjected to small in-plane deformations.For example, in the design stage of crystalline nanocelluloses (CNCs) composite, the overall responses and deformation profiles can be predetermined by utilizing the proposed model.In addition, the underlined theory can also be extended to the deformation analysis of lipid membranes (e.g., budding and thickness distension) [30][31][32], where phospholipids (microstructure) are aligned to the normal direction of the membrane, and therefore, their deformations can be mapped and computed using the proposed model.
roughout the paper, we make use of a number of wellestablished symbols and conventions such as A T , A −1 , A * , and tr(A).ese are the transpose, the inverse, the cofactor, and the trance of a tensor A, respectively.e tensor product of vectors is indicated by interposing the symbol ⊗ , and the Euclidian inner product of tensors A and B is defined by . e symbol | * | is also used to denote the usual Euclidian norm of three vectors.Latin and Greek indices take values in {1, 2} and, when repeated, are summed over their ranges.Lastly, the notation F A stands for the tensor-valued derivatives of a scalar-valued function F(A).

Incremental Elastic Deformations and Equilibrium Equations
e incremental deformation is defined by (see, for example, [33] and the references therein) the following equation: where ( * ) o denotes configurations of * evaluated at ε � 0 and ( _ * ) � z( * )/zε.In the forthcoming derivations, we define _ χ � (zχ/zε) � u for the sake of convenience and clarity.From equation (1), the gradient of the deformation function χ(X) can be approximated up to the leading order as shown in the following equation: where _ F � (z _ χ/zX) � ∇u.In the above equation, we assume that the body is initially undeformed and stress free at ε � 0, that is, F o � I and P o � 0, where P is the Piola stress.us, equation (2) becomes the following equation: and successively yields which can be found by using the identity FF −1 � I. e determinant of F can be approximated similarly as where we evaluate � tr(gradu) � divu, we obtain from equation (5) that In addition, the Euler equilibrium equation can be expanded as the following equation: Dividing the above equation by ε and letting ε ⟶ 0, we obtain which serves as the linearized Euler equilibrium equation.e expression of Piola stresses in the case of fiber composites reinforced with fibers resistant to flexure is given by [27]: where W, p, and C are the energy density function, Lagrange multiplier, and material constant of fibers (C � constant), respectively.Also, g is the geodesic curvature of fibers' trajectory (i.e., g � G(D ⊗ D)) and D is the initial director filed of fibers' where D � (zX(S))/zS.Here, S is the arc length parameter on the reference configuration.In general, most of fibers are straight prior to deformations.Even slightly curved fibers can be regarded as "fairly straight" 2 Advances in Materials Science and Engineering fibers, considering their length scales with respect to that of matrix materials.us, from here and forthcoming derivations, it is assumed that zD zX � 0, Accordingly, we find Div(D) � 0 and thereby reduce equation (10) to the following equation: where the last term of the above equation is obtained by using the identity us, from equation (11), _ P can be evaluated as the following equation: In the case of an incompressible Neo-Hookean material, the energy density function is given by the following equation: and we find In view of equations ( 8) and ( 14), the linearized Euler equilibrium equation then satisfies e evaluation of equation ( 15) is essential to extract boundary value problems (BVPs); nonetheless, the details are often heavily omitted in the literature (see, for example, [26][27][28][29]).To see this, we first compute the following equation: In the above equation, caution needs to be taken, in which ∇( * ) and div( * ) are the operators in the reference frame.Although, there are no clear distinction between the reference and current configurations for "small" incremental deformations, the mathematical procedure should reasonably address their connections especially when dealing with tensors with mixed bases (e.g., ∇u, F).
e collapse of bases for the present problem will be discussed in the later sections.e second term in equation (15) becomes where we use the Piola identity (i.e., F * iA,A � 0) and However, in order to recover the initial stress free state at ε � 0, we require from equations ( 11) and ( 13) that and thus find p o � μ � constant.erefore, equation (18) becomes and Div(F * ) � 0 (Piola's identity).Also, Div(∇ _ g(D ⊗ D)) can be evaluated as where g i � G iAB D A D B (see [27]) and is the second gradient of deformations (i.e., ∇F � G).Consequently, by substituting equations ( 16), ( 17), (20), and (21) into equation ( 15), the linearized Euler equation can be obtained as follows: For a single family of fibers (i.e., D � E 1 , D 1 � 1, D 2 � 0), equation (22) reduces to

Boundary Conditions and Solution to the Linearized System
e corresponding boundary conditions are given by [27] t where t, m, and f are, respectively, the expressions of edge tractions, edge moments, and the corner forces.Further, N and T are unit normal and tangent to the boundary.
e "small" increment of boundary forces are then computed as follows: e expression of _ P can be obtained from equation ( 14) that Advances in Materials Science and Engineering In order to apply boundary tractions (e.g., _ P 11 ), it is necessary to compute equation ( 26) as a function of u � u(X 1 , X 2 ).For the purpose, we first find the following equation: where iA (e i ⊗ E A ), we use the chain rule in the form of the following equation: Here, the expression of (zF * iA /zF jB ) can be found via the connection [34]: us, at ε � 0, we have from the above equation that and thereby obtain Consequently, substituting equations ( 27), (28), and (32) into equation ( 26) furnishes from which the expression of boundary tractions is completely determined in terms of u.
Remark 1.In the above equation, δ iA δ jB u j,B is interpreted as δ iA u B,B resulting δ iA δ jB u j,B � δ iA Divu in the reference configuration (Eulerian).However, one may also find δ iA δ jB u j,B � u j,j and thus obtain δ iA δ jB u j,B � δ iA divu in the current configuration (Lagrangian).is confirms the wellknown result from the linear elasticity theory that the bases collapse in the event of small deformations superposed on large (i.e., Divu � divu and E A � e i ).In the present problem, the result allows one to formulate linearized Euler equations without conflicting bases mismatch especially when operating linear transform of mixed basis tensors.
Although there are no clear distinctions between the current and deformed configurations, caution needs to be taken that the Euler equation (Div _ P) is, in principal, defined in the reference frame together with the boundary conditions.Now, in view of equations ( 5) and ( 6), the constraints of bulk incompressibility reduce to the following equation: Equation ( 34) serves as the linearized bulk incompressibility condition (i.e., u i,i � 0), which needs to be solved together with equation (23).In addition, since Divu � Divu for small deformations (see Remark 1), we find Divu � Divu � 0 and thereby reduce equation (33) to the following equation: which serves as an explicit form of the linearized Piola stress for elastic solids reinforced with single family of fibers.In particular, if the fibers' directions are either normal or tangential to the boundary (i.e., and therefore, we compute, for example, where N is assumed to be parallel to the fiber's director field (i.e., N � D � E 1 ).Lastly, from equations ( 27) and (36), the expression of boundary moments are obtained by 3.1.Solutions to the Linearized System.In the case of single family of fibers (i.e., D � E 1 ), the linearized system of partial differential equations (PDEs) is given by e second one in the above equation can be automatically satisfied by introducing the following scalar field φ: An analytical solution of the above exists (see [27]) and is completely determined by imposing admissible boundary conditions as discussed in equation (25).For example, the symmetric bending can be imposed as which serves as the boundary conditions for the equation (42).

Extensible Fibers
e Piola stress in the case of initially straight and extensible fibers is given by Zeidi and Kim [28]: where E is the elastic modulus of fibers (extension).Further, the fibers' stretch λ can be computed as In the case of Neo-Hookean type materials (see equation ( 13)), equation (45) becomes us, equations ( 1) and (47) can be approximated as Since F o � F * o � I, the above equation further reduces to the following equation: where To obtain the desired expression, it is required to compute the following equation: In the above equation, the initial director field (D) is represented by the current frame (i.e., e i ).However, D is, in principal, a vector in the reference coordinate (i.e., D � D A E A ). is is due to the collapse of bases as discussed in Remark 1. Caution needs to be taken when applying the Einstein summation.us, we obtain from equations ( 32), ( 34), (49), and (50) that where 51) can be used as the expression of boundary tractions in the case of extensible fibers.For example, if D � E 1 (single family of fibers), the above equation yields the following equation: Further, from equations ( 8) and ( 49), the corresponding equilibrium equation then satisfies the following equation: e evaluation of the above equation is well discussed through equations ( 16)-( 21) except the stretch term which is now equated as (54) erefore, we obtain the following equation: e linearized boundary conditions in the case of extensible fibers remain intact (see [28]), except the expression of _ P where the explicit expression is obtained in equation (52).us, for example, in the case of unidirectional fiber composited where a boundary vector N is parallel to the fiber's director field (i.e., N � D � E 1 ), the complete set of equations can be found as Also, the corresponding boundary conditions are given by Advances in Materials Science and Engineering where from equation (52), _ . Lastly, we note here that the value of p o , in the case of extensible fibers, can be obtained by evaluating the corresponding stresses (see equation ( 47)) at ε � 0 which is again found to be p o � μ.By applying the similar scheme as applied in Section 3.1, we rearrange equation ( 56) and thereby obtain e solution of equation ( 58) is not accommodated by conventional methods such as the separation of variables method, Fourier transform, and polynomial solutions.Instead, a particular form of solution (φ � X(x)sin(my)) can be proposed, inspired by the modified Helmholtz equation.

Extensible Bidirectional Fibers
In the forthcoming derivations, we confine our attention to the case of initially orthogonal fiber families.
e bidirectional fibers can be accommodated by using the following decompositions of the deformation gradient: where L and M are the director fields of each fiber family in the reference configuration and l and m are their counterparts in the current configuration.e fibers' stretches are computed similarly as in equation ( 46) that Equations ( 59) and (60) complete the deformation mapping, for example, L � L A E A and l � l i e i to yield Further, by employing the variational principles, the Piola stress for the bidirectional and extensible fiber composites can be found as where E i and C i are the material constants of fiber families, respectively, to extension and flexure.Applying the similar scheme as in Section 4, we approximate equation (62) and successively obtain erefore, equations ( 8) and (63) furnish where we evaluate Compatible results can be also found directly from equations ( 17) and (20) (i.e., Div( _ pF * o ) � _ p i e i and Div(p o _ F * ) � 0).In the case of bidirectional and extensible fibers, the corresponding boundary conditions are given by (see [26]) the following equation: (66) Now, we approximate the above boundary conditions and thus obtain If the fibers' directions are aligned with the axes of Cartesian coordinates (i.e., L � E 1 and M � E 2 ) and are either normal or tangential to the boundary (i.e., 64) and (67) further reduces to the following equation: Advances in Materials Science and Engineering which together with the incompressibility condition (u i,i 0) constitutes the complete set of the PDEs system to solve the nal deformation pro les.Similarly, equation (63) now becomes For example, if the unit normal of a boundary is N E 1 , the boundary traction can be computed as follows: Again, by introducing the scalar eld φ (see equation ( 40)) and employing the compatibility conditions (i.e., _ p ij _ p ji ), we reduce the above equation to the following equation: A complete analytical solution for equation ( 71) is available via the methods of iterative reduction and the principle of eigenfunction expansion (see, more details, [35][36][37]).Figures 1-4 illustrate the deformation pro les for the cases discussed through Sections 3-5.As the equations become mathematically compact, the corresponding deelds experience less oscillatory behaviors.is clearly demonstrated by results in Figures bidirectional cases, show a close agreement the models.We also mention here that equation (71) reduces to equations (42) and (58) in the limit of vanishing C 2 and E 2 , respectively.is, in turn, suggests that unidirectional cases can be assimilated within the systems of equation ( 71) by setting C 2 E 2 0 (see Figure 5).In fact, the latter model is recommended to use, since the bidirectional and extension bers case is the most general and compact form (with minimal singular behaviors) and therefore produces more accurate prediction results.
To elaborate the proposed model, we present comparisons with the experimental data obtained from the 3-point bending test of CNC ber composite (C 150 GPa, μ 1 GPa). is is a particular case of the proposed model, when C/μ 150 and vanishing E (i.e., equation (20)).Using the method presented in Section 3.1, we nd the solution of equation (20) as where Using equation (74), we compute the maximum deections of the CNC composite and illustrate the results in Figure 6.It is shown in Figure 6 that the proposed model predicts the normal de ections of the CNC composite strip.Since the slope of the graph in Figure 6 indicates the moduli of the composite (when divided by the corresponding length scale), it can also be used in the estimation of the overall material properties of the posite.due to the paucity of available data, we are not able to provide the quantitative analysis other than those presented in Figure 6 at this time.
e study is currently underway, and our intention is to present elsewhere when ready.We also mention here that the overall responses of composite materials can also be estimated using multiscale modeling methods (see, for example, Shahabodini et al. [8,9]).
Lastly, we assimilate the case of single family of ber composites subjected to the axial tension.e linear solution obtained from equation (56) demonstrates good agreement with the compatible nonlinear model [28] for small deformations superposed on large, while it shows discrepancies in the prediction of large deformations (See Figure 7).It is also noted here that we reserve the details in solving the corresponding di erential equations (equation ( 56)) for the sake of conciseness.However, the procedures for the particular case of the present example can be found in [28].

Further Considerations
For the analysis of soft materials-based composites, such as carbon rubber-ber composites and polymer composites, a di erent type of energy potential may be suggested instead of the Neo-Hookean model discussed in the previous sections.e Mooney-Rivlin model is one of the most commonly used energy potentials for the aforementioned cases (see, for example, [33]).In the forgoing development, we present a compatible linear model for the Mooney-Rivlin types of materials for the desired applications.
e expression of the Mooney-Rivlin potential is given by the following equation: e derivative of equation ( 75) with respect to the deformation gradient tensor then yields the following equation: Substituting equation (76) into equation ( 9) furnishes the following expression of the Piola stresses: (77) Now, by applying the same schemes as adopted in the previous sections, we nd the following equation: where the second term of the right side of equation (78) can be evaluated at ε 0 as us, from the results in equations ( 2), ( 27), ( 28), ( 33), (78), and (80), we nd the following stress expression for the cases of unidirectional ber composites as where _ F BB u B,B 0 from the conditions of bulk incompressibility (see equation (34)).Consequently, the corresponding Euler equilibrium equation satis es the following equation: Using the compatibility condition of u A,iA (i.e.u A,iA u A,Ai ) and the constraints of bulk incompressibility, we nd that u A,Ai (u A,A ) i 0 and thereby reduce equation (83) to the following equation: e resulting equilibrium equation (equation ( 83)) implies that there is no clear distinction between the Neo-Hookean model and Mooney-Rivlin model as far as the superposed incremental deformations with the augmented bulk incompressibility conditions are considered.is is due to the fact that the divergence of high-strain terms identically vanishes upon the imposition of linearized bulk incompressibility conditions (i.e., Divu divu 0).However, the expressions of Piola stresses are a ected by the introduction of Mooney-Rivlin type materials (see equations (33) and ( 81)).us, for example, if N D E 1 , we compute from equation ( 25) that which is the boundary tractions for the case of the Mooney-Rivlin energy potential.Lastly, from equation (77), we nd at ε 0 that and thereby obtain μ p o .Equation (83) together with the imposed boundary conditions (i.e., equations ( 25) and ( 84))  constitutes a complete set of PDEs system which can be solved using the similar schemes as presented in Section 3.1.Further, because the resulting equilibrium equations are of the same form (see equations (22) and (84)), equations ( 81) and (84) may be used in the determination of the parameters associated with high-strain terms (λ) by comparing stresses in equations ( 33) and (81).Investigations in this respect (including the implementations of the developed linear theory) are currently underway.Our intention is to report elsewhere.

Conclusions
We present complete linear models for the mechanics of an elastic solid reinforced with fibers resistant to flexure and extension.Within the prescription of the superposed incremental deformations, the first and second gradient of deformations is approximated through which the kinematics of fibers is explicitly determined.
e linearized Euler equilibrium equations and the Piola stress are then formulated for the Neo-Hookean types of materials.We also derive the corresponding boundary conditions and the conditions of bulk incompressibility for the sake of completeness.In addition, the cases of the bidirectional fiber composites are considered via the fiber decomposition of the deformation gradient tensor.It is found that the systems of PDEs for the bidirectional and extensible fibers reduce to those from single family of fibers in the limit of vanishing material parameters of fibers.
In particular, we demonstrate the well-known result from the linear theory of elasticity that the merging of the bases remains valid even with the incorporation of fibers' bending and extension into the models of deformations.A Mooney-Rivlin material is also elaborated, where we show that there is no clear distinctions in the resulting Euler equations obtained from the two different strain energy models (i.e., Neo-Hookean and Mooney-Rivlin) as long as the small deformations are concerned.However, the resulting Piola stresses are of different forms in that the one from the Mooney-Rivlin model demonstrates clear dependency on both the first and second invariant of deformation gradient tensors, whereas the Piola stress, in the case of the Neo-Hookean model, depends only on the first invariant.Lastly, we mention here that the present model can be used in the approximation of the nonlinear theory of strain gradient elasticity arising in finite plane elastostatics.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Figure 5 :Figure 4 :
Figure 5: Reduction of solutions in the limit of C 2 E 2 0.