A coupled thermo-hydro-mechanical finite element formulation of one-dimensional beam elements for three-dimensional analysis

Finite element (FE) analysis in geotechnical engineering often involves entities which can be represented as onedimensional elements in three-dimensions (e.g. structural components, drains, heat exchanger pipes). Although structural components require an FE formulation accounting only for their mechanical behaviour, for the latter two examples, a coupled approach is necessary. This paper presents the first complete coupled thermo-hydromechanical FE formulation for one-dimensional beam elements for three-dimensional analysis. The possibility of deactivating each of the systems enables simulation of both coupled and uncoupled behaviour, and hence a range of engineering problems. The performance of these elements is demonstrated through various numerical simulations.


Introduction
Many geotechnical problems involve structural components such as retaining walls, tunnel linings, anchors, props, etc. which must be included in the finite element (FE) analysis. Although these components can be represented by continuum finite elements, this approach may result in a large number of elements increasing the computational effort, or in elements with unacceptable aspect ratios due to their geometries. Instead, they can be modelled using specific finite elements which are formulated by reducing one or more dimensions of the structural component to zero, and hence represent it as a line or a surface with zero thickness [1].
This study focuses on structural components that can be modelled as line (i.e. one-dimensional) elements in three-dimensional (3D) analysis. Examples of such components include anchors, props, beams, columns, piles, sand columns as well as drains. Moreover, the recent increase in the use of systems that utilise geothermal energy (i.e. ground source energy systems) has created the need for accurate prediction of the temperature field in the ground. This requires the inclusion in the numerical model of the heat exchanger pipes through which a fluid is circulated facilitating the transfer of thermal energy between the building and the ground. Due to their small diameters, they can also be represented as one-dimensional elements in 3D. Clearly, in the case of structural elements such as anchors, props, beams or columns, their mechanical behaviour is of interest to geotechnical engineers, and therefore, the FE analysis must solve the force equilibrium equations. However, in order to model drains or sand columns, a coupled hydromechanical (HM) FE formulation is required, whereas a coupled thermo-hydraulic (TH) formulation is necessary for heat exchanger pipes. Furthermore, thermal drains, which have been shown to enhance the consolidation rates in soft soils by imposing high temperatures (e.g. [2,3]), require the solution of coupled thermo-hydro-mechanical equations (THM). It should be highlighted that in all the above examples, the one-dimensional elements are used within a soil mass simulated using 3D solid elements.
The FE software PLAXIS 3D [4] uses one-dimensional elements, referred to as beam elements, for modelling structural components in 3D. These are 3-noded elements with five degrees of freedom per node, i.e. three displacements (one axial and two transverse) and two rotations. As the rotation degree of freedom around the element's axis is not considered, these beams cannot sustain torsional moments. One-dimensional beam elements with an additional rotational degree of freedom, which allows for torsion, are available in the FE programs ABAQUS [5] and CRISP [6]. Moreover, one-dimensional elements (referred to as truss elements in ABAQUS and bar elements in CRISP), which have only displacement degrees of freedom, and hence, transmit only axial forces, have also been developed. However, it should be noted that these FE formulations are purely mechanical, and therefore, the one-dimensional elements cannot be used in problems where their thermal or hydraulic behaviour is important. In order to model drains, Hird et al. [7] and Russell [8] developed a 3-noded line drainage element for use in plane strain and axisymmetric consolidation analysis in CRISP. However, as this element has a zero cross-sectional area, it is unclear whether its use can be extended to the modelling of the mechanical behaviour of structural components. Furthermore, the element cannot simulate the fully coupled HM behaviour as the off-diagonal coupling terms in this HM formulation are zero. Finally, the formulation has not been extended for 3D analysis. In PLAXIS, drains are simulated using a seepage boundary which permits water to leave the mesh at atmospheric pressure [4].
The conductive-advective flows in heat exchanger pipes (i.e. coupled TH problems) have been simulated in 3D with COMSOL Multiphysics (e.g. [9][10][11]) using a pipe interface which solves the continuity, fluid momentum and energy balance equations [12]. However, to the authors' knowledge, the coupled THM FE formulation of one-dimensional elements in 3D analysis has not been developed.
This paper presents a complete fully coupled THM FE formulation of one-dimensional beam elements for 3D analysis which has been implemented into the bespoke FE softwarethe Imperial College Finite Element Program (ICFEP, Potts and Zdravković [1]). The elements may be either 2-or 3-noded with each node possessing several degrees of freedom: three displacements, three rotations, as well as pore fluid pressure and temperature. The mechanical formulation of these new elements is based on that for Mindlin's beam elements for two-dimensional (2D) analysis developed by Day and Potts [13] which has been shown to be valid for both straight and the more problematic curved beams. As the extension of this formulation to 3D involves additional displacement and rotation degrees of freedom, the new beam elements can transmit axial and shear forces, as well as bending and torsional moments. The coupled THM formulation is based on that for continuum elements presented by Cui et al. [14], and therefore, is compatible with other finite element types in ICFEP, allowing modelling of the interactions between the one-dimensional elements and other structural components and/or the surrounding medium. Additionally, in order to overcome problems associated with modelling of highly advective flows, a Petrov-Galerkin finite element method (FEM) instead of the Galerkin FEM may be used when advection dominates heat transfer along the element. It should be noted that any of the three systems of the coupled THM formulation (i.e. mechanical, hydraulic and thermal) may be disabled such that the beam elements can be used in a variety of problems involving coupled HM, TH or TM, or uncoupled mechanical, hydraulic or thermal behaviour. To demonstrate the performance of the new one-dimensional beam elements, a series of validation exercises is presented in this paper. Finally, it should be noted that one-dimensional bar elements have also been developed, although, as they differ from beam elements only in terms of their mechanical behaviour, their formulation is presented briefly in Appendix A. The sign convention adopted throughout this paper is such that tensile stresses, strains and forces are positive.

Mechanical behaviour
Mechanically, the one-dimensional beam elements have six degrees of freedom per node (three displacements and three rotations), and transmit not only bending moments, axial and shear forces but also a torsional moment. Fig. 1 shows a 3-noded beam element of length l in 3D and defines its six degrees of freedom: displacements u, v and w in the global x, y and z directions, respectively, and rotations θ x , θ y and θ z about the global x, y and z axes, respectively. The sign convention for rotations adopted in this paper follows the right-hand rule as illustrated in Fig. 1.
For convenience, the global nodal displacements and rotations are defined as two separate vectors: Additionally, the following vectors are defined:

Constitutive equations
Under isothermal conditions, the relationship between the total strains and element forces and moments is described by the constitutive equation [1]: where for one-dimensional beam elements: σ is the vector of incremental mechanical (or total under isothermal conditions) strains whose components are: ε a σ ,axial strain, γ 1 and γ 2shear strains, χ Ttorsional strain, χ 1 and χ 2bending strains, whereas σ {Δ } is the vector of incremental forces and moments whose components are: F aaxial force, S 1 and S 2shear forces, M Ttorsional moment, M 1 and M 2bending moments. These components are related to three local axes, one of which is along the beam and the remaining two are perpendicular to the beam. For a beam with linear elastic behaviour, the constitutive matrix, D [ ], is defined as: where E is the Young's modulus, A c is the cross-sectional area, k s1 and k s2 are the shear correction factors, G is the shear modulus, J T is the torsional constant and I 1 and I 2 are the second moments of area. Note that for beams with cross-sections which are symmetric with respect to the two axis (i.e. circular or square), = = k k k s s s 1 2 and = = I I I 1 2 in Eq. (8). The torsional constant, J T , describes the torsional stiffness of the beam, and its value depends on the shape of the cross-section of the beam.
In a coupled THM problem involving a two-phase porous medium (e.g. soil), the mechanical behaviour of the material is affected by changes in pore fluid pressure, as well as temperature. In order to formulate the equations governing the mechanical behaviour, the following assumptions are introduced: (1) Changes in temperature, as well as pore fluid pressure, affect only the axial strain in a beam element. (2) There is an instantaneous temperature equilibrium between the solid and fluid phases, therefore, only the overall temperature of the medium, T , is considered here. (3) Given a change in temperature, the solid particles and the solid skeleton are assumed to have the same thermal elastic strain, ε Δ T .
Under non-isothermal conditions the incremental total strain, ε {Δ }, can be divided into two components: the incremental mechanical strain, ε {Δ } σ , and the incremental thermal strain ε {Δ } T , such that: The incremental thermal strain is defined as where α T is the linear coefficient of thermal expansion of the material and T Δ is the incremental temperature. Note that in the case of a porous material, such as soil, being modelled, α T refers to the soil skeleton.
In soil mechanics, it is often convenient to express the constitutive equation in terms of effective stresses, such that: where ′ σ {Δ } is the vector of effective forces and moments, whereas ′ D [ ] is the effective constitutive matrix. However, under non-isothermal conditions, the effect of the thermal expansion/contraction of the material must be taken into account. Furthermore, in a coupled THM formulation for a porous medium, the constitutive equation must also include the effects of the pore fluid pressure. Therefore, substituting Eq. (9) into Eq. (10), and adopting the principle of effective stress results in:

Finite element formulation
The relationship between the nodal degrees of freedom, d {Δ } n , and total strains, ε {Δ }, is described by: where B [ ] is a matrix containing derivatives of the interpolation functions. Eq. (12) can also be split into four components (where ε Δ a is the incremental total axial strain): The strains in a beam element are defined in terms of the global displacements and rotations as: By combining Eqs. (13)-(33), the B [ ] matrices associated with each integration point, i, can be defined as: where N i is the isoparametric interpolation function of the displacement degree of freedom, whereas ′ N i is its derivative with respect to the natural coordinate S. For a 2-noded element shown in Fig. 3(a), the isoparametric shape functions are: The isoparametric shape functions for a 3-noded element shown in Fig. 3 It is important to note that, when full integration (i.e. 2-point for 2noded elements or 3-point for 3-noded elements) is used, shear force locking, which manifests itself as oscillating axial and shear forces, may occur [1]. To eliminate this problem, Day [15] proposed lower order interpolation functions derived from a least squares smoothed approximation to the standard isoparametric interpolation functions. For onedimensional beams, these substitute functions, N i , are used for the interpolation of θ { } in the shear strain equations (Eq. (35)). For a 2-noded elements, the substitute interpolation functions are: Whereas for 3-noded elements: In terms of global matrices, the finite element equation associated with equilibrium can be obtained by minimising the potential energy with respect to the incremental nodal displacements and rotations, resulting in:

Fluid flow
The fluid flow along a one-dimensional beam element is described by the continuity equation for compressible pore fluids: where v f is the seepage velocity along the beam element, n is the porosity, K f is the bulk modulus of the fluid, t is time and Q f is any fluid source or sink. Since in the formulation for this type of element it is assumed that no strains in the transverse direction are calculated, the volumetric mechanical strain, ε v σ , , is reduced to the axial strain: . However, in a coupled THM formulation, the effect of temperature changes on the pore fluid pressure must be taken into account. Therefore, an additional term representing the change in volume of the fluid due to a temperature change is included in Eq. (53), such that: where α Tf is the linear thermal expansion coefficient of pore fluid. The seepage velocity is governed by Darcy's law which is expressed as: where k f is the permeability, and h i is the hydraulic head defined as: where γ f is the specific weight of the pore fluid and is the unit vector parallel, but in the opposite direction, to gravity. Substituting Eqs. (55) and (56) into the continuity equation (Eq. (54)), applying the principle of virtual work and minimising the potential energy with respect to the incremental nodal fluid pressure, results in the following finite element equation associated with fluid flow under non-isothermal conditions: where the global matrices are defined as: In order to solve Eq. (57), a time marching scheme must be adopted. Provided that the solution at time t 1 is known, the solution at time can be found assuming that: where β 1 is the time marching parameter for fluid flow, and p ({ } ) f nG 1 is the pore fluid pressures at time t 1 . It should be noted that the time matching scheme is stable provided β 1 is between 0.5 and 1.0 [16]. Substituting Eq. (64) into Eq. (57) results in:

Heat transfer
The equation governing the transfer of heat in coupled THM problems is formulated based on the following assumptions: (1) Heat may be transferred along the beam by conduction and advection. Heat transfer by radiation is assumed negligible [17] and, therefore, is not taken into account in this formulation. (2) The previously made assumption that the solid particles and the solid skeleton have the same thermal elastic strain, ε Δ T , leads to the assumption that the thermal volumetric changes of the solid fraction in a mechanically unrestrained and free draining medium do not affect the void ratio. This implies that during thermal expansion or contraction of the solid particles, the volumes of both the solid skeleton and the voids change proportionally to their initial volumes and by the same factor. Therefore, only the volumetric changes caused by changes in effective stresses can alter the void ratio.
The flow of heat along a one-dimensional beam element is described by the law of conservation of energy: where Q T is the total heat flux, Q T represents any heat source or sink, dV is the volume of the beam and Φ T is the heat content of the beam per unit volume which, when using this type of finite element to model fully saturated porous media, can be calculated as: where n is the porosity, ρ f and ρ s are the densities of pore fluid and solid particles, respectively, C pf and C ps are the specific heat capacities of pore fluid and solid particles, respectively, and T r is a reference temperature. The total heat flux, Q T , can be divided into two contributions: heat conduction, Q d , and heat advection, Q a , which are defined as: where k T is the thermal conductivity of the material. Substituting Eqs.
(67)-(69) into Eq. (66) results in: However, the porosity can also be expressed in terms of the void ratio, e, such that = + n e e /(1 ), whereas the total volume can be written in terms of the volume of the soil particles, dV s , such that = + dV ed V (1 ) s . Substituting these expressions into Eq. (70) leads to: According to assumption (2) above, dV s varies with temperature such that: where dV s0 is the initial volume of the soil particles, and ε T is the thermal strain. Substituting Eq. (72) into Eq. (71) and eliminating dV s0 results in: The assumption that the thermal strain does not affect the void ratio leads to the following relationship between the changes in void ratio and the changes in volumetric strain in a small displacement analysis: where e 0 is initial void ratio. Substituting Eq. (74) into Eq. (73) and assuming that the fluid density does not change results in: Using the principle of virtual work and minimising the potential energy with respect to the incremental nodal temperature results in the following finite element equation governing the transfer of heat: where the global matrices are defined as: In order to solve Eq. (76), a time marching scheme must be adopted. Similar to the fluid flow governing equations, it can be assumed that: where β 2 and α 1 are the time marching parameters which must be between 0.5 and 1.0 for the scheme to be stable [18]. Substituting Eqs.

Finite element formulation
The matrix form of the coupled THM formulation for one-dimensional beam elements can be obtained by assembling Eqs. (48), (65) and (85), resulting in:

Petrov-Galerkin FEM
The finite element formulation presented so far adopts the widely used Galerkin finite element method (GFEM) which is capable of simulating most problems encountered in engineering. However, it is widely known that the GFEM, where the weight functions coincide with the shape functions, may produce erroneous solutions, characterised by an oscillating temperature distribution, when modelling highly advective fluid flowsi.e. flows where advection dominates heat transfer. It has been shown that the magnitude of these unnatural oscillations increases with increasing Péclet number (e.g. [18][19][20][21]), which describes the ratio between the advective and conductive heat fluxes and can be defined as: where L is the characteristic length in the direction of flow. In the case of the finite element method, L is the element length in the direction of fluid flow. Clearly, for a problem with specific material properties (ρ f , C pf , k T ) and fluid velocity (v f ), the Péclet number can be reduced by reducing the finite element length (L). However, in problems such as those involving heat exchanger pipes, this approach may result in a large number of very small elements, and hence significantly increase the computational effort.
To overcome these issues, a new finite element method called the Petrov-Galerkin finite element method (PGFEM), which modifies the nodal weighting functions to weigh the upstream nodes more heavily than the downstream ones, has been developed. Numerous weighting functions for both one-dimensional and two-dimensional elements have been proposed in the literature (e.g. [22][23][24][25][26]). Some of these formulations resulted in oscillating and/or over-diffused solutions which reduced the reliability of the PGFEM. However, as was argued by Brooks and Hughes [27] and also found by Cui et al. [28], the PGFEM produces accurate solutions provided that appropriate weighting functions are adopted and they are applied to all terms of the advection-diffusion equation.
In order to allow modelling of highly advective flows ( > Pe 1) along the one-dimensional beam elements presented in this paper, the PGFEM proposed by Cui et al. [28] was adopted instead of the conventional GFEM to obtain the FE formulation based on the law of energy conservation (Eq. (66)). This involves using weighting functions, W { }, which are different from the temperature interpolation functions, N { } T , (while, as previously mentioned, the GFEM assumes The resulting global matrices from Eq. (85) are presented in Appendix B. The FE equations governing the mechanical equilibrium and the fluid flow remain in the conventional Galerkin form. The PG formulation for one-dimensional elements presented by Cui et al. [28] adopts weighting functions for linear (i.e. 2-noded) elements developed by Huyakorn [23], whereas those for quadratic (i.e. 3-noded) elements were originally proposed by Heinrich and Zienkiewicz [29]. For completeness, the weighting functions are presented in Appendix B.

Validation
This section presents a number of validation exercises carried out to demonstrate the performance of the newly developed one-dimensional beam elements. The results obtained using ICFEP were compared against analytical solutions or other numerical solutions found in the literature. It should be noted that the PGFEM was adopted in the analyses involving highly advective flows (Section 3.4), whereas in all other analyses (Sections 3.1-3.3), the conventional GFEM was used.

Mechanical behaviour
In order to validate the mechanical behaviour, two validation exercises, which were originally performed by Day and Potts [13] with beam elements for 2D analysis, were carried out using the one-dimensional beam elements developed here. These particular tests were chosen as Day and Potts [13] demonstrated that previous formulations for beams for 2D analysis found in the literature were not able to simulate the behaviour of curved beams correctly. In both tests, 3-noded elements with reduced integration were used in order to prevent locking problems.
The first test involved rigid body displacement of the curved beam shown in Fig. 4. Points 1 and 3 were restrained in the y-direction and, given that this was a 3D analysis, an additional restraint preventing displacements along the z-direction (i.e. perpendicular to the plane shown in Fig. 4) Computers and Geotechnics 104 (2018) [29][30][31][32][33][34][35][36][37][38][39][40][41] x-direction was applied at point 1. The analysis was repeated with 1 and 10 elements. In both cases, the horizontal displacement was 1 m at all nodes, whereas all other degrees of freedom were zero. Moreover, all strain components were equal to zero, as expected from a rigid body movement. It should be noted that this test was repeated with different orientations of the beam in 3D space, yielding the exact solution every time.
In the second test, the cantilever shown in Fig. 5 was subjected to point loads of −10 kN in the xand y-directions at point 3. The beam's material properties are listed in Table 1. At point 1, all displacements and rotations were prevented. Analyses with either 2 or 10 elements were performed. Fig. 6 plots the nodal rotations (around the z -axis in the 3D case) along the cantilever beam, whereas Table 2 lists the reactions at point 1 and displacements at points 2 and 3. For comparison, the results of the equivalent 2D plane strain analyses performed by Day and Potts [13] are also shown. It is clear that the results of all numerical analyses are in excellent agreement with the exact solution from Roark and Young [30], even when only two elements are used. It should also be noted that the results of 2D and 3D analyses with the same number of elements are exactly the same.

Isothermal consolidation
The exercise used for validation of isothermal coupled consolidation along one-dimensional beams is based on the tests firstly carried out by Aboustit et al. [31] and Aboustit et al. [32], and subsequently used to validate other finite element software (e.g. [14,[33][34][35][36]). Although the original test was performed with 2D solid finite elements, the same problem is represented here in 3D with the new one-dimensional beam elements. The adopted mesh is shown in Fig. 7. The material properties used (listed in Table 3) are the same as in the original test, with the exception of the Poisson's ratio, which is zero for the one-dimensional beam elements (i.e. there is no Poisson's ratio effect in the one-dimensional beam formulation). The displacement boundary conditions      were prescribed such that only 1D deformation was allowed. A downward vertical load of 2 N (equivalent to the stress of 1 Pa used in the original analysis) was applied at point 1. Additionally, the pore fluid pressure at point 1 was not allowed to change throughout the analysis and a no fluid flow boundary condition was prescribed at point 2. As in the original test, gravity was set in the out-of-plane direction. 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 Cui et al. [37] and Cui et al. [14].
The vertical displacement of point 1 obtained using ICFEP was compared to the analytical solution presented in Biot [38]. Fig. 8 shows an excellent agreement between the numerical and analytical results demonstrating the validity of the coupled hydro-mechanical formulation.

Non-isothermal consolidation
The coupled THM behaviour of one-dimensional beams was validated through an exercise involving non-isothermal consolidation. As there is no exact solution to this problem, the ICFEP results can only be compared with the solutions of other FE programs. Cui et al. [14] verified the coupled THM equations for 2D solid elements in ICFEP by carrying out an analysis of non-isothermal consolidation which was first performed by Aboustit et al. [31] and Aboustit et al. [32], and subsequently by Noorishad et al. [33], Lewis et al. [34], as well as Gatmiri and Delage [35]. The ICFEP solution was then compared to the results presented in these publications showing a good agreement with Aboustit et al. [31], Noorishad et al. [33] and Lewis et al. [34]. It should be noted that the results of Gatmiri and Delage [35] did not agree with the other FE software in either the coupled isothermal consolidation exercise described previously nor this coupled non-isothermal consolidation analysis.
The coupled THM behaviour of one-dimensional beam elements was therefore validated against the solutions obtained using 2D and 3D solid elements in ICFEP. The mesh used in the analysis is shown in Fig. 7. The same material properties, listed in Table 4, as those in the original analysis by Aboustit et al. [31] were adopted. However, it should be noted that the Poisson's ratio of the solid was set to zero in order for their behaviour to be more comparable to that of the one-dimensional beam elements presented in this paper (since the present formulation assumes that no strains in the transverse direction are calculated). Similar to the previously described isothermal consolidation test, the mesh in this exercise was restrained to allow only 1D deformation. At point 1, a downward point load of 2 N and no change in pore fluid pressure were prescribed. In order to include the thermal effects, a constant temperature of 50°C was applied at point 1. No fluid or heat flow was allowed at point 2 and gravity was set in the out-of-plane direction. The time-step size was 1 s for the first 10 increments and 10 s for the remainder of the analysis as recommended by Cui et al. [14]. Fig. 9 presents the vertical displacement at point 1 observed in ICFEP analyses with different element types. In order to illustrate the effect of temperature, the results of the isothermal consolidation exercise (Fig. 8) are also plotted. As time passes, the volume of material subjected to a rise in temperature increases. The volumetric expansion associated with this temperature increase leads to an upward movement which partially offsets the settlement due to the consolidation process. This is especially visible at larger times, when the solution of  Table 3 Material properties for validation of isothermal coupled consolidation behaviour [31].  ). Due to the fact that only 1D deformation was allowed in all analyses, the vertical thermal strain is the same (i.e. ) and the horizontal total strain is zero for all element types. This implies that different elements experience different horizontal mechanical strains which are equal to − α T 2 Δ T and 0 for solid and 3D beam elements, respectively. Therefore, as the temperature increases and the material tries to expand, the horizontal restriction results in a reduction in the horizontal total stress and pore fluid pressure (i.e. these stresses become more compressive), in solid elements. In the vertical direction, the reduction in pore fluid pressure leads to an increase in the vertical effective stress (i.e. vertical effective stress becomes less compressive), and hence, an additional upward displacement. Conversely, the one-dimensional beam elements do not experience changes in pore fluid pressure due to temperature changes because of the assumption of the same thermal expansion coefficient of solid particles and pore fluid adopted in this problem. These mechanisms are further illustrated in Fig. 10 which plots the change in pore fluid pressure ( p Δ f ), the change in vertical effective stress ( ′ σ Δ y ) and the vertical mechanical strain ( ε Δ σ y , ) at point 3 (see Fig. 7) over time. Initially, the pore fluid pressure reduces by 1 kPa in all cases due to the application of the point load. As time passes, these excess pore fluid pressures dissipate such that the vertical effective stress reduces (i.e. becomes more compressive). The rate of pore fluid pressure dissipation (and therefore rate of reduction in vertical effective stress) is lower in the solid elements due to the additional excess pore fluid pressures resulting from lateral confinement. These additional excess pore fluid pressures lead to swelling of the material, and therefore, a lower rate of reduction of the vertical mechanical strain (which is compressive due to the application of the point load). Hence, by taking the displacement during isothermal consolidation as a baseline, the largest upward displacement is observed in the solid elements and the smallest in the onedimensional beam elements. Finally, it is interesting to note that, as all excess pore fluid pressures dissipate, the same final vertical displacement is measured in all analyses, with the thermally-induced component being given by α TL Δ T (where L is the height of the column).

Highly advective flow
The final validation exercises were performed in order to demonstrate the ability of the new one-dimensional beam elements to simulate highly advective flows, and therefore, its applicability to the modelling of heat exchanger pipes. In the following analyses, both the GFEM and the PGFEM were adopted. Fig. 11 shows the beam used in the 3D analyses. The conductiveadvective flow of heat from left to right was simulated by applying pore fluid pressures, p f , at both ends of the beam to induce a fluid flow with a constant and uniform velocity, as well as a heat source in the form of a fixed temperature of 10°C at point 1. At point 2, a fixed temperature of 0°C (i.e. no change from the initial temperature) was prescribed. The   material properties used in all analyses presented in this section are listed in Table 5. Analyses with Péclet numbers of 10, 100 and 10,000 were performed which was achieved by varying the fluid velocity and number of elements, n e , in the mesh. Full integration (i.e. 2-point for 2noded and 3-point for 3-noded elements) was chosen in all analyses. Gravity was set in a direction perpendicular to the beam so that it does not interfere with the fluid flow. Cui et al. [18] observed that the temperature oscillations produced by the GFEM occur only once the heat front reaches the boundary of the mesh (i.e. point 2). This is also demonstrated in Fig. 12 which shows that the results of the transient stage are free of oscillations. It should be noted that the numerical results of the transient stage (obtained using both the GFEM and the PGFEM) were compared against the approximate solution of the advection-diffusion equation by van Genuchten and Alves [39] showing an exact match, and hence verifying the accuracy of both the GFEM and the PGFEM for the transient stage of the analysis. When a fixed temperature is prescribed at point 2, the oscillations increase once the heat front reaches the mesh boundary until a steady state is achieved with maximum amplitude of oscillations. It is also clear that adopting the PGFEM eliminates this problem and results in an exact solution of a uniform temperature along the beam except for the last node where a no change from the initial condition is prescribed. The results of this validation exercise clearly demonstrate the effectiveness of the adopted PGFEM in the context of advection-dominated heat flux problems, such as those where the proposed element is used to simulate heat exchanger pipes.

Conclusions
This paper presents a coupled THM FE formulation of a one-dimensional beam element for 3D analysis. The key conclusions can be summarised as follows: (1) The mechanical behaviour of the new beam element was formulated based on that of Mindlin's beam elements for 2D analysis proposed by Day and Potts [13]. It was shown that the new onedimensional beam performs well as either a straight or curved beam. (2) The coupled THM formulation is based on that for continuum elements developed and implemented into ICFEP by Cui et al. [14]. Therefore, compatibility of the one-dimensional beam element with other finite element types available in this code is ensured. This allows modelling the interactions between one-dimensional beam elements and other structural elements, as well as the surrounding medium (e.g. the soil mass). (3) In order to enable modelling of highly advective flows such as those present in heat exchanger pipes, the Petrov-Galerkin FEM proposed by Cui et al. [28] was adopted for the solution of the energy conservation equation instead of the conventional Galerkin FEM. This method adopts weighting functions which are different from the temperature interpolation functions and produces solutions which are accurate and free of numerical oscillations. (4) The one-dimensional beams may behave as a two-phase porous material or a single-phase material depending on the properties chosen for the two phases. Moreover, as any of the three systems in the THM formulation may be disabled, the beam's behaviour can be modelled as coupled HM, TH or TM, or uncoupled mechanical, hydraulic or thermal. (5) The validation exercises presented showed an excellent match between the computed results and the exact or approximate solutions. In cases where such solutions were not available, the numerical results compared well with other numerical solutions obtained using other element types.
where β PG1 and β PG2 are the PG weighting factors and the function g S ( ) was defined as: The optimal value of the PG weighting factors β PG1 and β PG2 can be obtained from: The above weighting functions were first proposed by Heinrich and Zienkiewicz [29] who successfully simulated 1D steady state advection-diffusion problems with 3-noded quadratic elements. It should be noted that behaviour of these weighting function in solving 1D transient advection-diffusion problems has not been shown elsewhere in the literature.   B1 presents the results of the validation exercise described in Section 3.4, with Péclet numbers of 10 and 10,000 where a fixed temperature boundary condition was prescribed at point 2 (see Fig. 11). Although the oscillations at the steady state in the results obtained with the GFEM increase with Péclet number, the PGFEM produces the correct solution.