A mixed finite element formulation for Timoshenko fiber models based on the local stress field

1. Introduction

Numerous efforts are done to improve upon the accuracy of conventional displacement-based beam-column finite elements for geometric and material nonlinear analysis of frame structures by many researchers. Crivelli and Felippa [1] and Eduardo et al. [2] developed an incremental Total Lagrangian formulation for beam elements that included the effect of large rotation and the effect shear influence. Zhang et al. [3] presented the small strain beam element stiffness by combining Timoshenko beam theory with fiber model in deep beams. Also, Kabeyasawa et al. [4] proposed a macro model called Three Vertical Line Element Model (TVLEM) for numerical modeling of reinforced concrete shear walls and subsequent development of TVLEM was Multiple-Vertical-Line Element Model (MVLEM) presented by Vulcano et al. [5]. Numerical simulations investigating fiber approaches have been carried out by several researchers [6, 7]. In this paper, total Lagrangian approach is chosen among main aspects of kinematic descriptions and an approach is improved to handle the problem of field inconsistency found in 2D Timoshenko planar frame elements for large displacements with small strains and separate uniaxial fibers. Shear behavior of each element is covered with interface elements based on the local stress field model which can be developed for general cases. In the proposed methodology, combination of material nonlinearity and geometric nonlinearity is done through fiber arbitrary constitutive models and local stress field models assigned to interfaces as shown in the following. The present model is a simple applicable approach which is almost accurate and time saving, also it holds proper convergence compared to micro modeling methods.

2. Element formulations

The approach can be applied to finite elements consisting of fiber beam-column elements and interface elements simultaneously. The model has been constructed by combining the constitutive law for each fiber with a set of parallel fiber series. The main purpose of using fiber elements through a set of these series is to improve seismic behavior of the proposed solution with employing shear constraint elements.

2.1. Kinematics of total Lagrangian Timoshenko fiber-based model

In this research the general formulation of total Lagrangian (TL) kinematics is used to derive the finite element equations of a two-node Timoshenko plane beam element. Formulations are implemented based on Timoshenko model with local first-order shear deformation effects. Also, theory can account for geometrically nonlinear behavior due to large displacements and rotations beside material nonlinearity due to constitutive behavior of each fibers. Under the assumptions of the Euler-Bernoulli beam, the following relations can be obtained as Eq. (1). $u_X$ and $u_Y$ are the displacements of the projection of ${\mathrm{\Omega }}_0$ on the neutral axis in the $X-Y$ system. In addition, $\theta$ represents the rotation of the cross-section, $s$ as the arc length along the neutral axis in the current configuration and primes denote the derivatives with respect to $X$ as shown in Fig. 1.

Fig. 1Coordinates and elemental kinematics in the TL formulation

The deformation gradient matrix $\left(F\right)$ and displacement gradient matrix $\left(G\right)$ are derived from Eq. (1):

1

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}X+u_x-Y\sin \theta \\ u_y+Y\cos \theta \end{array}\right]\therefore F=\left[\begin{array}{cc}1+u_x^{\text{'}}-Y\kappa \cos \theta & -\sin \theta \\ u_y^{\text{'}}-Y\kappa \sin \theta & \cos \theta \end{array}\right]\therefore G=F-I.$

Strain-displacement relations can be written according to Eq. (1) and the plane portion of the Green-Lagrange (GL) strain tensor (${}_0{}^2{\epsilon }_{ij}$) is obtained. In the following, a consistent-linearization technique is used by definition an orthogonal matrix $(\mathrm{\Lambda }=\mathrm{\Gamma }+I)$. Therefore, the Green strain matrix can be rewritten and replaced by the simpler form as Eq. (2):

2

$2\epsilon =\left({\overline{F}}^T\overline{F}-I\right)=\left({\left(\mathrm{\Gamma }+I\right)}^T\left(\mathrm{\Gamma }+I\right)-I\right)={\mathrm{\Gamma }}^T+\mathrm{\Gamma }+{\mathrm{\Gamma }}^T\mathrm{\Gamma }\stackrel{{\mathrm{\Gamma }}^T\mathrm{\Gamma }\approx 0}{\Rightarrow }\epsilon =\mathrm{\Gamma }.$

Hence, the axial strains $\left({\epsilon }_{xx}\right)$ and shear strains $\left({\epsilon }_{yx}\right)$ can be arranged in a strain vector $\left(\stackrel{-}{\epsilon }\right)$ (Eq. (3)):

3

$\stackrel{-}{\epsilon }=\left\{\begin{array}{c}{\epsilon }_{xx}\\ {2\epsilon }_{yx}\end{array}\right\}=\left\{\begin{array}{c}e-Y\theta \text{'}\\ \gamma \end{array}\right\}\therefore \delta \mathcal{H}={\left[\begin{array}{c}\delta e\delta \gamma \delta \kappa \end{array}\right]}^T=\stackrel{-}{B}{\left[\begin{array}{c}\delta u_{X1}\delta u_{Y1}\delta {\theta }_1\delta u_{X2}{u}_{Y2}\end{array}\delta {\theta }_2\right]}^T.$

In the following, strain displacement relation can be obtained by nodal displacement variations:

4

$\stackrel{-}{B}=\frac{1}{L_0}\left[\begin{array}{cccccc}-\mathrm{c}\mathrm{o}\mathrm{s}\omega & -\mathrm{s}\mathrm{i}\mathrm{n}\omega & \gamma {\phi }_{4,x}\left(x\right)L_0& \mathrm{c}\mathrm{o}\mathrm{s}\omega & \mathrm{s}\mathrm{i}\mathrm{n}\omega & \gamma {\phi }_{6,x}\left(x\right)L_0\\ \mathrm{s}\mathrm{i}\mathrm{n}\omega & -\mathrm{c}\mathrm{o}\mathrm{s}\omega & -(1+e){\phi }_{4,x}\left(x\right)L_0& -\mathrm{s}\mathrm{i}\mathrm{n}\omega & \mathrm{c}\mathrm{o}\mathrm{s}\omega & -(1+e){\phi }_{6,x}\left(x\right)L_0\\ 0& 0& -1& 0& 0& 1\end{array}\right].$

Using the shape functions of the reference element, the discretization of the variables is performed. The elements with distributed nonlinearity were formulated with the classical stiffness method using cubic Hermitian polynomials to approximate the deformations along the element. According to Hermite interpolation the degrees of freedom for each element are the displacements in orthogonal directions and slopes at the two nodes:

5

$\theta \left(x\right)=\left[\begin{array}{cccc}{\phi }_{3,x}\left(x\right)& {\phi }_{4,x}\left(x\right)& {\phi }_{5,x}\left(x\right)& {\phi }_{6,x}\left(x\right)\end{array}\right]*{\left[\begin{array}{cccc}{v}_1& {\theta }_1& v_2& {\theta }_2\end{array}\right]}^T,$

where ${\mathrm{\Phi }}_z$ is presented as shear slenderness which defined as Eq. (6):

6

${\mathrm{\Phi }}_z=1+\frac{12E{I}_z}{{\kappa }_{y}GA{L}^2}\therefore E{I}_z={\sum }_{n_f=1}^{n_{tf}\left(x\right)}{E}_{n_f}{y}_{n_f}^2,$

where ${\kappa }_y$ is defined as shear stress correction factor which assumed 0.833 for rectangular sections. The variable $G$ is shear modulus obtained by a subroutine program based on local stress field direct analysis which detailed explanation is given in Section 3.

The $\stackrel{-}{B}$ is the matrix giving the generalized strains according to nodal displacements in each loading step. It is essential to compute sectional resistance force based on the given stress of the fiber state. The section resistance forces $\left({\mathcal{F}}_s\right)$ are computed from the fiber stress distribution and section stiffness $\left(K_s\left(x_i\right)\right)$ is assembled from the fiber stiffnesses. This force is implemented according to Eq. (7), and the following stages lead to the calculation of element resistance force:

7

$K_s\left(x_i\right)=\left[\begin{array}{ccc}{\sum }_{n_f=1}^{n_{tf}\left(x\right)}{E}_{n_f}{A}_{n_f}& 0& -{\sum }_{n_f=1}^{n_{tf}\left(x\right)}{E}_{n_f}{A}_{n_f}{y}_{n_f}\\ 0& {\sum }_{n_f=1}^{n_{tf}\left(x\right)}{G}_{n_f}{A}_{n_f}& {\sum }_{n_f=1}^{n_{tf}\left(x\right)}{G}_{n_f}{A}_{n_f}{y}_{n_f}\\ -{\sum }_{n_f=1}^{n_{tf}\left(x\right)}{E}_{n_f}{A}_{n_f}{y}_{n_f}& {\sum }_{n_f=1}^{n_{tf}\left(x\right)}{G}_{n_f}{A}_{n_f}{y}_{n_f}& {\sum }_{n_f=1}^{n_{tf}\left(x\right)}{E}_{n_f}{A}_{n_f}{y_{n_f}}^2\end{array}\right].$

Based upon an Total-Lagrangian approach, the incremental matrix equilibrium equation for a system can be formulated using the variation of the increment of total potential energy and governed equations are simplified as:

8

$\delta F_{int}=\delta \left({\int }_{L_0}{\stackrel{-}{B}}^T{\mathcal{F}}_{s}dX\right)={\int }_{L_0}\left({\stackrel{-}{B}}^T{K}_s\stackrel{-}{B}\delta u+\frac{\partial {\stackrel{-}{B}}^T}{\partial u}\delta u{\mathcal{F}}_s\right)dX=\left(K_{Mat.}+K_{Geo.}\right)\delta u,$

where $K_{Mat.}$ and $K_{Geo.}$ are the stiffeness of fiber planar Timoshenko frame element due to material nonlinearity and geometrical nonlinearity, respectively:

9

$K_{Mat.}={\int }_{L_0}\left({\stackrel{-}{B}}^T{K}_s\stackrel{-}{B}\right)dX=\frac{L_0}{2}{\sum }_{i=1}^{n_{GP}}{w}_i{\left({\stackrel{-}{B}\left(x\right)}^T{K}_s\left(x\right)\stackrel{-}{B}\left(x\right)\right)}_{{x=\frac{L_0}{2}(1+x}_i)},$

10

$K_{Geo.}={\int }_{L_0}\left(\frac{\partial {\stackrel{-}{B}}^T}{\partial u}{\mathcal{F}}_s\right)dX={\int }_{L_0}\left(\frac{\partial {\stackrel{-}{B}}^T}{\partial u}\left\{\begin{array}{c}\begin{array}{c}N\\ V\end{array}\\ M\end{array}\right\}\right)dX.$

3. Interface element implementation

A basic interface element which is consisted three springs at each pairing of mating nodes is based on beam column degrees of freedom. Springs are separated in classification of axial, rotating and shear stiffness which is represented the interface behaviors. The stiffness matrix has the following shape:

where $k_t$, $k_n$, and $k_s$ are stiffness of springs in shear, axial, and moment, respectively. The shear stiffness is earned through relative tangential displacements between the two nodes $(x\mathrm{},\mathrm{}y)$ from displacement vector in the local coordinate and local stress field direct analysis approach (See Section 3).

4. Constitutive law in total strain form

Under uniaxial tension-compression or biaxial compressive stress state of the material, nonlinear material characteristics are considered [8-10]. In the RC structures, the constitutive relation in the uncracked state is restricted to linear elasticity. The constitutive shear model is adapted according to Li [11] which was developed for modeling the nonlinear behavior of concrete elements. In this paper, local stress field theory based on fixed smeared crack approach is implemented in the subroutine of the program. Some details can be found in [12]. Nonlinear element analysis is used in the form of a subroutine in the main program, and these local interface elements are assembled along with fiber beam column elements. Local analysis is performed in each computational step state. The nonlinear behavior of the element can be obtained using an iterative solution method by applying in-plane incremental stresses or strains with secant elasticity modulus $\left({\left[E\right]}^G\right)$ as Eq. (12). The analytical process is executed by implementation of direct displacement control that was introduced by Jirásek and Bazant [13]:

Based on mentioned scheme, the applicability of chosen model for different types of structures with arbitrary constitutive behavior is implemented.

5. Numerical implementation and solution

Some numerical examples are used to verify the accuracy and show the efficiency of the proposed geometrically nonlinear frame element as well as the solution marching schemes. The iterative-incremental method (Arc-Length method) with a variable stiffness scheme was applied to analyze structures. Afterward, several numerical investigations were performed with the proposed model in order to study the effects of nonlinear shear deformations and flexural responses, simultaneously.

Fig. 2Load-vertical displacement responses of cantilever beam

Fig. 3Load-vertical displacement responses for Williams’ Toggle frame

The results of the nonlinear computer analyses are compared with the observed data and analytical results. These examples are considered as “classic” benchmark since some of them have been widely used by several researchers to evaluate their nonlinear frame models. At first, A cantilever beam subjected to an end-point load [14] shown in Fig. 2, Williams’ Toggle frame [15] (Fig. 3) and Lee et al. frame [14] (Fig. 4) are employed.

The ability of the proposed method to predict the nonlinear behavior of the structures is validated in this section by comparing the analytical and experimental results of the tests performed by Yaghoubifar’s rehabilitations [16] and Chun and Kim’s RC joint connection [17]. The SSBW2 was strengthened on one side of the masonry wall with concrete.

Fig. 4Load-displacement responses for Lee et al.’s frame: a) load-vertical displacement diagram and b) load-horizontal displacement diagram

a)

b)

Fig. 5Load-displacement responses for strengthened masonry walls a) SSBW2 and b) Chun and Kim’s RC joint connection

a)

b)

6. Conclusions

In this paper, the ability and accuracy of the proposed element is illustrated to model a fiber beam-column experiencing considerable large displacements compared to its initial length and showing unlimited hardening behavior which can account for different kinds of fiber constitutive materials (Steel, Masonry, Concrete). The fiber beam-column model proposed in this paper is based on distributed plasticity techniques. The advantage of the distributed-plastic over the plastic-hinge method is its accuracy in defining the plastic zone during the analysis. Moreover, a new formulation to attain the shear response along with flexural behavior for any load cases is suggested as local stress field model implemented in interface elements based on the fixed smeared crack approach. The proposed model was calibrated and validated through a comparison with experimental results and various numerical analyses were performed in order to study the influence of nonlinear flexural-shear interaction. Thereafter, this method could yield accurate and convergent results in agreement with the problems.