Elastica of Non-Prismatic and Nonlinear Elastic Cantilever Beams under Combined Loading

: This study presents the elastica of non-prismatic cantilever beams with rectangular cross-sections that are subjected to combined loading. The considered beams are nonlinearly elastic and obey Ludwick’s constitutive law. The combined loading system used in this study provides uniform loading, tip point loading, and tip couple loading individually or in combination. This loading system can create a total of seven loading cases that have not been covered in the literature. Ordinary differential equations governing the large deformed shapes of the elastica of the beams are derived and solved numerically. The effects of beam parameters on elastica behavior, including tip responses and strains and stresses loaded onto the cross-sections, were studied. It was observed that the exponential constant of the mechanical properties is very sensitive to elastica behavior. The results also demonstrate that the stresses near the neutral axis are larger with a larger exponential constant.


Introduction
Elastica problems are frequently encountered in many areas of modern structural engineering. The first study on elastica was published by Euler [1], who studied the large deformed shapes of a slender rod. The elastica problem of geometric and material nonlinearity has been discussed by many researchers in recent years. Additionally, non-prismatic beams are commonly used in engineering practice for their economic, aesthetic, and optimization benefits. A number of studies have been published on these topics. For example, the elastica behavior of non-prismatic cantilever beams fabricated from nonlinear elastic materials, such as Ludwick-type materials, under combined loading has been studied over the past few decades.
Representative works in the same area as this study are listed and briefly discussed below. Oden and Childs [2] studied the elastica of buckled bars represented by moment-curvature relationships based on the hyperbolic tangent law. Prataph and Varadan [3] investigated elastica under tip point loading, where the beam material followed the Ramberg-Osgood stress-strain law. Lewis and Monasa [4,5] computed the elastica of prismatic beams subjected to point loading and tip couple loading, respectively, for beams obeying Ludwick's constitutive law. Varadan and Joseph [6] computed the elastica of prismatic beams fabricated from Ramberg-Osgood-type materials under tip couple loading. Fertis and Lee [7] conducted considerable research on the elastica of prismatic and non-prismatic flexible beams using the method of equivalent systems. Lee [8] studied the elastica of prismatic beams obeying Ludwick's constitutive law that were subjected to combined loading consisting of uniform loading and tip point loading. Eren [9] computed the elastica of rectangular cantilever beams fabricated from Ludwick-type materials under combined loading based on different arc length assumptions. Brojan et al. [10] investigated the elastica of non-prismatic beams under tip couple loading that obeyed the generalized Ludwick's constitutive law. Borboni et al. [11] studied the large deflection of a nonlinear, elastic, asymmetric, Ludwick cantilever beam subjected to point loading at the free end. Liu et al. [12] investigated the large deflection of curved elastic beams fabricated from Ludwick-type materials subjected to uniform loading and concentrated vertical loading at the free end. Changizi et al. [13] studied a close-form solution for the nonlinear deflection of non-straight, Ludwick-type beams subjected to point loading at the free end using Lie symmetry groups. Lee and Lee [14] computed the elastica of non-prismatic cantilever columns fabricated from Ludwick-type materials that were subjected to axial point loading.
In all the works discussed above, uniform loading, tip point loading, and tip couple loading were not examined simultaneously in a combined loading system. The elastica problem considered in the present paper does not absolutely satisfy the superposition principle. Therefore, it is important to include as many individual loading cases as possible in a combined loading system because the net response to multiple stimuli is not equal to the sum of the responses to individual stimuli.
This study aimed to examine the elastica of a cantilever beam subjected to a combination of three kinds of loads, namely uniform loading, tip point loading, and tip couple loading, simultaneously. The loading system used in this study can generate a total of seven loading cases that have not been covered in the literature. The cross-section of the non-prismatic beam is solid and rectangular with a width and height that vary linearly along the beam axis. The nonlinear beam material obeys Ludwick's constitutive law.
For the analyses of elastica in this study, the following assumptions were made. The neutral axis is inextensible, the effects of Poisson's ratio and transverse shear deformation are negligible, and the strain caused by bending is small. Additionally, the tip point and uniform loads were assumed to be vertical after deformation.
We will begin by formulating the linearly varying width and height of the rectangular beam. The differential equations governing the elastica of the deflected beam are then derived based on the Bernoulli-Euler beam theory and solved numerically using the Runge-Kutta and Regula-Falsi methods. The tip responses in this study are compared to those in the literature for validation purposes. Finally, for beams fabricated from steel, the NP8 aluminum alloy, the annealed copper, and the results of elastica behavior, including tip responses and strains and stresses, are discussed in detail.

Geometry of the Non-Prismatic Beam
In general, the configuration of a non-prismatic member is arbitrary. In this study, the geometry of a cantilever beam with a span length of l ( Figure 1) was selected as a representative non-prismatic member. The beam cross-section is solid and rectangular with a width and height that vary linearly in the z and y directions, respectively, along the axial length s. The width and height along s are denoted as B and H, respectively. Those at the clamped end c (s = 0) are denoted as B c and H c , respectively. Those at the free tip f (s = l) are denoted as B f and H f , respectively. It should be noted that the z axis is the bending axis.
To express B and H along s functionally, the two following non-dimensional section ratios are defined: Considering r B and r H as defined above, B and H are expressed as linear functions of s as follows:

of 13
The geometry of other types of non-prismatic beams other than the aforementioned linearly varying rectangular cross-section can be easily formed in a similar manner.
The geometry of other types of non-prismatic beams other than the aforementioned linearly varying rectangular cross-section can be easily formed in a similar manner.

Mathematical Modelling
The symbols and loads for the cantilever beam are defined in Figure 2. The beam was subjected to three types of loading, either individually or in combination: a uniform load , tip point load , and tip couple load . The straight line and solid curve are the axes of undeflected and deflected (elastica) beams, respectively. The material point of the elastica at is defined in Cartesian coordinates ( , ), where the beam's rotation is and the shear force and bending moment are and , respectively. At the free tip ( = ), the vertical displacement is ∆ , horizontal displacement is ∆ , and rotation angle is . In this study, the beam material was nonlinearly elastic with a stress-strain relationship defined by Ludwick's constitutive law as follows:

Mathematical Modelling
The symbols and loads for the cantilever beam are defined in Figure 2. The beam was subjected to three types of loading, either individually or in combination: a uniform load W, tip point load P, and tip couple load C. The straight line c f and solid curve c f are the axes of undeflected and deflected (elastica) beams, respectively. The material point of the elastica at s is defined in Cartesian coordinates (x, y), where the beam's rotation is θ and the shear force and bending moment are Q and M, respectively. At the free tip (s = l), the vertical displacement is ∆ v , horizontal displacement is ∆ h , and rotation angle is θ f .
The geometry of other types of non-prismatic beams other than the aforementioned linearly varying rectangular cross-section can be easily formed in a similar manner.

Mathematical Modelling
The symbols and loads for the cantilever beam are defined in Figure 2. The beam was subjected to three types of loading, either individually or in combination: a uniform load , tip point load , and tip couple load . The straight line and solid curve are the axes of undeflected and deflected (elastica) beams, respectively. The material point of the elastica at is defined in Cartesian coordinates ( , ), where the beam's rotation is and the shear force and bending moment are and , respectively. At the free tip ( = ), the vertical displacement is ∆ , horizontal displacement is ∆ , and rotation angle is . In this study, the beam material was nonlinearly elastic with a stress-strain relationship defined by Ludwick's constitutive law as follows: In this study, the beam material was nonlinearly elastic with a stress-strain relationship defined by Ludwick's constitutive law as follows: where σ and ε represent the stress and strain, respectively, and E and n represent the Young's modulus and exponential constant of the material, respectively. Note that the sign of σ follows that of ε. The Bernoulli-Euler bending moment-curvature relationship for the rectangular beam with B and H fabricated from a Ludwick-type material can be written as follows [8]: where κ is the curvature dθ/ds of the elastica. The shear force Q is expressed using the free body shown in Figure 2 as follows: Substituting Equations (3) and (4) into Equation (6) yields where The first derivative of Equation (7) with respect to s is obtained as follows: Considering dM/ds = Q, dM/ds is obtained from Equation (8) as Substituting Equation (10) into Equation (9) yields the following ordinary differential equation: The geometric relationships of the Bernoulli-Euler nonlinear beam are defined as follows: Now, we consider the boundary conditions. At the clamped end (s = 0), no deflection occurs, resulting in Note that the boundary condition of dθ/ds at the clamped end is unknown, meaning Equations (11) and (12) cannot be solved as initial value problems. To calculate the unknown dθ/ds as a boundary value problem, the boundary condition at the free tip (s = l) is considered. At the free tip, M in Equation (5) is equal to the tip couple loading C (i.e., M − C = 0), meaning At the free tip, B f and H f are obtained from Equations (2) and (3) as follows: Appl. Sci. 2019, 9, 877 5 of 13 Combining Equations (14) and (15) yields the following boundary condition at the free tip: To facilitate numerical analysis, the load and geometry of the beam are cast into non-dimensional forms. First, the load parameters are defined as The arc length s and coordinates (x, y) are normalized by the span length l as The width and height at the clamped end are normalized by the length l as The tip deflections ∆ v and ∆ h are normalized by the length l, and the tip angle parameter ϕ is defined as follows: Combining Equations (11) and (12) with the non-dimensional parameters in Equations (17)-(19) yields where The boundary conditions in Equation (13) at the clamped end (λ = 0) become The boundary condition in Equation (16) at the free tip (λ = 1) becomes Equations (21)-(25) represent fourth-order differential equations and boundary conditions that govern the elastica (i.e., large deformed shapes of beams) considered in this study.

Solution Methodology
Based on the theories developed above, a FORTRAN computer program was written to compute the elastica. The input beam parameters are the span length l, mechanical properties of E and n, beam geometries of B c , H c , B f , and H f , and individual loads of W, P, and C. These parameters are given in dimensional units but are transformed into nondimensional parameters (i.e., b c , h c , r B , and r H , and p, w, and c). To integrate the differential equations, the Runge-Kutta method [15] was adopted.
To determine the unknown dθ/dλ at the clamped end, the Regula-Falsi method [15] was adopted. The proposed solution methodology can be summarized as follows:

1)
Set the input parameters of n, b c , h c , r B , r H , w, p, and c. 2) Consider a trial (dθ/dλ) t at the clamped end (λ = 0). The first trial is numbered as zero.

3) Integrate Equations (21) and (22) using the boundary conditions in Equation (24) and a trial
(dθ/dλ) t . After completing integration, the trial solution includes (ξ, η, θ, dθ/dλ) in 0 ≤ λ ≤ 1. 4) Calculate D (i.e., the trial boundary condition at the tip (λ = 1)) for Equation (25) as follows: If D = 0, the trial solution computed in step three is the characteristic solution. If |D| ≤ ER 1 , the trial solution is the numerical solution. Here, the first convergence criterion is set to ER 1 = 1 × 10 −8 . 5) If the criterion in step four is not satisfied, increment ∆ to the previous trial (dθ/dλ) 1 to calculate the next trial as (dθ/dλ) 2 = (dθ/dλ) 1 + ∆, and repeat steps 2-5. 6) Note the sign of D 1 × D 2 , where D 1 and D 2 are the values of D in the previous and current steps, respectively. If the sign changes, the characteristic solution of dθ/dλ lies between (dθ/dλ) 1 and (dθ/dλ) 2 . The next trial (dθ/dλ) 3 is defined using the Regula-Falsi method, and a solution for the nonlinear equations is defined as 7) Once the numerical process is entered into the Regula-Falsi scheme, repeat step six until the following criterion is satisfied: where D 1 × D 3 < 0 and the second convergence criterion is set to ER 2 = 1 × 10 −5 .
If the two convergence criteria are satisfied, output the characteristic solution (ξ, η, θ) with the tip responses (δ v , δ h , ϕ) and end calculation.

Results and Comparisons
To analyze typical elastica behaviors, three beam materials were considered in this study: steel, the NP8 aluminum alloy, and annealed copper. For these materials, the values of the Young's modulus E and exponent constant n to be introduced into Ludwick's constitutive law σ = E|ε| 1/n are listed below.
Steel: E = 207 GPa and n = 1.0 NP8 aluminum alloy: E = 455.8 MPa and n = 4.79 Annealed copper: E = 458.5 MPa and n = 2.16 First, the results of this study are compared to those from the literature for degenerate cases, including non-prismatic, nonlinearly elastic, and combined loading effects, for validation purposes. The steel and annealed copper materials are considered in these comparisons. Figure 3 presents comparisons between the tip responses (∆ h , ∆ v ) and (δ h , δ v ) from this study and those from [4,5,8,16]. The beam parameters used for comparison are presented. The tip responses from this study ( ) and those from the literature (•) are presented in the context of the elastica obtained in this study. Both sets of responses are in good agreement for the four degenerate cases. These comparisons validate the elastica theories and solution methods presented in this paper.
Next, we present parametric analysis of the elastica for various combinations of w, p, and c. The results for the annealed copper material (n = 2.16) are presented in Figures 3-5. The beam parameters considered are also presented.

7)
Once the numerical process is entered into the Regula-Falsi scheme, repeat step six until the following criterion is satisfied: where × < 0 and the second convergence criterion is set to = 1 × 10 .
If the two convergence criteria are satisfied, output the characteristic solution ( , , ) with the tip responses ( , , ) and end calculation.   and at two-cycle flip-over (i.e., ϕ = 8) when c = 4.69 × 10 −5 (see loading combinations five and six). Additionally, the elastica under the negative couple load c = −1.79 × 10 −5 5 are presented and show upward deflection with a negative curvature κ (see loading combination seven). For the combined loading case (d), the deflections of the elastica for all seven loading combinations are presented. Load superposition increases the deflection compared to the sums of individual deflections. However, the superposition principle is not satisfied for the nonlinear elastic beams considered in this study. When comparing loading combinations four (w + p), seven (w + p + c), and eight (w + p − 1.45c), one can see that a positive c increases the deflection, and a negative c considerably reduces the deflections caused by the vertical loads w and p. Therefore, the elastica for loading combination eight show upward deflection beyond the horizontal axis.
Next, we present parametric analysis of the elastica for various combinations of , , and . The results for the annealed copper material ( = 2.16) are presented in Figures 3-5. The beam parameters considered are also presented. Figure 4 presents the tip responses of ( , , ) versus the load curves (i.e., equilibrium paths) for three loading cases: (a) uniform load , (b) tip point load , and (c) tip couple load . The trends appear as expected. Specifically, the relationships are strongly nonlinear and ( , , ) all increase as the load increases. For the two loading cases of and , the angle parameter = ( 2 ⁄ ) ⁄ converges to a value of = 1 because the tip rotation cannot physically reach 2 ⁄ based on the vertical loading direction. However, for load can overcome this limitation and exceed = 1 as increases. Additionally, the physical meaning of = 4 is one-cycle flip-over at the tip end, that of = 8 is two-cycle flip-over, etc.  The effects of the cross-sectional ratio r B (= r H ) on the three tip responses of (δ v , δ h , ϕ) and the elastica are presented in Figure 6 for the combined loading case with all of w + p + c. The beam parameters are included in the figure. In Figure 6a, all the values of (δ v , δ h , ϕ) decrease as r B (= r H ) increases. The reduction rate is larger than that with the smaller cross-sectional ratio r B (= r H ). From Figure 6a, one can clearly see that the relationship between deflection and the cross-sectional ratio is strongly nonlinear, as indicated by the equilibrium paths in Figure 4. Additionally, in Figure 6b, the deflections of the elastica become deeper as r B (= r H ) decreases. superposition principle is not satisfied for the nonlinear elastic beams considered in this study. When comparing loading combinations four ( + ), seven ( + + ), and eight ( + − 1.45 ), one can see that a positive increases the deflection, and a negative considerably reduces the deflections caused by the vertical loads and . Therefore, the elastica for loading combination eight show upward deflection beyond the horizontal axis. The effects of the cross-sectional ratio (= ) on the three tip responses of ( , , ) and the elastica are presented in Figure 6 for the combined loading case with all of + + . The beam parameters are included in the figure. In Figure 6a, all the values of ( , , ) decrease as (= ) increases. The reduction rate is larger than that with the smaller cross-sectional ratio (= ). From Figure 6a, one can clearly see that the relationship between deflection and the cross-sectional ratio is strongly nonlinear, as indicated by the equilibrium paths in Figure 4. Additionally, in Figure 6b, the deflections of the elastica become deeper as (= ) decreases. In a third series of experiments, the effects of beam materials on elastica were investigated. Three different beam materials, namely steel, the NP8 aluminum alloy, and annealed copper, were considered. The mechanical properties of and for these materials were provided above. For these parametric experiments, the numerical results are presented in dimensional units, rather than the nondimensional forms, because the load parameters of ( , , ) are different for the same load In a third series of experiments, the effects of beam materials on elastica were investigated. Three different beam materials, namely steel, the NP8 aluminum alloy, and annealed copper, were considered. The mechanical properties of E and n for these materials were provided above. For these parametric experiments, the numerical results are presented in dimensional units, rather than the nondimensional forms, because the load parameters of (w, p, c) are different for the same load magnitudes with the same beam geometry for different materials. The consistent beam geometry and loading conditions for all three beam materials are listed below. Based on the beam parameters listed above, the elastica, bending moments, strains, and stresses of the clamped beam ends (s = 0 m) were computed in the dimensional forms. The results are compared in Figures 7-9. Figure 7 presents the deflections of the elastica for the three beam materials. One can see that deflections increase in the order of the annealed copper to the aluminum alloy to the steel. In particular, the deflections of the copper and aluminum beams are very different, despite the E values of both materials being nearly identical. From these results, we can conclude that the property of n, rather than E, has a dominant effect on the elastica.
Bending-moment diagrams along the beam axis for the three beam materials are presented in Figure 8. One can see that the bending moments decrease from the steel to the aluminum to the copper. It is noteworthy that the bending moments of the steel and aluminum alloy are nearly identical, despite the two materials having very different values of E and n. In contrast, the bending moments of the aluminum alloy and annealed copper are very different. The E values of these two materials are nearly the same, but the n values differ significantly. From these results, we can conclude that both E and n have significant effects when calculating bending moments. Additionally, the maximum bending moments ( ) occur at the clamped ends, whereas the minimum bending moments (•) occur at the free tips. These moments correspond to the tip couple loading c for all three materials. In a third series of experiments, the effects of beam materials on elastica were investigated. Three different beam materials, namely steel, the NP8 aluminum alloy, and annealed copper, were considered. The mechanical properties of and for these materials were provided above. For these parametric experiments, the numerical results are presented in dimensional units, rather than the nondimensional forms, because the load parameters of ( , , ) are different for the same load magnitudes with the same beam geometry for different materials. The consistent beam geometry and loading conditions for all three beam materials are listed below. Based on the beam parameters listed above, the elastica, bending moments, strains, and stresses of the clamped beam ends ( = 0 m) were computed in the dimensional forms. The results are compared in Figures 7-9. Figure 7 presents the deflections of the elastica for the three beam materials. One can see that deflections increase in the order of the annealed copper to the aluminum alloy to the steel. In particular, the deflections of the copper and aluminum beams are very different, despite the values of both materials being nearly identical. From these results, we can conclude that the property of , rather than , has a dominant effect on the elastica.  materials are nearly the same, but the values differ significantly. From these results, we can conclude that both and have significant effects when calculating bending moments. Additionally, the maximum bending moments (■) occur at the clamped ends, whereas the minimum bending moments (•) occur at the free tips. These moments correspond to the tip couple loading for all three materials.   Figure 8. Bending-moment diagrams for steel, aluminum alloy, and annealed copper beams. Figure 9 presents diagrams of the strains and stresses along the beam height at the clamped ends for the three beam materials. The strain and stress at a height are computed as follows: where is the distance measured from the centroid along the vertical beam height. By using the bending moments (■ in Figure 8) at the clamped ends ( = 0 m) for the three beam materials, the strains and stresses were computed using Equations (29) and (30). The results for and are presented in Figure 9a, b, respectively. For the reader's convenience, the extreme values and on both sides of the beams (i.e., (= 2 ⁄ ) = ±5 cm) in Figure 9a, b are tabulated in Table 1 along with the corresponding values of , curvatures , and mechanical properties and . The strains along the height for the three beam materials are presented in Figure 9a. The values of decrease from the annealed copper to the aluminum alloy to the steel. Some unexpected patterns are visible in these results. For example, the strains of the steel and aluminum are nearly identical despite the two materials having very different and values. In contrast, the strains of the aluminum and copper are very different despite both materials having similar values. Figure 9b presents the stresses along the beam height for the three beam materials. The values of decrease from the steel to the aluminum to the copper. The stress distributions of the copper and aluminum near the neutral axis are relatively larger than that of the steel because their exponent constants are larger. However, one can see that the relationships between the stress , and and are reversed compared to the strain relationships shown in Figure 9a, which is an unexpected result (see Table 1). From the results in Figure 9 and Table 1, we can conclude that the effects of and play important roles in computing the strain and stress of nonlinear elastic beams.  Figure 9. Values of (a) strain and (b) stress as functions of the distance from the neutral axis for steel, aluminum alloy, and annealed copper beams at clamped ends.

Conclusions
This paper presented the computation of the elastica of non-prismatic and nonlinear elastic cantilever beams under combined loading. The tested beams had solid rectangular cross-sections with widths and heights varying linearly along the beam axes. The beam materials obey Ludwick's constitutive law. The loading system used for testing can provide uniform loading, tip point loading, and tip couple loading individually or in combination, allowing for a total of seven loading combinations that have not been covered in the literature. We derived ordinary differential equations that govern the elastica of such beams and solved them numerically using the Runge-Kutta and Regula-Falsi methods. To analyze typical elastica behaviors, three beam materials were considered: steel, the NP8 aluminum alloy, and annealed copper. From our parametric study of elastic behaviors, Figure 9. Values of (a) strain ε and (b) stress σ as functions of the distance H from the neutral axis for steel, aluminum alloy, and annealed copper beams at clamped ends. Figure 9 presents diagrams of the strains ε and stresses σ along the beam height at the clamped ends for the three beam materials. The strain ε and stress σ at a height H are computed as follows: where H is the distance measured from the centroid along the vertical beam height. By using the bending moments M c ( in Figure 8) at the clamped ends (s = 0 m) for the three beam materials, the strains ε and stresses σ were computed using Equations (29) and (30). The results for ε and σ are presented in Figure 9a, b, respectively. For the reader's convenience, the extreme values ε e and σ e on both sides of the beams (i.e., H(= H c /2) = ±5 cm) in Figure 9a, b are tabulated in Table 1 along with the corresponding values of M c , curvatures κ c , and mechanical properties E and n. The strains ε along the height for the three beam materials are presented in Figure 9a. The values of ε decrease from the annealed copper to the aluminum alloy to the steel. Some unexpected patterns are visible in these results. For example, the strains of the steel and aluminum are nearly identical despite the two materials having very different E and n values. In contrast, the strains of the aluminum and copper are very different despite both materials having similar E values. Figure 9b presents the stresses σ along the beam height for the three beam materials. The values of σ decrease from the steel to the aluminum to the copper. The stress distributions of the copper and aluminum near the neutral axis are relatively larger than that of the steel because their exponent constants are larger. However, one can see that the relationships between the stress σ, and E and n are reversed compared to the strain relationships shown in Figure 9a, which is an unexpected result (see Table 1). From the results in Figure 9 and Table 1, we can conclude that the effects of E and n play important roles in computing the strain and stress of nonlinear elastic beams.

Conclusions
This paper presented the computation of the elastica of non-prismatic and nonlinear elastic cantilever beams under combined loading. The tested beams had solid rectangular cross-sections with widths and heights varying linearly along the beam axes. The beam materials obey Ludwick's constitutive law. The loading system used for testing can provide uniform loading, tip point loading, and tip couple loading individually or in combination, allowing for a total of seven loading combinations that have not been covered in the literature. We derived ordinary differential equations that govern the elastica of such beams and solved them numerically using the Runge-Kutta and Regula-Falsi methods. To analyze typical elastica behaviors, three beam materials were considered: steel, the NP8 aluminum alloy, and annealed copper. From our parametric study of elastic behaviors, the following conclusions can be drawn:

1)
The tip response of (δ v , δ h , ϕ) increases with an increasing load along the equilibrium path.
Regarding the relationship between tip response and cross-sectional ratio, the tip responses of (δ v , δ h , ϕ) decrease as the cross-sectional ratio r B (= r H ) increases. The tip responses in this study are in good agreement with those in the literature. 2) The mechanical properties of E and n were both observed to have important effects when calculating bending moments. The exponential constant n, rather than the Young's modulus E, dominates elastica calculations. 3) For two loads of w and p, the angle parameter ϕ converges to ϕ = 1 as w, p, and both increase.
However, ϕ for a load c can exceed a value of one. The relationship between deflection and cross-sectional ratio is strongly nonlinear, as indicated by the equilibrium paths. 4) The stress values of copper and aluminum near the neutral axis are relatively larger than that of steel. The stress near the neutral axis increases as the exponent constant n increases.
Author Contributions: J.K.L. proposed the idea, derived the governing equations, and drafted the paper; B.K.L. coded the computer program, obtained the calculations, and assisted the writing of the paper.