A coupled thermo-hydro-mechanical finite element formulation for curved beams in two-dimensions

To enable the use of beam elements in the modelling of coupled thermo-hydro-mechanical (THM) geotechnical problems, a fully coupled and robust THM formulation is required. This paper presents such a formulation which allows both fluid flow and heat transfer along a 2D curved beam, while ensuring compatibility with coupled THM solid elements commonly used to discretise soils. Verification exercises and application with the proposed coupled beam element are carried out to demonstrate its satisfactory behaviour. The results of these analyses are compared against closed form solutions, solutions obtained using solid elements, and field measurements, showing an excellent agreement.


Introduction
Coupled thermo-hydro-mechanical (THM) interactions in porous media (e.g. soils) have been the focus of recent research as many geotechnical engineering problems, such as geothermal energy systems or storage of high-level radioactive waste, involve temperature variations [1]. To adequately simulate the coupled THM behaviour of porous media in a finite element (FE) analysis, a number of coupled THM formulations have been developed in the literature [2][3][4][5][6][7][8][9][10][11][12]. However, most of the existing formulations have only been applied to solid elements, which are commonly used to discretise the porous medium. In theory, it is possible to use these coupled THM solid elements to model structural components such as thermo-active retaining walls, piles, tunnel linings or heat exchanger pipes. In practice, however, this may result in either a very large number of elements or elements with unacceptably large aspect ratios [13]. To overcome these shortcomings and adequately model the interaction between the porous medium and the structure in a coupled THM FE analysis, a special type of beam element, which can accommodate both consolidation and heat transfer, is required. Additionally, beam elements are frequently a preferred element type for discretising structural components, as they are formulated directly in terms of structural forces, rather than the component stresses, hence being easier to use.
Russell [14] adopted a three-noded line element, described as a drainage element, to model a vertical drain in an isothermal consolidation FE analysis. However, as this coupled hydro-mechanical (HM) line element has a zero cross sectional area, it becomes problematic when applying it to model the mechanical behaviour of structural components. Moreover, the off-diagonal coupling terms in the HM formulation of the line element are zero, which prevents using the line element to simulate the fully coupled HM behaviour. To the authors' knowledge, coupled THM formulations for beam elements have not been published elsewhere.
The mechanical formulation for a two-dimensional (2D) curved beam element developed by Day and Potts [15] is used here as a basis from which to extend the formulation to account for the effects of both pore fluid pressure and temperature changes, resulting in a coupled THM formulation. The authors have recently developed a new fully coupled THM formulation for 2D and 3D solid elements [12]. In order to ensure compatibility between solid elements and beam elements in an FE analysis, the new THM formulation for 2D beam elements adopts the same development approach and assumptions. Therefore, the proposed beam elements can be used to model the behaviour of either a two-phase porous medium intrinsically, or a single-phase material (e.g. concrete structures) by setting appropriate material properties to each phase. Moreover, as the formulation of each coupling system is developed individually, it is possible to disable any of the three systems (thermal, hydraulic and mechanical) if they are not active in the analysis, thus enabling the coupled THM beam element to be used in coupled HM, TM or TH, or uncoupled mechanical, hydraulic or thermal analyses. This flexible formulation for the coupled THM 2D curved beam element has been successfully implemented into the bespoke software ICFEP [13], to demonstrate the element's performance in a number of verification exercises, as well as case studies. In the present formulation, displacement, pore fluid pressure and temperature are the nodal degrees of freedom.

. Strain definition
The beam element developed by Day and Potts [15] was designed as a special type of zero-thickness line element for modelling structural components, e.g. retaining walls and tunnel linings, in a 2D analysis. It is formulated directly in terms of bending moments, axial and shear forces and their associated strains, under the consideration that the dimension of the structural component in the in-plane transversal direction is zero. Fig. 1 shows a general curved 2D beam element in the global coordinate system x-y. The corresponding degrees of freedom are marked as displacements u and v in the x and y directions respectively, and the cross section rotation θ.
The strains of the 2D curved beam element have to be defined in a local coordinate system, with directions tangential and normal to the beam [15]: Axial (longitudinal) strain Bending strain Shear strain where l is the distance along the beam, R is the radius of curvature, u l and v l are the displacements tangential and normal to the beam, respectively. The angle α b in Fig. 1 is the rotation from global to local coordinate directions. For axisymmetric analysis, the definition of additional terms is required: Circumferential membrane strain Circumferential bending strain where r 0 is the circumferential radius.

Mechanical equilibrium
When only the mechanical behaviour of the beam is taken into account, the constitutive relation between the strain terms and the element forces and bending moments is defined in the same manner as that in [15]: Ψ represents the incremental forces and moments vector. Note that, since the beam is represented as a line (zero-thickness) in the formulation of this type of element, no strains or forces related to the in-plane transverse direction (i.e. perpendicular to the line representing the beam) are defined in these vectors. It should be noted that his assumption has been widely used in the literature (e.g. [16]). In axisymmetric analysis, the force/moment components are defined as: F Δ a is the incremental axial force, M Δ a is the incremental bending moment, S Δ is the incremental shear force, F Δ Ψ is the incremental circumferential force and M Δ Ψ is the incremental circumferential bending moment. For plane strain analysis, F Δ a , S Δ and M Δ a are the incremental in-plane axial force, shear force and bending moment respectively, F Δ Ψ and M Δ Ψ are the incremental out of plane force and bending moment respectively, and [ ] is the constitutive matrix of the beam element, the form of which depends on the assumed material behaviour (e.g. linear elastic, non-linear elastic, elasto-plastic). For an isotropic linear elastic material in plane strain and axisymmetric analysis, D [ ] takes the form: where A is the cross sectional area of the beam, I is the second moment of area, E is the Young's modulus, μ is the Possion's ratio, G is the shear modulus and k is a shear correction factor which is dependent on the shape of the cross section. The values of I and A are specified in plane strain and axisymmetric analysis for a unit depth in the out of plane direction or unit angle, respectively. For a coupled problem involving both fluid and heat transfer, as the effects of both pore fluid pressure and temperature changes are isotropic, they are assumed to affect only the direct strains of the beam element (i.e. ε Δ l and ε Δ Ψ ). Naturally, since in the in-plane transverse direction the dimension of this element is assumed to be zero, no direct strain in this direction is evaluated. This is consistent with the mechanical behaviour as noted above. Moreover, the temperature of the solid phase is assumed to be the same as that of the fluid phase, which implies instantaneous equilibrium in temperature between the two phases. Therefore, the incremental total strain ε {Δ } can be written as the sum of the incremental strain due to stress changes (mechanical strain), ε {Δ } σ , and the incremental strain due to temperature changes (thermal strain), ε {Δ } T : By combining Eqs. (6) and (8) and adopting the principle of effective stress, the force equilibrium equation for the coupled THM beam element can be written as: The finite element equations related to the force equilibrium of the coupled THM beam element can be obtained by applying the principle of minimum potential energy as: where B [ ] is a matrix which contains derivatives of the shape functions based on local to global coordinate transformation [13] {Δ } nG are the global vectors of incremental nodal displacement, pore fluid pressure and temperature respectively, which in the FE analysis are the primary unknowns (degrees of freedom), F {Δ } is the vector of applied external forces and the subscripts G and nG denote the global matrix or vector and the global nodal values respectively. It should be noted that the form of Eq. (10) is similar to that for coupled THM solid elements [12], while the matrix K [ ] G , which represents the mechanical constitutive behaviour, is the same as that proposed by [15]. Additionally, L [ ] G is the HM coupling term accounting for the effect of the pore fluid on the mechanical response, M [ ] G is the TM coupling term accounting for the effect of temperature on the mechanical response and R {Δ } G denotes any externally applied forces and/ or moments.

Fluid flow along a curved beam
Similar to the formulation for saturated soils outlined in [12], compressible fluid flow along a beam element under isothermal conditions is assumed to be governed by the continuity equation, which can be written as: where v f l , is the velocity of the fluid flowing along the beam element, n is the porosity of the beam skeleton, K f is the bulk modulus of the fluid, Q f represents any pore fluid sources and/or sinks, ε v is the volumetric strain of the beam which can be expressed as = + ε ε ε v l Ψ as the dimension of the structural component in the transversal direction is not taken into account, and t is time.
Under non-isothermal conditions, the volume of fluid flowing into or out of the beam element induced by the difference in the thermal expansion coefficients of the solid and the fluid phases also needs to be considered. In this case, Eq. (15) can be rewritten as: where α T f , is the linear thermal expansion coefficient of the fluid. The first term on the left-hand side of Eq. (16) represents the flow of fluid along the beam element, the second term denotes the changes in the volume of the fluid due to its compressibility, while the third term relates to the excess fluid flow generated by the difference between the thermal expansion coefficients of the two phases. The term on the righthand side of Eq. (16) represents the change in the mechanical volumetric strain, which reflects the changes in the pore space available for the fluid within the porous medium to flow. Darcy's law defines the velocity of the fluid flow along a beam element as: where k f l , is the hydraulic conductivity of the beam, and h f l , is the hydraulic head defined as: with γ f being the bulk unit weight of the fluid and the vector T being the global unit vector parallel, but in the opposite direction, to gravity. Substituting Eqs. (17) and (18) into Eq. (16) and applying the principle of virtual work, as well as the divergence theorem, yields: Adopting the isoparametric formulation for beam elements, similar to that for solid elements [13], the coordinates and global degrees of freedom at any point on the beam element can be expressed in terms of the nodal values using the shape functions, as: where s is the natural ordinate that varies from −1 to +1 over the element length (see Fig. 2), the subscript i indicates the elemental node number, and the prime denotes the derivative with respect to s. It should be noted that the above formulations also applies to a 2-noded beam element which only contains the end nodes (nodes 1 and 2 in Fig. 2). The transformation of the hydraulic equation from local to global coordinates can be performed with the help of Fig. 3 as: where the determinant of J is calculated as: Noting that (see Fig. 3): Eq. (19) can then be rewritten as: Approximating dp dt / Coordinates of a three-noded isoparametric beam element.
G represent the effect of stress-strain and temperature changes, respectively, on the hydraulic behaviour, and G and [S ] G are the hydraulic conductivity and compressibility matrices of the pore fluid, respectively. When a time marching scheme is adopted, such that for the time interval t t [ , ] 1 , where t 1 is the initial time and t is the current time, the following is obtained: where β 1 is the time marching parameter for pore fluid flow, and p ({ } ) f nG 1 is the initial pore fluid pressure at the beginning of the given time interval. Finally, Eq. (29) becomes:

Heat transfer along a curved beam
Similar to fluid flow, heat is assumed to transfer along the beam element. Adopting the law of conservation of energy results in: where Q T is the total heat flux per unit volume, Q T represents any heat source/sink, dV is an infinitesimal volume of the beam and Φ T is the heat content of the beam per unit volume which, for a fully saturated porous medium, can be calculated as: with n being the porosity of the beam skeleton, ρ f and ρ s being the densities of the pore fluid and beam skeleton, respectively, C pf and C ps being the specific heat capacities of the pore fluid and beam skeleton, respectively, and T r being a reference temperature. Ignoring heat radiation, the total heat flux, Q T , can be written as the sum of heat conduction, q d , and heat advection, q a , which are expressed respectively as: where k Tl is the thermal conductivity of the beam. Substituting Eqs. (38)-(40) into Eq. (37) results in: Noting that = + n e e /(1 ), and = + dV ed V (1 ) s , where dV s is the infinitesimal volume of the beam skeleton, Eq. (41) can be rewritten as: Under isothermal conditions, dV s is generally assumed to be constant, regardless of the changes in effective stresses. Under non-isothermal conditions, however, dV s is temperature dependent and can be written as: where dV s0 is the initial volume of the beam skeleton, which is assumed to be constant for a coupled THM problem, and ε vT is the thermal volumetric strain. Substituting Eq. (43) into Eq. (42), eliminating dV s0 , and applying mass conservation of the material, yields: It should be noted that the effect of the mechanical volume change on heat transfer, i.e. effect of a variation in void space, ∂ ∂ e t / , on heat storage, is retained in Eq. (44), while this term is often ignored in the literature when dealing with solid elements (e.g. [10]). Adopting the assumption that the change in void ratio is only a result of the mechanical strain [12], the following relationship in a small displacement analysis can be established: Substituting Eq. (45) into Eq. (44), and applying the coordinate transformation process as for the hydraulic equations, the finite element formulation governing heat transfer along the beam element is expressed as: Adopting a time marching scheme which is similar to Eq. (35) results in: where β 2 and α 1 are time marching parameters for heat transfer, and T ({ } ) nG 1 is the initial temperature at the given time interval.

Finite element formulation for a fully coupled THM problem
Assembling Eqs. (10), (36) and (52) leads to the finite element formulation for a fully coupled THM beam element:

Verification and application of the coupled THM 2D curved beam element
The formulation presented above for a coupled THM beam element has been implemented into the bespoke FE code ICFEP. To demonstrate the new capabilities of the coupled beam element, a series of benchmark analyses (i.e. T, TH, HM and THM analyses) are carried out. The results of these exercises are compared against either existing analytical solutions, or solutions obtained using solid elements. It should be noted that the validity of the THM formulation for solid elements has been demonstrated in [12] by comparisons against both closed form solutions and numerical studies found in the literature. In addition, an application involving a coupled HM problem was carried out where the numerical solutions using the proposed beam elements are compared to both field measurement and numerical results found in the literature where a HM coupled line element was employed.

Conductive heat transfer along a curved beam (T)
Both fluid and heat are assumed to transfer along the beam element. As noted earlier, to model this behaviour, local axes in the tangential and normal directions to the beam are adopted in the governing equations. To verify the associated transformation between the local and global coordinate systems, an exercise of heat transfer along a curved beam is carried out. An exercise involving fluid flow has also been performed following a similar procedure, but is not included here for brevity.
A plane strain analysis of heat conduction along a 0.1 m long curved beam is carried out using 2-noded beam elements with an element length in the axial direction of 0.01 m (i.e. 10 elements). The geometry of the problem is shown in Fig. 4(a) and the material properties are listed in Table 1. The initial temperature is 20°C and constant temperatures of 40°C and 20°C are prescribed at points 1 and 2, respectively. The time-step size is chosen arbitrarily as 60 s with a time marching parameter of = α 0.8 1 . A value of 1.0 m 2 was adopted for the area of the beam in the analysis for simplicity.
To verify the results using beam elements, the same problem is modelled with 4-noded solid elements which are also assigned the material properties shown in Table 1. The geometry of the dashed curved line in Fig. 4(b), which lies in the middle of the curved plane generated by the 2D solid elements, is the same as that of the curved beam shown in Fig. 4(a). Constant temperatures of 40°C and 20°C are prescribed at the sides between points 1 and 2 and between points 3 and 4, respectively, while a no heat flux boundary condition is prescribed at the remaining boundaries. Fig. 5 compares the results of the temperature distribution obtained using 2-noded beam elements and 4-noded solid elements at a transient state (t = 600 s) and the steady state, showing a perfect agreement between the two analyses.

Plane strain analysis
A plane strain analysis of a one-dimensional (1D) conductive-advective flow along a 1 m bar of soil with element length of 0.01 m was conducted to verify the thermo-hydraulic coupling of the proposed beam element. The geometry of the considered problem is shown in Fig. 6 and the material properties used in the simulations are listed in Table 2. To include the conductive-advective heat transfer, constant pore water pressures of 10 kPa and 0 kPa were prescribed at the left side  (i.e. node 1) and right side (i.e. node 2) of the bar, respectively, resulting in a constant velocity of 1.02 × 10 −5 m/s from left to right. A constant temperature boundary condition of 10°C was applied at node 1, while the initial temperature was 0°C, generating heat transfer along the bar. A coupled thermo-hydraulic boundary condition, which is detailed in [18], was prescribed at node 2, to balance the change of energy associated with the water flow through the boundary. The time-step was chosen arbitrarily as 60 s, larger than the minimum value established by [17] to avoid thermal-shock in this type of problem. Fig. 7 compares the results of the temperature distributions along the bar obtained at different time instants (t = 3600 s, 12,000 s and 24,000 s) in the analysis against those obtained using the analytical solutions from van Genuchten and Alves [19] where the same boundary conditions were prescribed, showing very good agreement.

Axisymmetric analysis
To verify the behaviour of the proposed beam elements in modelling fluid flow and heat transfer in an axisymmetric analysis, conductiveadvective flow along a horizontal disc as shown in Fig. 8, was simulated. The same material properties as those listed in Table 2 were adopted and the same initial conditions and boundary conditions as

Table 2
Material properties for simulation of conductive-advective flow.     those prescribed in the plane strain analysis were employed to generate the conductive-advective flow, though in this case in the opposite direction (i.e. from right to left). The fluid velocity increases along the inward radial direction with a maximum at the left side due to the axisymmetric nature of the problem. This problem was also modelled with solid elements adopting the same thermal and hydraulic properties as those listed in Table 2. To ensure the fluid flow is only in the radial direction with a zero vertical component due to the gravity effect, a hydrostatic initial pore fluid pressure condition was specified for the solid elements, while an increase of 10 kPa and a zero change of pore fluid pressure were prescribed at point 2 and point 1 respectively. Fig. 9 compares the results of the fluid velocity distribution obtained using beam elements and solid elements, while Fig. 10 compares the temperature distributions at a transient state (t = 20,000 s) and the steady state, showing a perfect agreement between those two analyses.

1D elastic isothermal consolidation (HM)
To validate the hydro-mechanical coupling of the beam element, a coupled consolidation analysis is carried out. The exercise is based on the numerical tests firstly performed by [2,3], and subsequently used to validate other finite element software (e.g. [4,5,8]). While the original test was performed with 2D solid elements, the same problem is represented here with the mesh discretised with the newly developed coupled beam elements, as shown in Fig. 11. The original material properties used in the literature for solid elements are listed in Table 3. Clearly, to impose 1D conditions when using solid elements in a plane strain analysis, horizontal displacements need to be restrained along the lateral boundaries of the mesh. Such mechanical restriction will result in additional axial strains due to the Poisson's ratio of the material not being zero. Conversely, when using beam elements, these lateral restraints are not required since the element is considered to have zero thickness and, hence, no stresses or strains are evaluated in the in-plane transversal direction. As a result, if identical results are to be obtained, the axial stiffness of the beam elements needs to be adjusted to compensate for the absence of the effect of the Poisson's ratio in the in-plane transversal direction when these elements are used. This can be achieved by considering that, for a 1D elastic consolidation problem, the compressibility coefficient defined by the analytical solution of Biot [20] can be written as: Eq. (54) is valid for the analysis with solid elements, however, for the analysis with beam elements, it is necessary to modify the expression to account for the difference in the way the elements are formulated, yielding: Equating Eqs. (54) and (55) leads to the equivalent material properties for beam elements in Fig. 11, which are summarised in Table 3, assuming that the Poisson's ratio remains unchanged. Clearly, if the material would be characterised by a Poisson's ratio of zero, no adjustment of the Young's modulus would be required. The displacement boundary conditions are prescribed such that only 1D deformation is allowed. A downward vertical load of 2 N (equivalent to the stress of 1 Pa used in the original analysis) is applied at point 2. Additionally, the prescribed hydraulic boundary conditions ensure that pore fluid can only leave the beam through the top boundary (i.e. a zero change of pore fluid pressure boundary condition is prescribed at point 2 and a zero flow boundary condition is prescribed at point 1). The initial pore water pressure is hydrostatic between points 2 and 1. The time-step size was chosen as 1 s for the first 10 increments and 10 s for the remainder of the analysis, as suggested by [12]. The predicted vertical displacement at point 2 is compared to the analytical solution of Biot [20] in Fig. 12, showing an excellent agreement between the numerical and analytical results.

Application to the modelling of vertical drains
To further demonstrate the capability and behaviour of the presented coupled beam element, a coupled analysis of a case study of the Porto Tolle experimental embankment in Italy [21] is carried out here. The embankment, which had a crest width of 30 m and a length of over 300 m, was constructed at a constant rate of height increase to a  Table 3 Material properties for isothermal 1D consolidation analysis.
Solid element [18] Beam elements maximum height of 5.5 m (equivalent to a surface load of 99 kN/m 2 ) over a period of 4 months. Then, the embankment was left to consolidate for 10 months before it was removed. Prefabricated drains were installed to accelerate the consolidation of the 20-23 m thick layer of normally consolidated soft silty clay beneath the embankment (see Fig. 13). The case study was adopted by [22] to demonstrate the performance of procedures (i.e. geometry or permeability matching methods) used to model the effect of vertical drains in plane strain analysis on consolidation of the foundation soil beneath embankments.
To investigate the effect of the vertical drains beneath the Porto Tolle embankment, a plane strain analysis is performed with a mesh of a half-width unit cell (see Fig. 14), which has been shown in the literature (e.g. [14,21;23]) as appropriate for predicting vertical movements and pore pressures. In the unit cell, only the layer of soft silty clay is discretised as in [22]. The limitation of the unit cell analysis, compared to the analysis using a complete mesh of the embankment, is that the lateral movements cannot be predicted. The vertical drain is modelled with the coupled beam element (the line between points 1 and 4) using its coupled HM capability, while the surrounding normally consolidated soft clay is represented by coupled HM solid elements. The input initial stress profiles are plotted in Fig. 13 and material properties for both beam and solid elements are listed in Table 4. A hydrostatic initial pore fluid pressure condition along the depth was adopted. The properties for the solid elements have been shown by [22] to match the soil behaviour well, as demonstrated from in-situ and laboratory testing. The adopted constitutive model in both studies is the Modified Cam-clay model. It should be noted that the smear zone is neglected and the permeability matching scheme proposed by [24] is adopted here, resulting in a spacing of 4.0 m (2.0 m for half unit cell) and a horizontal hydraulic conductivity of 6.91 × 10 −5 m/day for the clay (a horizontal hydraulic conductivity of 3.54 × 10 −4 m/day was measured in the in-situ test). With respect to the hydraulic response of the beam element, which is characterised by a cross sectional area of 0.003 m 2 , a hydraulic conductivity of 254.1 m/ day was adopted in order to achieve the overall discharge capacity of 140 m 3 /year estimated by [25]. In this calculation, an in-situ hydraulic gradient of 0.5 was assumed, as suggested by [26].

Medium dense sand
Soft silty clay Dense silty sand  To simulate the embankment construction, the equivalent construction surcharge of 99 kN/m 2 was applied on the top boundary of the mesh over 50 increments, with a time-step size of 2.4 days (120 days in total). Subsequently, the consolidation process was modelled employing a further 50 increments with a time-step of 6.0 days (300 days in total). Displacement boundary conditions were prescribed so that only vertical displacement was allowed along the lateral boundaries (see Fig. 14). As sand layers existed both above and below the clay stratum, a no change of pore pressure boundary condition is prescribed at both the top and bottom boundaries of the unit cell, to ensure that drainage is possible. 8-noded solid elements and 3-noded beam elements were used in the analysis and each node in the mesh had both displacement and pore pressure degrees of freedom. [22] used the computer program CRISP in their analysis, in which the vertical drain was represented by a three-noded line element (i.e. so called "drainage" element in CRISP). As such, this element, as explained by [14], does not have a real cross sectional area and therefore the off-diagonal matrices L [ ] G in Eq. (53) become zero in a 2D plane strain consolidation analysis. Figs. 15-17 compare the results from this study against field measurements and those from [22]. It is evident that both numerical predictions of excess pore pressures do not match exactly the measured changes in pore pressures. The likely reason for this is that the analyses consider only the so called unit cell and not the realistic extent of the embankment and foundation soil domains. However, it is noted that the predictions of both the average surface settlement and the peak excess pore pressure are more satisfactory with the coupled beam element used in this study compared to those using the drainage element in the analysis by [22]. The differences between the predictions may be caused by two aspects: (1) The FE formulations are different between the coupled beam and drainage elements. While the formulation of the former is consistent with that of the solid elements in this study, the cross-coupling term was missing in the formulation of the latter. (2) The employed critical state model was different between ICFEP and CRISP. A Mohr-Coulomb hexagon was adopted for the shape of the yield surface in the deviatoric plane of the modified Cam-clay model in this study, while a circular shape was used in CRISP [27].
The dissipation of excess pore pressures measured in the field was slower than that computed in either of the numerical analyses, which [22] attributed to the presence of the organic gas in the soil beneath the embankment. This was not simulated in any of the analyses.

1D non-isothermal consolidation (THM)
A 1D thermo-elastic consolidation (i.e. fully coupled THM) analysis of a fully saturated porous medium, which was performed in the literature with solid elements (e.g. [2,4,5;8]), is also carried out here with beam elements. Cui et al. [12] verified the coupled THM formulation for solid elements in ICFEP by performing the same exercise, demonstrating a very good agreement between ICFEP predictions and the results reported in the literature. Therefore, the coupled THM formulation for beam elements in ICFEP can be validated against the existing solid elements in ICFEP.
The same mesh as that shown in Fig. 11 is used and the same material properties (see Table 5) as those in the original analysis by [2] are adopted. However, due to the difference in the definition of both mechanical and thermal strains, Poisson's ratio of both the solid and the beam elements is set to zero for simplicity in order for the comparison not to be affected by this parameter. It should be noted that this also leads to the same Young's modulus being adopted for solid and beam    elements. Similar to the previously described isothermal consolidation test, the mesh in this exercise is restrained to allow only 1D deformation. At point 2, a point load of 2 N, as well as no change in pore fluid pressure, is prescribed. In order to include the thermal effects, a constant temperature of 50°C and 0°C are applied at point 2 and point 1 respectively with an initial temperature of 0°C in the mesh. Fig. 18 compares the surface settlement (i.e. vertical movement of point 2) versus time, which is a consequence of variation in both pore pressure and temperature, computed using beam and solid elements. A slight difference can be observed as the thermal volumetric strain, ε Δ vT , is defined differently for these two elements types. Solid elements deform thermally in three dimensions (i.e. = − ε α T Δ 3 Δ vT T ), while in beam elements, temperature change can only induce thermal strains in the in-plane and out-of-plane axial directions (i.e. = − ε α T Δ 2 Δ

vT T
). Since lateral total strain is restricted in both cases, there is an additional mechanical strain which is larger for solid elements ( α T 2 Δ T ) than for beam elements (α T Δ T ). As a result, there is a larger pore water pressure increase in the former case, which translates into a slight heave in the vertical direction due to the reduction of effective stresses. However, this difference is clearly a transient phenomenon, meaning that in the long-term, once drainage is complete, there is no difference between the two types of elements.

Conclusions
This paper presents the first complete finite element formulation for a fully THM coupled 2D beam element. The full derivation process, as well as the adopted assumptions, is clearly stated before the final FE form of the formulation is presented. Conclusions can be summarised as follows: (1) The proposed THM formulation for the beam element is fully coupled and compatible with an existing coupled formulation for solid elements, thus enabling the use of both beam and solid elements in an FE mesh discretisation to simulate soil-structure interaction in coupled THM problems (e.g. thermo-active piles or retaining walls). (2) The flexibility of the proposed FE formulation allows the beam elements to be employed in the FE modelling of various types of coupled problems (e.g. HM, TM or TH) for discretising either porous or single-phase materials.