1 Introduction

Classical mechanics states that engineering materials have the same elastic moduli in tension and compression. However, a myriad of experimental tests show that there are also such materials as concrete^[1], metal alloy^[2], biological materials^[3], rubber^[4], rock^[5], and nuclear graphite^[6] having different tensile elastic and compressive elastic moduli, which are called bi-modulus materials. Specially, graphene, a novel two-dimensional (2D) nanomaterial with the thinnest thickness but the highest strength, is such a kind of materials^[7-8]. Under this circumstance, a relatively high calculating error may occur if the classical modulus theory is still adopted. In contrast, the bi-modulus theory is able to reflect the constitutive law of materials and the mechanical response of structures accurately.

The concept of bi-modulus materials was firstly proposed by Timoshenko when conducting the mechanical characteristic research of a pure-bending beam^[9]. The extension of the bi-modulus theory centers on the constitutive model, the finite element method (FEM), and the analytical method.

The bi-modulus constitutive model remedies the strain error of materials when subjected to the tensile stress and the compressive stress with the same absolute values. Ambartsumyan^[10] introduced two straight line models to describe the constitutive relationship of bi-modulus materials. Jones^[11] proposed two modified models considering Possion's effect by introducing the arithmetic mean value and weighting coefficient^[12], respectively. Besides, Jones^[13] and Ye^[14] proposed a principal strain model and a bilinear biaxial principal stress model, respectively.

The FEM is initially introduced to study the mechanical response of the structures composed of bi-modulus materials. Several dominant methods have been proposed in the last five decades, e.g., the initial stress method^[15], the accelerating convergence factor, the equivalent shear modulus method^[16-17], the assumption µ⁺/E⁺=µ^-/E^-[18], the smoothing function method^[19], and the double shear modulus method^[20].

The analytical method is able to reflect the accurate relationship among the sensitive parameters, e.g., the tensile modulus, the compressive modulus, and the mechanical response of bi-modulus structures. Ambartsumyan^[10] established numerous analytical models for the bi-modulus structures subjected to a single load. Yao and Ye demonstrated that the position of the neutral layer was unrelated to the shear stress, and derived the position solutions of the neutral layer, stress, and displacement of a bi-modulus bending beam^[21] } and a bending-compression column^[22] under a random load and different boundary conditions. They provided methods for the internal force of the bi-modulus statically indeterminate structures^[23] and composite structures under compress stresses^[24] with nonlinear iterative processes. He et al.^[25] adopted the simplified equivalent section method, and obtained the approximate elastic mechanical solution of a bi-modulus bending-compression column and a bending beam. Qu^[26] studied the deformation of the geocell, considering the bi-modulus property by introducing the concept of an equivalent elastic modulus. Leal et al.^[27] offered a compressive strength equation of a bi-modulus high performance fiber. He et al.^[28] derived the analytical expressions of the internal force and deformation of a bi-modulus bending thin plate. Yao et al.^[29] and Yao and Ma^[30] introduced several dimensionless parameters and a compiled nonlinear iterative program, and analyzed the bi-modulus buckling stability of constant and variable cross-sections under the infinitesimal deformation theory. He et al.^[31] deduced the analytical expressions of bi-modulus curve bending beams with different sections.

However, few researches have been conducted to study the mechanical response of the bi-modulus structures placed in the temperature field. Merely, Kamiya^[32] derived the analytical expressions of the stress and deformation of the bi-modulus cylinder in the nuclear reactor. Yao and Ye^[33] studied the temperature stress field in the bi-modulus single frame, and found that the concrete bi-modulus property deduced by solar radiation might have a significant effect on the mechanical response of the pavement and bridge. In view of this, a bi-modulus beam placed on the Winkler foundation is studied. A linear temperate field is distributed along the beam height, and the mechanical response of the bi-modulus beam on the Winkler foundation is studied to identify the calculation error between the bi-modulus theory and the classical same modulus theory.

2 Material model for different moduli

In this paper, Ambartsumyan's bi-modulus constitutive model is adopted to describe the relationship between the stress and the strain of the material^[10]. As shown in Fig. 1, the tensile modulus and the compressive modulus of the material are E_p and E_n, respectively.

3 Structural model 3.1 Basic assumption

The bi-modulus beam studied in this paper is assumed to be continuous, homogeneous, and isotropic. Subjected to a certain load, merely small elastic deformation will occur in the beam. Therefore, the equilibrium equations, geometric equations, and deformation continuity equations are the same as those in classical mechanics. However, the constitutive equations are totally different from the same modulus theory.

Simultaneously, it is assumed that the beam is placed on the Winkler foundation and the temperature varies linearly along the height of the beam.

3.2 Structural model

As shown in Fig. 2, a bi-modulus beam (2l ×b×h) is placed on the Winkler foundation with the resistance coefficient k. The thermal expansion coefficient of the material is α. Let the origin be on the middle of the top surface of the bi-modulus beam, and establish a rectangular coordinate system. In the foundation beam, the temperature distribution function can be described as follows:

(1)

where T₀ is the temperature at the bottom surface of the bi-modulus beam, ΔT is the temperature difference between the top surface and the bottom surface, and G is the temperature gradient.

4 Theoretical analysis and analytical derivation 4.1 Establishment of governing equations

Considering that the symmetry of two sides is along the y-axis, only half of the foundation beam can be simplified. According to the plane cross-section assumption, the normal strain of the random cross-section of the bi-modulus beam can be obtained as follows:

(2)

where α is the thermal expansion coefficient of the material.

Therefore, the constrained strain of the bi-modulus beam can be obtained as follows:

(3)

It is clear that the bi-modulus beam will be divided into two regions along the y-axis, i.e., the tensile region and the compressive region. The interface of the two regions is called the neutral surface or neutral axis, on which the strain is zero. Set the expression of the constrained strain to be zero. Then, we have

(4)

As shown in Fig. 2, y₀ is indeed the position of the neutral axis. Therefore, the upper region of the bi-modulus beam is the compressive region, and the bottom region is the tensile region. The stress in the bi-modulus beam can be written as follows:

(5)

where A(x), B(x), and y₀ are three unknown quantities. However, only Eq. (4) is available, and the condition is insufficient. Therefore, the equilibrium equations of the axial force and bending moment are required.

According to the bi-modulus Winkler foundation beam model established in this paper, the following equilibrium equations in the horizontal orientation can be obtained:

(6)

(7)

4.2 Derivation of undetermined parameters of normal strain Substituting Eq. (4) into Eq. (6)

Substituting Eq. (4) into Eq. (6), we have

(8)

Then, substituting Eq. (8) into Eq. (7), we have

(9)

As shown in Fig. 3, an infinitesimal segment dx is selected from the bi-modulus beam. In line with the equilibrium condition, we have

(10)

Considering the Winkler assumption, we get

(11)

where v is deemed to be positive when it points down.

The geometric equations can be written as follows:

(12)

Then, we have

(13)

Substituting Eq. (11) and Eq. (13) into Eq. (10) and simplifying it yield

(14)

Taking the derivative of both sides of Eq. (14), we get

(15)

Substituting Eq. (9) into Eq. (15) and making some simplifications, we can get the following differential equation of the undetermined parameter B(x):

(16)

where

Significantly, when E_n=E_p=E, we get . Equation (16) can be returned to the differential equation of the undetermined parameter B(x) under the classic same modulus theory.

Solving Eq. (16), we can get

(17)

where C₁, C₂, C₃, and C₄ are undetermined constants.

Integrating Eq. (14) and invoking the boundary conditions

(18)

we have

(19)

The unknown constants can be solved as follows:

(20a)

(20b)

(20c)

(20d)

Substitute Eq. (20) in sequence into Eq. (17) and Eq. (8). Then, we can deduce the undetermined parameters of the normal strain, i.e., A(x) and B(x), as follows:

(21)

(22)

4.3 Derivation of normal stress, bending moment, and displacement

Substituting Eq. (21) and Eq. (22) into Eq. (5) yields the normal stress as follows:

(23)

Substituting Eq. (21) and Eq. (22) into Eq. (9) yields the bending moment as follows:

(24)

Substituting Eq. (24) into Eq. (10) and Eq. (11) yields the displacement as follows:

(25)

5 Development of numerical calculation procedure

In this paper, a new numerical calculation procedure is developed for calculating the stress field in the bi-modulus structure subjected to a temperature effect by adopting the ABAQUS user subroutine UMAT^[34]. The task of the procedure is to finish the updating of the Jacobian matrix and stress in the end of each increment.

In the temperature field, the total strain of the bi-modulus structure in the increment form can be written as follows:

(26)

where dε is the total strain increment, dε^e is the elastic strain increment, and dε^T is the thermal strain increment. The corresponding expressions are

(27)

(28)

(29)

The stress increment of the bi-modulus structure can be written as follows:

(30)

Considering the constitutive law of bi-modulus materials, the relationship between the stress and the strain in the increment form can be obtained as follows:

(31)

where D is the Jacobian matrix defined by

(32)

In the above equation,

(33)

where E₁, E₂, E₃, v₁, v₂, and v₃ will be identified, considering the symbols of the three principal stresses σ₁₁, σ₂₂, and σ₃₃. Since one principal orientation is random, if the sign of the principal stress in this orientation is positive, the elastic modulus in this orientation is the tensile modulus E_p, and Possion's ratio is the tensile Possion's ratio v_p; otherwise, the elastic modulus is the compressive modulus E_n, and Possion's ratio is the compressive Possion's ratio v_n. Simultaneously, the equation v_p/E_p=v_n/E_n will be satisfied.

The equivalent shear modulus is calculated by introducing the accelerating convergence factor η^[16]. When σ₁₁ < 0, σ₂₂ < 0, and σ₃₃ < 0, η =0. When σ₁₁ >0, σ₂₂ >0, and σ₃₃ >0, η =1. When σ₁₁σ₂₂σ₃₃ < 0, η is the ratio of the sum of absolute values of positive σ_ii (i=1, 2, 3) to the sum of absolute values of σ₁₁, σ₂₂, and σ₃₃, ranging from 0 to 1. For example, if σ₁₁ >0, σ₂₂ >0, and σ₃₃ < 0, then .

Therefore, we can get the expression of the equivalent shear modulus as follows:

(34)

During the whole loading process, E_p, E_n, and α are all constants with the time increment change. Therefore, we have

(35)

Substituting Eq. (35) into Eq. (31), we get

(36)

In the procedure, five parameters are defined, i.e., the tensile modulus E_p, the compressive modulus E_n, the tensile Possion's ratio v_p, the compressive Possion's ratio v_n, and the coefficient of thermal expansion α. Besides, a status variable is defined to store and update the increment of the thermal strain at the end of each incremental step. Through the continuous updating of the Jocabian matrix and the stress of each increment step, the stress field in the bi-modulus structure can be finally obtained. The flow chart of the calculation procedure is shown in Fig. 4.

6 Examples and discussion

Consider a bi-modulus beam placed on the Winkler foundation (see Fig. 2). The detailed dimensions of the beam are 2l×b×h =4.8 m×0.3 m×0.6 m. The thermal expansion coefficient of the material is α =8 × 10⁻⁶ (℃)⁻¹, and Possion's ratio is v=0.18. The temperature at the bottom of the beam is T₀=60 ℃. The temperature of the top surface is lower than that at the bottom, and the difference is ΔT₀=45 ℃. The resistance coefficient of the Winkler foundation is k=20 N/cm³.

A different elastic modulus is assumed for the following three cases:

Case Ⅰ E_n=2.4 × 10⁷ kN/m², m_n ∈ [1,5];

Case Ⅱ E_p=2.4 × 10⁷ kN/m², m_n ∈ [1,5];

Case Ⅲ ,

where .

Four methods are adopted to calculate the temperature stress of the bi-modulus beam, i.e., the analytical method based on the same modulus theory, the analytical method based on the bi-modulus theory proposed in this paper, the numerical simulation processed by the traditional FEM software ABAQUS, and the numerical method processed by the FEM procedure developed in this paper. The results of the neutral axis position and the maximum stress of the beam are listed (see Table 1, partial results list only).

As shown in Table 1, the bi-modulus analytical solution derived in this paper can completely recover the classical same modulus solution when E_p=E_n. Simultaneously, the computation error between the analytical solution and the numerical solution is less than 3%. The computation error may be caused by the grid density, the iterative process, the round-off error of the terminal value, and the density and stiffness of the spring element adopting to simulate the Winkler foundation. Therefore, it is verified that the computation is acceptable, and the computing methods adopted in this paper are valid and rational.

When the tensile and compressive moduli change, there is a regular variation of the neutral axis position of the bi-modulus beam. The neutral axis lines are, respectively, labeled in Figs. 5 and 6 with the lines σ =0 to certify the positions of the neutral axis in the bi-modulus beam in each case. When E_n increases, the height of the compressive region decreases, and vice versa. At the same time, the deflection rate of the neutral axis in the bi-modulus beam decreases when E_n/E_p increases.

As shown in Figs. 5 and 6, after including the bi-modulus property of the material, the variation law of the normal stress on the cross-section of the beam behaves discrepantly to the same modulus solution. Although the maximum tensile and compressive stress occur still at the bottom and top of the beam, normal stress distributes symmetrically beside the neutral axis in the same modulus theory. However, according to the bi-modulus theory, the normal stress distribution behaves unsymmetrically. Therefore, the neutral axis divides the beam into a tensile region and a compressive region. The stress values are not the same at the corresponding points in the tensile and compressive regions.

The bi-modulus property has a significant effect on the maximum stress in the foundation beam. For the concrete material widely used in the road and bridge engineering, E_n/E_p is approximately 2.5. Figure 7 shows the relationships among E_n/E_p, the maximum tensile stress difference between the values calculated by the bi-modulus theory and the classical same modulus theory δσ_max^p, and the maximum compressive stress difference between the values calculated by the bi-modulus theory and the classical same modulus theory δσ_maxⁿ. As shown in Fig. 7, σ^p _max calculated by the bi-modulus theory is 21% larger than that calculated by the classical same modulus theory. Therefore, excessive deformation may appear, and more and more oversize cracks may generate correspondingly.

As depicted in Fig. 8, σ_max^p decreases when E_n/E_p increases, which indicates that increasing E_n will deduce the tensile stress. Analogously, raising E_p will reduce the compressive stress (see Fig. 9).

In Fig. 10, when E_n/E_p aggrandizes from 1 to 5, E_n fortifies and σ_maxⁿ increases accordingly. When E_n/E_p∈ [0.2, 1.0], E_n dwindles, and σ_maxⁿ abates. Significantly, the irregularity of the section stiffness reduces at the same time, which leads to the attenuation of the integrated average stiffness of the section, and then σ_maxⁿ diminishes. When the two kinds of damping effects on σ_maxⁿ superpose, the variation curves of σ_maxⁿ against E_n/E_p with the limits of 0.2 to 1 become steeper than the curves when E_n/E_p changes from 1 to 5. In Fig. 11, when E_n/E_p ranges form 1 to 5, when E_n/E_p magnifies, E_p declines, while σ_max^p lessens, but the curves behave more placid compared with the cases when E_n/E_p ranges from 0.2 to 1. However, the cumulative curves illustrate the change rules of σ_maxⁿ and σ_max^p against E_n/E_p, i.e., when E_n/E_p increases, σ_maxⁿ decreases, while σ_max^p increases, which is absolutely concordant with the law of the stiffness regulating internal force.

As shown in Figs. 8-11, the maximum of the normal stress emerges at the middle cross-section and along the beam span, and the normal stress reduces gradually to 0.

Since the total stiffness of the cross-section is constant, E_n/E_p changes. In Figs. 12 and 13, the displacement of the beam at x=0 increases when E_n/E_p increases. It indicates that when the bi-modulus characteristic is considered, compared with the same modulus theory, the displacement increases as a result of the non-uniformity of the cross-section stiffness, and the resistance to the deformation of the foundation beam decreases.

7 Conclusions

In this paper, based on the bi-modulus theory and the Winkler foundation assumption, the governing equations about the neutral axis position and undetermined parameters of the normal strain of the foundation beam subjected to the linear temperature are firstly established. Then, the calculating formula of the normal stress, displacement, and bending moment of the bi-modulus beam are derived. A temperature stress analysis procedure of bi-modulus structures is developed, and the numerical solution of the temperature stress of a bi-modulus beam placed on the Winkler foundation is obtained. The analytical solution can completely recover the same modulus solution, and the computational difference between the analytical solution and the numerical solution is less than 3%, which indicates that the analytical solution and the temperature stress procedure are both valid and rational. In the end, the discrepancy between the bi-modulus analysis and the classical same modulus analysis is discussed. The detailed conclusions are as follows.

(ⅰ) When E_n>E_p, as E_n/E_p increases, the neutral axis of the foundation beam offsets to the bottom, leading to the attenuation of the compressive region and the enlargement of the tensile region, and vice versa. Along the height of the section, the stress distribution turns to be bilinear.

(ⅱ) For the representative concrete material widely used in road and bridge engineering, σ_max^p calculated by adopting the bi-modulus theory increases by 21% compared with the classical same modulus solution. In practical engineering, it is worth noting that the cracks and damage of pavement might be a matter of ignoring the bi-modulus properties of the concrete structures. Therefore, it is necessary to introduce the bi-modulus theory to analyze the mechanical responses of the bi-modulus structures under specific operating conditions so as to modify the error generating according to the classical same modulus theory.

(ⅲ) Enhancing the compressive modulus can effectively decrease the tensile stress in the bi-modulus beam placed on the Winkler foundation. Thus, the rational selection of materials with a certain difference between the tensile modulus and the compressive modulus in practical engineering can improve the stress state of structures.

(ⅳ) Compared with the classical same modulus theory, the structure stiffness discreteness increases, and the resistance to the deformation of the foundation beam decreases accordingly after considering the bi-modulus characteristic of the material.