Special Finite Elements with Adaptive Strain Field on the Example of One-Dimensional Elements

Abstract

The conception of special finite elements called multi-area elements for the analysis of structures with different stiffness areas has been presented in the paper. A new type of finite element has been determined in order to perform analyses and calculations of heterogeneous, multi-coherent, and layered structures using fewer finite elements and it provides proper accuracy of the results. The main advantage of the presented special multi-area elements is the possibility that areas of the structure with different stiffness and geometrical parameters can be described by single element integrated in subdivisions (sub-areas). The formulation of such elements has been presented with the example of one-dimensional elements. The main idea of developed elements is the assumption that the deformation field inside the element is dependent on its geometry and stiffness distribution. The deformation field can be changed and adjusted during the calculation process that is why such elements can be treated as self-adaptive. The application of the self-adaptation method on strain field should simplify the analysis of complex non-linear problems and increase their accuracy. In order to confirm the correctness of the established assumptions, comparative analyses have been carried out and potential areas of application have been indicated.

1. Introduction

The first step in the finite element method (FEM) is the discretization of the domain into finite elements. The meaning of this process is essential especially if the influence of the results is concerned. To ensure that the simulation of the original domain is suitable to perform the calculations, such parameters as shape, size, type, and configuration of the elements have to be taken into consideration. Another parameter is the number of the finite elements, which has significant influence of the computational effort.

In practice of structural modeling using FEM [1,2,3] very often it is necessary to define domains with areas of significantly different geometrical and material characteristics. The obvious solution is to divide these areas into subregions (sub-areas) and describe them with separate finite elements. This approach is natural for FEM, however it can be computationally expensive (a large number of finite elements, complicated discretization algorithms, etc.). In particular it refers especially to materials with high heterogeneity.

A typical example of the occurrence of sub-areas with significant stiffness differences is the analysis of physically non-linear processes. In such a case, stiffness in the analysed points can even reach the extremal values, e.g., yielding or cracking. Simultaneously the changes in stiffness are not followed by the adjustment of the strain field. This is one of the main methodological reasons of non-linear analysis error generation. In this case, it could be reasonable to use special type of elements called “multi-area elements”, where sub-areas with different stiffness parameters can be combined and completed with a self-adaptive procedure of deformation field. It means that one special multi-area element can be used to describe larger (wider) parts of structure even with different geometry or material parameters.

Special multi-area elements have been developed to deal with problems in analysis of structures with different geometrical and stiffness parameters for example inclusions or accidental damages in form of various holes, etc. Moreover multi-area elements can be used in analysis of constructions with openings, for example wall panels with openings for doors or windows. The main goal of using such an element is the limitation of element and nodes number in comparison with the fine mesh finite element model [4,5].

Currently, in the literature the significant number of scientific papers is dedicated to searching for a new finite elements that ensure the appropriate accuracy of the results and simultaneously reduce the number of unknowns in the computational model. They are often called dedicated elements, developed in order to analyse a specific physical phenomenon. Such type of elements are exactly the special multi-area elements, formulated in order to analyse structures and structural systems containing areas of different geometrical and stiffness characteristics.

Generally, searching for dedicated elements used in the numerical analysis of building structures is carried out at different levels of observation scales. These issues concern the macroscale and microscale. For example, macromodeling technique is used for large structures and consists in describing by finite element a fragment of a construction with a complex structure, for example part of a brick wall containing several bricks and a mortar between them, like in [6,7]. Similar approach was represented by Choi and Bang [8]. Their plate rectangular element has been used in analysis of shear walls with openings. The stiffness matrix of such an element has been developed by subtracting the stiffness of opening from the stiffness of whole plane of stress element without opening. The main disadvantage of this method was the increase of calculation error due to increase of the openings size. Another type of element used in the analysis of structures with different stiffness areas are super elements presented by Kim and Lee [9]. These elements have been also applied in structural analysis of shear walls with openings. Typical for super elements is elimination of nodes using the matrix condensation technique. That is why super elements have only four nodes at the corners of a wall panel, but to satisfy the compatibility conditions at the boundaries of the element, so called fictitious beams have to be added. A certain limitation of modelling the structure with use of macromodeling technique is the fact that homogenization of the material properties is required, which is performed on the basis of theoretical assumptions and the results of experimental research.

However, micromodeling is based on the analysis of small-sized parts of the constructions, including microscopic structures, and relates only to selected areas of building constructions. Similar to macromodelling, the final finite element is the result of homogenization of various sub-areas. Additionally, the interaction between the components of the microstructure is taken into account, introducing complex interface elements and/or describing non-linear phenomena such as plasticity or cracking [10]. Mostly analysis of microstructure is performed in order to assess influence of microdefects on bearing capacity of building components, on strength of materials and also on potential progress of failure mechanism [11]. A characteristic attribute of the micromodeling technique is the development of a computer model of a single cell of the structure, and then homogenization [12,13,14] Unfortunately, analysis of microstructures with this technique is computationally expensive, which often results in the lack of practical application of this approach in the calculations of large and complex structural models.

As heterogeneous materials are largely treated, the materials containing various types of inclusions like in works [15,16] or holes. Very often, the analysis of structures containing holes or inhomogeneities is the part of the so-called multi-scale modelling. In recent decades, there has been a significant development of methods that take into account the multi-scale nature of materials. An example can be the work [17], where the authors describe the developed element containing an opening and a crack simultaneously. In this work, according to the presented examples, the authors confirm the reduction of finite element number and prove that it has an influence on improved numerical calculations and obtaining a satisfactory accuracy of the results. This is in line with the current trend of searching for methods aimed at limiting the number of finite elements used in calculations. As another example, we mention the work [18], where the numerical implementation of the element with an elliptical hole and crack with a different number of nodes on the edge of the element was presented. In these elements, Trefftz functions were used as test functions for stresses and strains that meet the appropriate differential equations and boundary conditions in advance. These functions make it possible to calculate the stiffness matrix of a finite element by means of integration along the border (edge) of a finite element.

Wang and Qin [19] used a hybrid elements with elliptical holes in the calculations of plates in a plane stress state. A plate with a single horizontal, sloping elliptical hole and with two horizontal elliptical holes was analysed. In the work [20], the modified stress functions were presented, that were used in the calculation of the stress distribution in a constrained (finite) plate with a rectangular hole subjected to uniaxial tension. A finite area containing a rectangular opening has been mapped to a finite area containing a single circular opening using the so-called mapping function.

The aim of development and implementation of all mentioned elements and also presented in the paper multi-area elements is striving for reduction of elements and nodes number in calculation models.

The basic type of possible finite element in which the conception of developed multi-area element can be implemented is 1-D double-hinged linear element. An example of a one-dimensional element in which the multi-area elements can be applied has been shown in Figure 1b. This type of elements was chosen because it is the clearest and understandable way to explain the presented method.

To obtain the multi-area element, a simple technique of stiffness summation [21] has been applied, according to the Formula (1).

$$K_e={\sum }_{k=1}^n{K}_k(\chi ),\mathrm{where}\chi =f(\sigma ,E,t\dots )$$

(1)

Expression (1) shows that the stiffness redistribution inside a finite element is further modified by the parameter χ. The parameter χ for example can be a function of the state of effort σ, material parameter such as Young’s modulus E or element geometry such as thickness t. The explanations and determination of parameter χ have been described in Section 3. In Equation (1) index k denotes the number of particular sub-area, n means total number of sub-areas included in finite element.

At this point it is necessary to comment simply technique of the stiffness summation, which is often used in similar cases. As the authors show in the next section, this technique is inadequate, especially when the differences in stiffness of sub-areas are considerable. Hence the proposed method described by Formula (1) has been developed.

2. Research Problem

As an object of analysis, damaged posttensioned concrete girder KBOS-24 [22] has been chosen. In this case, cracks were observed between the column and the top chord of girder, which has been shown in Figure 2b.

Analyzing the problem, it was found that the valid girder model can and should take into account hinge connections at both ends of column. This means at the places of connection with rigid top and bottom chord. In Figure 2b, an extreme case of concrete discontinuity between the column and the top chord of the girder has been shown. The same damages (but much smaller) have been also noticed on the other columns of the girder KBOS-24 as well. According to the research, these cracks are mainly caused by the technology of realization. It should be mentioned that such a types of girders were produced and assembled in Poland in the 60’s of the last century, when the Polish economy was still recovering from the war damage.

Potential calculation models have been illustrated in Figure 3. First model, presented in Figure 3a, is a typical 3D model made in professional computer program. The second, shown in Figure 3b, is the computational model that was made using the “ORCAN” Structure Analysis System (Version 0.98, Bialystok University of Technology, Bialystok, Poland) [23]. The former model, apart from a lot advantages, also has disadvantages like quite problematic geometry or difficulties in finding the internal forces (shear or normal). The latter model is simplified but based on calculations made during the design process of posttensioned girders. In this model the columns were modeled using multi-area double-hinged linear elements. The theoretical assumptions of multi-area elements have been described in the next chapter.

In case of construction damage shown in Figure 2b, analysis of the crack propagation can be done in many different methods. The most obvious is the rearrangement of the mesh discretization in the calculation, computational model. This method of analysis is effective but computationally expensive. It also requires complex calculation algorithms regarding nonlinear analysis. Of course, nonlinear analysis can be modeled using the Gauss integration points in the finite element. However, in relation to the rearrangement of the mesh discretization method, this solution can be less precise. This is caused by the necessity of approximation for results obtained in the Gauss points, e.g., transmission of strain values to the nodes of finite element.

In the performed analysis, the authors decided to use one-dimensional multi-area elements as girder columns, shown in Figure 3b. This is justified as follows:

• The girder column can be modeled using single one-dimensional double-hinge element, because it is proved by the practical aspects concerning erection of the concrete posttensioned girders; lack of concrete continuity in the area of column connections with girder chords was observed very often.
• In the procedure of the integration in sub-areas (described in the next section), the width of crack can be easily adjusted, and the localization of crack occurrence can be chosen as well. This second aspect is especially relevant in transient processes, where the location of damage is unknown in advance and depends on the state of structure response.
• There is an opportunity whereby elastic elements with significant differences of stiffness characteristics can be connected in series and described by a single finite element.

The adequacy of the proposed method has been checked exactly on the example of the girder column. For this purpose, the single column was separated and the tension force P was applied, as shown in Figure 4. The value of the displacement in form of extension of such a rod (bar element) was calculated. Calculations have been made using Maxima program [24], in two variants:

• according to standard stiffness summation—integration in sub-areas.
• according to stiffness summation—integration in sub-areas with χ parameter—a multi-area element.

In both variants the calculations have been performed using FEM and only one finite element has been applied. The stiffness matrix for the element has been obtained using the integration in sub-areas procedure. The main difference between these two variants is the distribution of shape functions inside the finite element, which has been explained in Section 3.

The calculation results have been presented in the

Table 1. Comparison of the displacement values [mm], Figure 5.

Applied Method	Width of the Crack ls
	ls = 1 mm	ls = 5 mm	ls = 10 mm	ls = 20 mm	ls = 50 mm
Reference solution—FEM three separate elements	0.754	0.815	0.891	1.042	1.498
Stiffness summation from sub-areas—FEM one finite element integrated in sub-areas	0.739	0.742	0.746	0.753	0.776
Calculation error of stiffness summation from sub-areas [%]	1.9	8.9	16.3	27.8	48.2
Multi-area elements—FEM one finite element integrated in sub-areas with modification of strain field	0.754	0.815	0.891	1.042	1.498
Calculation error of multi-area elements [%]	0.0	0.0	0.0	0.0	0.0

. As a variable parameter the crack width l_s was established in the range from 1 mm to 50 mm. As a reference solution, the values of displacement obtained from FEM calculations with discretization in form of three separate finite elements have been taken. The diagram presenting calculation error has been shown in Figure 5.

According to the presented results, it can be noticed that the calculations with the application of only a standard stiffness summation procedure are insufficient. Therefore, multi-area elements should be used. In the case of 1-D elements, the use of multi-area elements completely eliminated the calculation error.

3. General Theoretical Assumptions of One-Dimensional Multi-Area Elements

The conception of multi-area finite element has been presented and explained in Figure 6, on the example of one-dimensional element. Comparison between standard solution (Figure 6a) and multi-area elements (Figure 6b) shows that the deformation field in both cases is described by the shape functions [25,26]. However, in multi-area elements, the shape function distribution is modified by the χ parameter. This modification (χ) causes that stiffness is redistributed inside the finite element. Therefore, the basic assumption of multi-area conception is appropriate adjustment of shape function distribution to changes of stiffness characteristics in sub-areas, as shown in Figure 6b. Shape functions modified in this way could be described for example using spline functions. Theoretical considerations concerning spline functions are included in chapters 14 and 15 of the work [27]. The use of spline functions in FEM calculations is presented in chapter 18 of work [28].

However, the authors propose quite uncomplicated solution in form of modified stiffness summation from sub-areas, as it is described by Formula (1). The advantage of this solution is that there is no need for homogenization process.

The stiffness matrix of a single sub-area is determined from the Formula (2).

$$K_k\left(\chi \right)=\underset{v_k}{{\displaystyle \int }}{B}_k(\chi )·D_k·B_k\left(\chi \right)·d{v}_k$$

(2)

Global coordinates x and generalized dimensionless coordinates ξ are enclosed in the following intervals:

$$x\in \langle -\frac{L}{2},+\frac{L}{2}\rangle ,\xi \in \langle -1,+1\rangle , \text{where:} \xi =\frac{2}{L}x.$$

(3)

Shape functions expressed in natural coordinates ξ have the following forms for i and j node, respectively:

$$\begin{array}{c}{N}_k^i=\frac{{\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l}{{\xi }_{l+1}-{\xi }_l}·\left(\xi -{\xi }_l\right)+{\chi }_l·m_l\\ N_k^j=1-N_k^i\end{array}$$

(4)

The strain field:

$$\begin{gathered} \begin{array}{c}{B}_k=L·N_k=\left|\frac{d}{d\xi }\right|·\left|\begin{array}{cc}{N}_k^i& N_k^j\end{array}\right| \\ =\\ =\left|\begin{array}{cc}\frac{{\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l}{{\xi }_{l+1}-{\xi }_l}& -\frac{{\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l}{{\xi }_{l+1}-{\xi }_l}\end{array}\right|=\left|\begin{array}{cc}b& -b\end{array}\right|\end{array} \end{gathered}$$

(5)

After substituting to the Equation (2), final formula has been obtained:

$$K_k\left(\chi \right)=\underset{v_k}{{\displaystyle \int }}{B}_k(\chi )·D_k·B_k\left(\chi \right)·d{v}_k=\underset{{\xi }_l}{\overset{{\xi }_{l+1}}{{\displaystyle \int }}}\left|\begin{array}{c}b\\ -b\end{array}\right|·\left|E_k\right|·\left|\begin{array}{cc}b& -b\end{array}\right|·J·d\xi$$

(6)

where:

• B_k—strain matrix of sub-area k,
• D_k—material matrix, D_e^k = [E_k],
• N_k—shape function matrix of sub-area k,
• J—Jakobian matrix equals 2/L,
• L—length of the whole element.

After transformations, stiffness matrix of the component sub-area has the following explicit form:

$$\begin{gathered} K_k\left(\chi \right)={\left(EA\right)}_k·{\left(\frac{\left({\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l\right)}{\left({\xi }_{l+1}-{\xi }_l\right)}\right)}^2·\underset{{\xi }_l}{\overset{{\xi }_{l+1}}{{\displaystyle \int }}}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]·J·d\xi \\ =\left[\begin{array}{cc}\frac{2{\left(EA\right)}_k}{L}\frac{{\left({\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l\right)}^2}{\left({\xi }_{l+1}-{\xi }_l\right)}& -\frac{2{\left(EA\right)}_k}{L}\frac{{\left({\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l\right)}^2}{\left({\xi }_{l+1}-{\xi }_l\right)}\\ -\frac{2{\left(EA\right)}_k}{L}\frac{{\left({\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l\right)}^2}{\left({\xi }_{l+1}-{\xi }_l\right)}& \frac{2{\left(EA\right)}_k}{L}\frac{{\left({\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l\right)}^2}{\left({\xi }_{l+1}-{\xi }_l\right)}\end{array}\right] \end{gathered}$$

(7)

or the same stiffness matrix can be written in a different way:

$$K_k\left(\chi \right)=\frac{2{\left(EA\right)}_k}{L}\left[\begin{array}{cc}c& -c\\ -c& c\end{array}\right], \mathrm{where}:c=\frac{{\left({\chi }_{l+1}·m_{l+1}-{\chi }_l·m_l\right)}^2}{\left({\xi }_{l+1}-{\xi }_l\right)}$$

(8)

The essence of the presented multi-area elements is the use of modified shape functions described by proper parameters. These parameters are of course the dimension of the sub-area ξ and the distribution of the deformation field—parameter χ.

If the stiffness matrix of one-dimensional multi-area element is concerned, the strain field distribution is adapted to stiffness changes of sub-areas inside the finite element exactly by the use of parameter χ. Formulation of this parameter χ is based on the assumption that the sub-areas create a system of springs connected in series, as shown in Figure 7.

The value of adaptive shape function ${\chi }_l·m_l$ is calculated according following formula:

$${\chi }_l·m_l=\frac{\overline{k}}{\overline{k_l}}$$

(9)

where:

$\overline{k}$—resultant stiffness for total spring system (Figure 7):

$$\overline{k}=\frac{1}{{\sum }_{ii=1}^n\frac{1}{k_{ii}}}$$

(10)

and $\overline{k_l}$—resultant stiffness for the part of spring system (Figure 7) to the point where value $m_l$ is calculated:

$$\overline{k_l}=\frac{1}{{\sum }_{ii=l}^n\frac{1}{k_{ii}}}$$

(11)

In the Formulas (10) and (11), $k_{ii}$ denotes stiffness of single sub-area and index n is total number of sub-areas in the finite element. In Figure 7, the stiffness of single sub-area is marked respectively as: k₁, k₂, k₃ and the total number of sub-areas is n = 3.

4. Directions of Potential Application

Calculation Test—Subsidence of the Ground Subsoil

Apart from presented in the Section 2, previous example of concrete girder KBOS-24 column calculation, the developed method of one-dimensional special finite elements has been implemented in the calculation of ground subsoil subsidence. This concerns mainly the calculation of the multi-layered subsoil, where each layer is characterized by different thickness and strain modulus [29]. One of the examples where subsidence values have been calculated is the analysis of a layered subsoil under a building wall. In the presented case the calculations are limited to one-dimensional elastic model of subsoil [30].

In the computational test, the subsidence values for the multi-layered subsoil under the wall of a 5-storey building were compared, as shown in Figure 8. For calculation a reinforced concrete wall measuring 9.6 m × 14.0 m × 0.15 m, loaded by the floor ceilings and dead load has been taken. Calculations were made in three variants, using different calculation methods.

In the variant I the calculations were carried out according to hypothesis of the elastic subsoil parameters (Winkler’s hypothesis) and the formulas given in [30]. The problem of modeling the multi-layered subsoil using the hypothesis of the elastic subsoil is to determine the stiffness of the layer (as a spring). The solution is to use exactly the Winkler’s hypothesis, according to which the subsidence of the elastic subsoil s is proportional to the acting load q, as in the following Expression (12):

$$q= k_z·s$$

(12)

The value of flexibility coefficient k_z for homogeneous soil to the depth z [30] can be determined according to the formula:

$$k_z=\frac{E_0}{\omega ·B·\left(1-{\nu }^2\right)}$$

(13)

On the other hand, in case of the multi-layered subsoil, the coefficient k_z is the sum of the flexibility coefficients for the individual layers of the subsoil. For a single layer i, this coefficient is determined from the following expression:

$$\begin{gathered} k_z^i=\frac{E_0^i}{\mathsf{\Delta }{\omega }_i·B·\left(1-{\nu }^2\right)} \\ \mathrm{and}\mathsf{\Delta }{\omega }_i={\omega }_i-{\omega }_{i-1}, \end{gathered}$$

(14)

where:

• q—acting load,
• E₀—strain modulus of the subsoil,
• B—the width of the loaded area (foundation),
• ν—coefficient of the lateral expansion of the ground,
• ω_i (Δω_i)—influence coefficient, depending on the shape of the loaded area (foundation) determined according to the corresponding tables and nomograms [30].

In the variant II, calculations have been made using a two-dimensional FEM model in order to compare and control the results. For plane elements, that were implemented in ground subsoil modelling, a transformation of the strain modulus has been applied in order to use the flexibility of the subsoil as in Winkler’s hypothesis:

$$E^i=E_0^i·\frac{h_i}{\mathsf{\Delta }{\omega }_i}$$

(15)

where:

• $E_0^i$—strain modulus for particular layer of the ground,
• h_i—the thickness of the particular layer,
• Δω_i—influence coefficient, depending on the shape of the loaded area (foundation) determined according to the corresponding tables and nomograms [30].

In the variant III, calculations have been performed using special one-dimensional finite elements integrated in sub-areas, presented in the paper and implemented in the authors’ analysis system “ORCAN” [23].

As a results of the calculations, the following subsidence values have been obtained.

Variant I: Winkler’s hypothesis:

$$\begin{gathered} s_1= \frac{q·\mathsf{\Delta }{\omega }_1·B·\left(1-{\nu }^2\right)}{E_0^1}= \frac{189.19 \mathrm{kN}/\mathrm{m} ·0.86·1.0\mathrm{m}·\left(1-{0.3}^2\right)}{80,000 \mathrm{kPa}}=1.851·{10}^{-3} \mathrm{m} \\ {s}_2= \frac{q·\mathsf{\Delta }{\omega }_2·B·\left(1-{\nu }^2\right)}{E_0^2}= \frac{189.19 \mathrm{kN}/\mathrm{m} ·0.337·1.0\mathrm{m}·\left(1-{0.3}^2\right)}{20,000 \mathrm{kPa}}=2.901·{10}^{-3} \mathrm{m} \\ {s}_3= \frac{q·\mathsf{\Delta }{\omega }_3·B·\left(1-{\nu }^2\right)}{E_0^3}= \frac{189.19 \mathrm{kN}/\mathrm{m} ·0.259·1.0\mathrm{m}·\left(1-{0.3}^2\right)}{130,000 \mathrm{kPa}}=0.343·{10}^{-3} \mathrm{m} \\ \mathrm{Total}\mathrm{subsidence}:s=s_1+s_2+s_3=5.095·{10}^{-3} \mathrm{m}=5.095\mathrm{mm} \end{gathered}$$

Variant II: FEM model for the entire ground cross-section with the use of plane two-dimensional elements. The following values of the strain modulus respectively for first, second and third layer (Formula (15)) has been applied:

$$E^1=E_0^1·\frac{h_1}{\mathsf{\Delta }{\omega }_1}=80,000 \mathrm{kPa}·\frac{2.0 \mathrm{m}}{0.86}=186046.51 \mathrm{kPa}$$

$$E^2=E_0^2·\frac{h_2}{\mathsf{\Delta }{\omega }_2}=20,000 \mathrm{kPa}·\frac{2.0 \mathrm{m}}{0.337}=118694.36 \mathrm{kPa}$$

$$E^3=E_0^3·\frac{h_3}{\mathsf{\Delta }{\omega }_3}=130,000 \mathrm{kPa}·\frac{2.0 \mathrm{m}}{0.259}=1003861.00 \mathrm{kPa}$$

$$\mathrm{Total}\mathrm{subsidence}:s=5.429·{10}^{-3} \mathrm{m}=5.429 \mathrm{mm}$$

$$\begin{gathered} \mathrm{Total}\mathrm{subsidence}\mathrm{after}\mathrm{blocking}\mathrm{of}\mathrm{the}\mathrm{ground}\mathrm{lateral}\mathrm{expansion}\mathrm{possibility}: \\ s=5.116·{10}^{-3} \mathrm{m}=5.116 \mathrm{mm} \end{gathered}$$

Variant III: Special one-dimensional multi-area elements (replacement of the entire ground cross section by single springs). After the discretization, along the foundation 50 nodes were established, so the length assigned to single spring is equal: $11.6 \mathrm{m}÷50 springs=0.232 \mathrm{m}$.

The flexibility coefficients for the individual layers of the subsoil are following:

$$\begin{gathered} k_s^1= \frac{E_0^1}{\mathsf{\Delta }{\omega }_1·B·\left(1-{\nu }^2\right)}·A= \frac{80,000 \mathrm{kPa}}{0.86·1 \mathrm{m}·\left(1-{0.3}^2\right)}·\left(1 \mathrm{m}·0.232 \mathrm{m}\right)=23715.82\frac{\mathrm{kN}}{\mathrm{m}} \\ {k}_s^2= \frac{E_0^2}{\mathsf{\Delta }{\omega }_2·B·\left(1-{\nu }^2\right)}·A= \frac{80,000 \mathrm{kPa}}{0.337·1 \mathrm{m}·\left(1-{0.3}^2\right)}·\left(1 \mathrm{m}·0.232 \mathrm{m}\right)=15130.27\frac{\mathrm{kN}}{\mathrm{m}} \\ {k}_s^3= \frac{E_0^3}{\mathsf{\Delta }{\omega }_3·B·\left(1-{\nu }^2\right)}·A= \frac{80,000 \mathrm{kPa}}{0.259·1 \mathrm{m}·\left(1-{0.3}^2\right)}·\left(1 \mathrm{m}·0.232 \mathrm{m}\right)=127964.70\frac{\mathrm{kN}}{\mathrm{m}} \end{gathered}$$

According to flexibility coefficients assigned to each ground layer, resultant stiffness for the single special one-dimensional multi-area element (single spring) was calculated:

$$\begin{gathered} {\overline{k}}_s= \frac{1}{\frac{1}{k_s^1}+ \frac{1}{k_s^2}+ \frac{1}{k_s^3}}=\frac{1}{\frac{1}{23715.82 \mathrm{kN}/\mathrm{m}}+ \frac{1}{15130.27 \mathrm{kN}/\mathrm{m}}+ \frac{1}{127964.70 \mathrm{kN}/\mathrm{m}}}=8615.25\frac{\mathrm{kN}}{\mathrm{m}} \\ \mathrm{Total}\mathrm{subsidence}:s=5.156·{10}^{-3}\mathrm{m}=5.156 \mathrm{mm} \end{gathered}$$

Comparison of the values of subsoil subsidence obtained from three, different calculation variants has been presented in the

Table 2. Comparison of the subsoil subsidence values [mm], Figure 8

Calculation Variant	Value of Subsoil Subsidence s
Variant I—Winkler’s hypothesis	$0.005095 \mathrm{m}=5.095 \mathrm{mm}$
Variant II—FEM model for the entire ground cross-section with the use of plane two-dimensional elements.	$0.005429 \mathrm{m}=5.429 \mathrm{mm}$
	$0.005116 \mathrm{m}=5.116 \mathrm{mm}$1
Variant III—special one-dimensional multi-area elements (replacement of the entire ground cross section by single springs)	$0.005156 \mathrm{m}=5.156 \mathrm{mm}$

¹ After blocking of the ground lateral expansion possibility.

.

5. Discussion

The aim of this paper was to present a new type of finite element, which allow to perform calculations of heterogeneous, multi-coherent structures using less number of finite elements and provide proper accuracy of the results.

The idea and the main assumptions of the special multi-area finite elements have been explained on the example of one-dimensional (linear) element. The results of the computational tests have been presented in order to confirm the suitability of these elements in calculations of structures containing areas with different physical and geometrical parameters.

The computational calculations and results presented in the Table 1, confirm that the use of elements integrated in sub-areas in a standard way (without strain field modification inside the element) can lead to incorrect results. The greater difference between the subdivision stiffness, the greater computational error is. Therefore, the displacement field modification is necessary. Using multi-area special elements it is possible to obtain the reference solution even for significant differences in the stiffness of the component sub-areas. Reference solution means the results obtained from FEM calculations with fine mesh discretization where each sub-area is described by separate finite element.

According to the presented examples and obtained results (Table 1 and Table 2), it can be noticed that the application of special one-dimensional (linear) finite elements in the calculations of structures with variable stiffness and geometric parameters, provides results that are consistent with reference solution or are very close to the expected solution. Simultaneously the number of finite elements used in calculations has been decreased.

The conception is based on the assumption that area of the structure with different stiffness and geometrical parameters is described by a single element integrated in subdivisions (sub-areas). However, the main idea of the special multi-area elements, that makes them different from the typical elements integrated in sub-areas with a linear distribution of the shape functions, is the application of the modification of the strain field distribution within the finite element, depending on the changes of geometrical and material parameters. For this purpose the original functional for modification of the linear distribution of the shape functions has been developed. The stiffness matrix for the multi-area elements has been determined in the explicit form. This increases the computational efficiency compared to the numerical integration.

6. Conclusions

The research presented in the paper allows to define the following conclusions:

• An original conception, which has been used in FEM calculations, for the analysis of multi-coherent structures consisted of areas with different geometrical and material properties was developed and presented.
• According to the developed method, a special stiffness matrix for one-dimensional multi-area finite element has been derived. As the finite element consists of parts with different stiffness parameters, in the formulation of the stiffness matrix, the integration in sub-areas procedure has been applied.
• The conducted research shows that integration in sub-areas is insufficient when linear continuous distribution of shape functions is assumed. The necessary is modification of shape functions distribution. The authors developed this modification in the form of adaptive shape functions. According this approach, the special multi-area elements have been formulated.
• The performed calculations confirm that the use of multi-area elements provides results that are consistent or are very close to the expected solution.
• The use of special multi-area elements helps to reduce the number of finite elements and thereby to reduce the computational calculation time.
• The obtained proper results in the case of one-dimensional finite elements provide the opportunity to continue research on 2D and 3D elements.