Hybrid finite elements for nonlinear thermal and hygral problems

Computers and Structures 182 (2017) 14–25

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Hybrid finite elements for nonlinear thermal and hygral problems J.A. Teixeira de Freitas a,⇑, P.T. Cuong a,b, Rui Faria c a

Departamento de Engenharia Civil, Arquitectura e Georecursos, CERIS, Instituto Superior Técnico, Universidade de Lisboa, Portugal Faculty of Civil Engineering, Ho Chi Minh City, University of Transport, Viet Name c Departamento de Engenharia Civil, Faculty of Engineering, University of Porto, Portugal b

a r t i c l e

i n f o

Article history: Received 27 July 2016 Accepted 17 November 2016

Keywords: Hybrid finite elements Heat transfer Thermal singularities Moisture transport

a b s t r a c t The performance of a hybrid finite element formulation developed for the solution of nonlinear transient problems is assessed. The formulation, already applied to the solution of heat transfer and moisture transport problems, develops from the uncoupling of the approximation of the state variables and the mapping of the geometry. The formulation qualifies as hybrid because the state variable and its gradient are approximated independently. Orthogonal bases are used to enhance numerical stability under highorder approximations. The resulting solving system is highly sparse and well-suited to adaptive refinement and parallel processing. Besides the benchmarks used in the assessment of conform elements, the rates and patterns of convergence of the hybrid element are defined and its sensitivity to shape distortion is analysed. Also illustrated is the simulation of singular heat flow fields in cracked plates and the hygro-thermo-chemical modelling of cement hydration in concrete structures. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Nonlinear transient problems are usually solved using the conventional (conform) formulation of the finite element method (FEM) ever since the applications of Zienkiewicz and Cheung [1], as well as of Bathe and Khoshgoftaar [2] to heat transfer. Although literature is rich in using alternative hybrid formulations of FEM in a wide range of application fields, their use in the solution of heat transfer problems has been very limited since the work of Fraeijs de Veubeke and Hogge [3] on steady-state applications. Besides its extension to transient analyses using the Trefftz variant [4], reports on the application of hybrid formulations to thermal problems remain limited in scope, restricted to thermo-elastic fracture, inverse heat conduction and convection-diffusion problems [5–7]. The modelling flexibility offered by hybrid formulations is substantially increased when the approximations of geometry and state variables are uncoupled. They can be implemented using unstructured, coarse meshes of high-degree elements, to yield sparse solving systems well-suited to adaptive refinement and parallel processing. Examples are the use of hybrid elements for simulating cement hydration in early age concrete structures and in fire protection of structural elements strengthened with CFRP laminates [8–10]. ⇑ Corresponding author. E-mail address: [email protected] (J.A. Teixeira de Freitas). http://dx.doi.org/10.1016/j.compstruc.2016.11.009 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

Those reports focus on assessing the use of conform and hybrid variants of FEM in the solution of nonlinear and coupled transient problems. The focus of the present paper is distinct: it aims to characterize the performance of the hybrid element in terms of rates and patterns of convergence, under both p- and hrefinement, and its sensitivity to gross mesh distortion. A benchmark test set in [2] is also assessed to illustrate an axisymmetric application. These tests are complemented with the modelling of temperature fields in cracked plates, to illustrate the extension of the finite element approximation basis to model singular heat flows. Numerical testing closes with the simulation of cement hydration in a concrete block to show that the formulation written here for heat transfer problems can be readily extended, in this particular case to hygro-thermo-chemical analyses. Besides defining the mathematical model, the paper addresses the procedures used to approximate the geometry and the state variables, to characterize the approximation bases and to discuss the properties of the resulting solving system. Comments on numerical implementation are limited to essential aspects, as the detailed information presented in [8,10] remains valid for the present extension to singular flux fields. For simplicity, the mathematical model, the finite element formulation and all testing problems are presented for heat transfer problems, with the exception of the closing application, which is used to comment the extension to moisture transport problems.

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

15

2. Mathematical model

3. Discretization in space and time

The domain conditions that govern the thermal field, Tðx; tÞ, in the domain V e of a finite element with boundary Ce , assigned to space and time systems of reference x and t, respectively, can be summarized as follows:

The hybrid formulation presented below is implemented on coarse meshes of high-degree elements. Most elements are regular in shape, as the domain decomposition is mainly constrained by the medium geometry and by inhomogeneity of its properties. Although no conceptual constraints limit the topology of hybrid elements, e.g. [11], domain decomposition has been implemented using parametric descriptions,

$T r þ cT_ ¼ Q_

in V e

ð1Þ

e ¼ $T in V e r ¼ ke in V

ð2Þ e

ð3Þ

x ¼ NðnÞc

ð16Þ

in V

In the thermal equilibrium Eq. (1), $T r is the divergence of the heat flow, c is the volumetric specific heat and Q_ represents the heat source. Eq. (2) defines the temperature gradient, e, and Eq. (3) represents the constitutive relation, as k is the local conductivity matrix. Conditions (1)–(3) are often combined to obtain the state equation:

where vector c defines the coordinates of the element control nodes, N is the shape function matrix and n represents the parametric co-ordinate system. It is noted that mapping (16) is independent of the approximation of the state variables of the problem, as shown in Section 4. Following the usual practice in the solution of the first-order problems, time integration of temperature T is based on trapezoidal rules with a fixed time increment, dt,

$T k$T þ Q_ ¼ cT_ in V e

T ¼ T 0 þ c0 dtT_ 0 þ cdt T_

ð4Þ

ð17Þ

The mathematical model is completed with the following set of boundary and initial conditions:

where c and c0 are integration constants (c ¼ c0 ¼ 12 in the tests presented below).

nT r ¼ nT k$T ¼ qN

4. Finite element approximations

T ¼ TD

on C

on CeN

ð5Þ

e D

T ¼ T 0 at t ¼ t 0

ð6Þ in V e

ð7Þ

Eqs. (5) and (6) assume that the element boundary, with a unit outward normal n, is uncoupled into complementary Neumann and Dirichlet parts,

Ce ¼ CeN [ CeD

ð8Þ

which are decomposed as follows, under the notation defined below:

CeN ¼ Ceq [ Cecr [ Cecd

ð9Þ

CeD ¼ CeT [ Cei

ð10Þ

The Neumann condition (5) is used to model prescribed heat flux fields, q, convection-radiation and conductance conditions:

qN ¼ q on Ceq

ð11Þ

qN ¼ hcr ðT  T a Þ on Cecr

ð12Þ

e cd

ð13Þ

qN ¼ hcd ðT  T k Þ on C

In the equations above, Ta is the ambient temperature and Tk the surface temperature of a connecting element. The radiation condi2

T 2a Þ

tion is set letting hcr ¼ hr ðT þ T k ÞðT þ in Eq. (12), with hr representing the radiation coefficient. The Dirichlet Eq. (6) is used to model prescribed temperature fields, T, and inter-element thermal continuity:

TD ¼ T TD ¼ Tk

on CeT on Cei

ð14Þ ð15Þ

The prescribed temperature and heat flux fields and the ambient temperature are defined as functions of time and space. The specific heat and the conductivity coefficients, as well as the convection, radiation and conductance coefficients, are defined in a similar way and may vary with temperature.

Either the temperature or the heat flow fields can be selected for primary approximation in a hybrid formulation of thermal problems. The first option, also used in the derivation of conform elements, is designed to locally enforce the thermal gradient and the conductivity conditions (2) and (3), while the second is designed to satisfy Eq. (3) and the thermal equilibrium condition (1). These are the approaches followed in the early contribution of Fraeijs de Veubeke and Hogge [3], with the aim of bounding the error of FEM solutions of steady-state heat conduction problems. The first option is chosen here because in most thermal applications the main concern of the analyst is to directly control the quality of the temperature approximation. It is written assuming separation of variables in time and space,

Tðx; tÞ ¼ HðnÞTðtÞ in V e

ð18Þ

where vector T defines the amplitudes of the temperature modes listed in row-vector H. The temperature approximation basis is defined as the product of (naturally hierarchical, orthogonal) Legendre polynomials,

H ¼ f    Li ðnÞLj ðgÞLk ðfÞ    g

ð19Þ

expressed in the natural element co-ordinate system. The dimension of this basis is,

ND ¼

D Y j¼1

(

1 ðd D!

þ jÞ

if dn ¼ dg ¼ df ¼ d

ðdj þ 1Þ

if independently set

ð20Þ

for two-dimensional (D ¼ 2) and three-dimensional (D ¼ 3), depending on whether the same or different degrees of approximation are implemented in each direction. The same basis is used to approximate the temperature rate and, as for conform elements, the thermal gradient and the heat flow rate fields are determined enforcing conditions (2) and (3) for the assumed temperature approximation:

T_ ¼ HT_

in V e

e ¼ ð$HÞT in V e

ð21Þ ð22Þ

16

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

r ¼ kð$HÞT in V e

ð23Þ

However, and in consequence of generalizing the temperature approximation basis H, which is fully independent of the geometry mapping basis, N in Eq. (16), the implicit enforcement of thermal continuity is no longer viable and implies the independent approximation of the heat flux on the Dirichlet boundary of the element, see definition (10),

qðx; tÞ ¼ uðnÞqðtÞ on CeD

ð24Þ

where n is the boundary co-ordinate system and vector q defines the amplitudes of the heat flux modes listed in row-vector u:



u ¼    Li ðnÞLj ðgÞ   



ð25Þ

The Legendre polynomials are defined in the surface coordinates of the element and the dimension of the basis on a given surface of the Dirichlet boundary of the mesh is still defined by Eq. (20) replacing D by D  1. The supporting code has the option of implementing the approximation bases (19) and (25) using Chebyshev polynomials, which are numerically more stable when (very) high degrees of approximation are used. Independently of the type of the approximation functions, the first option in definition (20) is used in structured meshes with elements regular in shape. The second option is useful in the implementation of unstructured meshes and elements of high aspect ratios, namely in applications involving elements with distinct geometries and thermal properties [8]. Bases (19) and (25) can be implemented with any type of element, as their definitions are independent of the parametrization technique used to define the geometry of the element. 5. Finite element solving system The alternative procedures to derive the weak form of the state Eq. (4) for conform elements,

H T_ þ ðK þ CÞT þ Q D ¼ Q_ þ Q k  Q

Z

HT cHdV Z

e

ð27Þ

ð$HÞT kð$HÞdV

K¼ Z

e

ð28Þ

Z HT hcr HdCecr þ

C¼

HT hcd HdCecd

ð29Þ

as well as for the heat source vector and for the prescribed and unknown heat flux terms on the Neumann and Dirichlet boundaries (9) and (10), the former under conditions (11)–(13):

Q_ ¼

Z

Z QD ¼

HT Q_ dV

e

ð30Þ

HT nT rdCeD

ð31Þ

HT hcd T k dCecd

ð32Þ

Z Qk ¼ Z Q¼

Z T

e q

H qdC 

T

H hcd T a dC

e cd

Q D ¼ Bq

ð33Þ

ð34Þ

where the heat flux matrix is defined as follows, under identification (10):

Z

HT udCeD

B¼

ð35Þ

The approximation basis is used to enforce on average the Dirichlet condition (6),

Z

uT ðT  T D ÞdCeD ¼ 0

ð36Þ

and the weak form of thermal continuity is obtained enforcing in the equation above the temperature approximation (18) under boundary conditions (14), (15) and definition (35),

BT T ¼ T þ Rk

ð26Þ

can be directly applied to establish the equivalent equation for hybrid elements, under approximations (18) and (21) for the temperature and temperature rate fields. Similar expressions hold for the specific heat, conductivity and convection/conductance matrices,

H¼

Consequent upon the choice for the approximation bases, which no longer involve nodal (Lagrange-type) functions, neither the quantities defined above represent equivalent nodal resultants, nor is it possible to enforce thermal continuity in strong form simply by equating the nodal temperatures of connecting elements: instead of nodal values, the coefficients of vector T represent now the weights of the functions used to set-up the approximation basis. Therefore, the elementary system expressed in Eq. (26) can no longer be assembled through nodal incidence of temperatures and equilibrium of the nodal resultants of the heat flux on the Dirichlet boundary (10) of the element. This is what justifies the need for the independent approximation of the heat flux defined by Eq. (24), and thus the hybrid nature of the finite element formulation used here. The heat flux resultants on the Dirichlet boundary (31) become explicit variables of the weak form of thermal equilibrium (26),

Z

ð37Þ

uT TdCeT

T¼ Z Rk ¼

ð38Þ

uT T k dCei

ð39Þ

where Tk, the surface temperature in connecting elements, is also approximated in form (18). The nodal summation operations that typify the assemblage of the conform elements are no longer present in the combination of the elementary Eqs. (26) and (37). Assemblage consists, simply, in organizing those equations according to the listing in vectors T and q of the weights of the temperature and heat flux modes in all elements and Dirichlet boundaries of the mesh, to yield the symmetric, nonlinear and first-order solving system:



H O

( _ )  K þC T þ O BT q_

O

B

  T

O

q

( ¼

Q_  Q T

) ð40Þ

Consequently, all matrices and vectors in system (40) are strictly element dependent (matrices H, K and C are blockdiagonal), with the exception of the flux matrix B, as its coefficients are shared by at most two connecting elements. 6. Numerical implementation The first-order nonlinear system (40) is solved combining the time integration rule (17) with the Newton-Raphson method, which is modified to preserve symmetry. The following gradient form is obtained for the n-th iteration,

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25



Sn BT

B O

(

T~ n ~n q

)

( ¼

R0en R0cn

) ð41Þ

where the tilde identifies the residual variable (for instance, ~ n ). The non-null coefficients of system (41) are calcuT n ¼ T n1 þ T lated once for linear terms (as for the heat flux matrix, B, and often for the convection matrix), every time-step for time-dependent terms (typically the ones associated with prescribed fluxes and temperatures), and at every iteration for temperature-dependent terms (for instance, those associated with radiation conditions). Particular attention is paid to scaling. The problem is first nondimensionalized using the characteristic values for length, temperature and thermal properties. The approximation bases are scaled at element level to normalize the conductivity and heat flow matrices (to stabilize numerical integration). Finally, the solving system is independently normalized (to stabilize the implementation of the solver). Under this normalization, the numerical zero is set to 106 in the control of the convergence of the Newton-Raphson cycles and time integration steps. There are no constraints on the design of the mesh. Besides the combination of different types of master-elements in the same mesh to obtain adequate representations of curved surfaces, the uncoupling of geometry and state variable approximations simplifies the extension of mapping (16) to generate unstructured meshes. This technique, illustrated in [8], is useful to reduce both the meshing effort and the number of degrees-of-freedom in the analysis of structures combining elements with distinct material and geometric profiles. It is also useful to maximize the number of regular rectangular/hexahedral elements necessary to exploit the full (odd-even) orthogonality of Legendre (Chebyshev) bases. Finite element meshes that evolve in time are implemented in the supporting computer code to simulate the staged construction of structures. At the onset of a new construction stage, the weighting vectors present in the temperature, temperature rate and heat flux approximations (18), (21) and (24) are known, and the corresponding information has to be determined for the elements used to simulate the new stage of construction. The implementation of this technique is illustrated in [8]. The user is free to select the degrees of the domain and boundary bases (18) and (24), in each element and on each boundary, and in each direction of the space, as implied by definition (20). There are, however, linear independence constraints on the pairing of the polynomial degrees for the temperature and heat flux fields. The thermal continuity condition, defined by the second equation in system (40), shows that the degree in the boundary (heat flux) approximation (24) cannot exceed the degree of the domain (temperature) approximation (18) mapped on that boundary. Under this condition and the basic requirements of completeness and linear independence of the individual element bases, spurious solutions (linear dependency at mesh level) may be avoided by ensuring that the element is thermally undetermined. According to the same continuity condition written for the extreme case of a single-element under Dirichlet boundary conditions, this implies that the dimension of its domain basis should exceed the sum of the dimensions of the approximations on its boundary, as defined by Eq. (20). The implementation of these criteria on the selection of the domain and boundary degrees of approximation has proved sufficient to avoid the mobilization of spurious solution modes [8]. A variable time step can be used throughout the analysis, as its duration may vary from seconds to months in different applications [8–10] with substantially different levels of nonlinearity in the response throughout that period. The process is considered to

17

be stable when a time increment is implemented with up to three Newton-Raphson iterations. To control the adequacy of the time-step, the integration rule (17) is used in the iterative process to update the temperature rate _ and the value obtained is compared, upon convergence, vector, T, with the solution of the first (thermal balance) equation of system (40), which is solved at element level for the current T and q fields. Because of the non-nodal nature of the approximation bases, the initial value of the temperature weight vector, T0, is determined enforcing on average condition (7):

S0 T 0 ¼ t0 Z HT HdV

S0 ¼

ð42Þ e

ð43Þ

e

ð44Þ

Z t0 ¼

HT T 0 dV

System (40) is solved at element level to determine the weights for the initial heat flux and temperature rate fields, q0 and T_ 0 . In what regards the structure of the elementary finite element matrices, definitions (35) and (43) show that matrices B and S0 are linear. The specific heat and convection/conductance matrices K and C, defined by Eqs. (27) and (29), are also linear in some applications, for instance in the modelling of cement hydration. When they are implemented on rectangular/hexahedral elements, matrices K, C and S0 are diagonal and matrix B is also highly sparse if Legendre bases are used. Moreover, their coefficients have simple analytic expressions. Thus the emphasis placed on the maximization of regular elements in the mesh. The fill-in caused by different sources of nonlinearity is weaker than that caused by shape distortion under linear conditions, which has a stronger effect on their condition numbers. The key property of the structure of the finite element solving system (40) is that the temperature degrees-of-freedom, T, are strictly element-dependent and the heat flux degrees-of-freedom, q, are only shared by connecting elements. This is why system (40) is well-suited to parallelization (the management of shared information is minimized) and adaptive p-refinement (it suffices to compute the additional coefficients of the matrices and vectors associated with the added degrees-of-freedom). A third consequence is that the system is highly sparse, even when irregular elements are used. Orthogonality of the temperature and heat flux approximation bases and, in particular, the option of not condensing the solving system in the temperature degrees-of-freedom leads to relatively high sparsity indices. They are in the 80–90% range in the tests presented below, because the dimension of the systems is relatively small. They increase substantially with that dimension and may easily reach the level of 99%, as shown in [8]. Naturally, the features summarized above features are relevant only in the solution of large modelling problems. They are seriously compromised when system (40) is pre-condensed, at element level, on the boundary degrees-of-freedom, the heat flux amplitudes q. 7. Numerical tests The solutions obtained for six testing problems are presented next, namely four performance tests and two applications, one of the latter on the modelling of singular heat flow fields and the other on the implementation of coupled problems. The first three tests assess the performance of the hybrid element in terms of rates and patterns of convergence, as well as sensitivity to mesh distortion. The fourth test is chosen to solve with hybrid elements

18

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

a benchmark for conform elements set by Bathe and Khoshgoftaar [2]. The first application problem is selected to illustrate two features of the approach followed here, namely to enrich the approximation bases to model singular heat flow fields in a cracked plate implementing approximations bases with overlapping supports. The closing test, on the simulation of the hygro-thermo-chemical process of cement hydration in a concrete block cast on a rock foundation, is used to show that the formulation can be readily applied and coupled to other than heat transfer problems. No smoothing is used in all graphs presented below.

7.1. Convergence

ordinate system. The effect it may have on the stability of the basis is assessed next.

7.2. Sensitivity to shape distortion Distortion sensitivity of hybrid-Trefftz elements for transient heat conduction problems has been reported in [4] using a regular and smooth problem. The option here has been to use the more demanding problem (45). The distortion scheme defined in Fig. 2a) for a 2  2 element mesh is applied, and sensitivity is measured by the energy variation with respect to the value obtained in the undistorted mesh:

Z

A one-dimensional test is used in [7] to illustrate the convergence of a hybrid-mixed finite element formulation for convection-diffusion applications under local Peclet numbers associated with boundary layer effects. This test is extended here into two dimensions,

E¼

coshðx=kÞ sinhðy=kÞ  Tðx; yÞ ¼ coshð1=kÞ sinhð1=kÞ

k ¼ 101 , using two approximation bases, with degrees d ¼ 7 and d ¼ 11, independently of the level of distortion. The error in energy,

ð$TÞT kð$TÞdV

ð46Þ

ed ¼ Abs½1  EFE ðdÞ=EFE ðd ¼ 0Þ

ð47Þ

The results presented in Fig. 2(b) are obtained for the test case

ð45Þ

and defined on a unit square domain ½x; y ¼ ½0; 1. This solution is used to prescribe Neumann conditions on all sides, except on side y ¼ 0, where a Dirichlet condition is assumed and the independent heat flux approximation (24) is implemented. Problem (45) is also used below to assess rates of convergence and sensitivity to mesh distortion. The problem is first solved using a single-element mesh with decreasing values of coefficient k in definition (45) to cause the formation of an exponential temperature variation in the vicinity of vertex x ¼ y ¼ 1, where the peak temperature occurs, T max ¼ 1. The results obtained are presented in Fig. 1a). The polynomial basis easily captures the shape of the exact solution, with small oscillations in the flat solution zone, which decrease with the level of p-refinement. This can be verified, for instance, solving test case k ¼ 102 with an approximation of degree d ¼ 11. According to Eq. (20), the degrees-of-freedom are N 2 ¼ 78 for temperature and N 1 ¼ 12 for heat flux. The difficulty is to model the gradient of the function and, consequently, to converge to its peak value, as shown in Fig. 1(b) for three test cases and approximation degrees in the range 1 6 d 6 31. Numerical stability under high-degree approximations is a direct consequence of the orthogonal basis being used. However, orthogonality of the approximation basis is destroyed by distorted geometric mappings (16), as the basis is defined in the natural co-

1.0

T

Section y = x

and eE  2  1010 for approximation d ¼ 11 (close to the converged solution, for double precision calculation). These relative levels of convergence can be assessed in terms of the peak temperature using the results presented in Fig. 1(b). The graphs in Fig. 2(b) recover results that have been wellestablished in the assessment of hybrid elements for solid mechanics applications, e.g. [11]: the highest values of the sensitivity parameter, ed , remain marginal for strongly distorted elements, with d ¼ 0:49 (when the areas of two elements of the mesh are close to vanish and the shape of the other two degenerate from squares to triangles); the level of sensitivity decreases substantially with the quality of the convergence of the solution (the maximum sensitivity has a two-order decrease when the degree of the approximation is improved from d = 7 to d = 11). In this application, sensitivity is higher when d ! 0:49 (lower when d ! þ0:49) because the region where the highest gradient of the function develops is modelled by an element that is increasing (decreasing) in area, under an approximation basis that remains unchanged during the test.

Tmax

Point x = y = 1 k = 10 2

k = 10 1

1

0.8

0.6

ð48Þ

is eE  1:5  104 for approximation d ¼ 7 (solution yet to converge)

1.0

Exact Hybrid

0.8

eE ¼ 1  EFE =Eexact

0.6

(d = d = d = 31)

0.4

0.4

k= 1

0.2

10

2

x

0.0 0.2

0.4

0.6

0.8

(d = d = d)

0.2

10 1

0.0

Exact solution: Tmax = 1

1.0

(a) Temperature profile on diagonal section

d

0.0 1

6

11

16

21

26

(b) Convergence to peak temperature

Fig. 1. Single-element mesh solution.

31

19

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

y

δ ≥0

k = 10 1

3.0 2.4

0.5

δ ≥0 0.5

1.8

102⋅

(d = d = d = 7)

1.2

104 ⋅

(d = d = d = 11)

0.6

0.5

x

0.5

0.0 -0.5

-0.3

(a) Distortion scheme

-0.1

0.1

0.3

0.5

(b) Variation in energy

Fig. 2. Sensitivity to gross mesh distortion.

7.3. Rates of convergence

Table 1 Rate of convergence,a r.

The pattern of convergence and the rate of convergence, r, as measured by relation eE ¼ r 0 N r between the error in energy (48) and the number of degrees-of-freedom of system (40), N, are shown in Fig. 3 and defined in Table 1. The degrees of approximation, with dn ¼ dg ¼ dn ¼ d as in the previous example, are tested in the range 3 6 d 6 11; h-refinement is implemented sub-dividing the single-element mesh into regular meshes with 2  2 and 4  4 elements. The range in the number of degrees-of-freedom is 14 6 N 6 1584. The results obtained for test case k ¼ 101 confirm that convergence under p-refinement is stronger than under h-refinement. This test can also be used to show that convergence is weaker (stronger) for ill- (well-) behaved problems, that is, for k ¼ 102 (k ¼ 1). This latter exponential layer problem can only be adequately solved using high levels of either h- or p-refinement in the boundary layer region, namely by enriching the basis with functions that model high-gradient fields. They can be implemented in the framework of hybrid finite element formulations, as illustrated below for singular heat flow fields. 7.4. Axisymmetric problem The tests used in Bathe and Khoshgoftaar [2] cover a broad range of problems and became benchmarks for conform elements. The axisymmetric test consists of a solid cylinder, with radius L, height 2h and initial temperature T 0 . The cylinder is allowed to cool in air under the temperature relation T a ¼ 0:55T 0 , in Kelvin.

ln( E) 4

3 d =d =d = 3 d =7

-6 Meshes: 1×1 2×2 ♦ 4×4

-16

7

d = 11 r=−

11

-26 p-refinement

h-refinement

d = 11 ln(N)

-36 1

3

5

7

Fig. 3. Convergence under p- and h-refinement.

11 5.3

22 10.1

44 16.2

3 0.9

7 4.0

11 7.6

h-refinement Degree (d) r a

Determined by linear fitting.

In a cylindrical co-ordinate system, the boundary conditions are adiabatic at the end surfaces z ¼ h, and convection and radiation conditions are applied on the lateral surface, r ¼ L. The problem is solved with two uniform hybrid elements (with lengths 0.5L) and a fifth-degree approximation. Variation of the temperature in time at points placed on the axis and on the curved surface of the cylinder shown is Fig. 4 (solid lines) recovers the reference (finite difference) solution given in [12]. A non-dimensional time is used,

s ¼ ðktÞ=ðcL2 Þ

9

ð49Þ

with k=L ¼ 2, k=ðT 30 LÞ ¼ 1:6  108 and smax ¼ 3. The convection coefficient is defined in terms of the Biot number, Bi ¼ hc L=k. The finite element solution reported in [2], which does not include the axial result for test Bi = 2.0, is obtained with a biased mesh of 12 elements (their lengths range from 0.1L to 0.025L). 7.5. Application to singular heat flow

Hs ðr; /Þ ¼ r n expðin /Þ in V s

d=7

ln (ε E ) ln (N )

Mesh r

The main motivation of one of the earliest applications of hybrid elements to transient heat transfer problems [5] was the enrichment of meshes of conform elements with hybrid elements embodying the singularity of the heat flow at the tip of cracked media. It consists, essentially, in enriching the temperature approximation basis (19) with functions,

k = 10 1

d =3

p-refinement

ð50Þ

with the origin of the polar co-ordinate system ðr; /Þ placed at the crack tip (see Fig. 5). This harmonic function, where i is the imaginary unit, is contained in the analytic solution of the linear, steadystate form of Eq. (4). Singular heat flow fields are usually modelled [5] with n ¼ 12 in Eq. (50), being the functions typically implemented in all elements of the mesh that contain the crack tip. A similar technique has been applied to the modelling of elastic fracture with hybrid-Trefftz ele-

20

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

= (T Ta) / (T0 Ta)

= (T Ta) / (T0 Ta)

= 1.0

1.0

= 1.0

1.0

0.8

0.8

Bi = 1

( ) Conform

0.6

Bi = 1

Surface (r = L): Bi = 1; 2; 3.33; 10

0.6 3.33

Axis (r = 0): Bi = 1; 2; 3.33; 10

0.4

( ) Conform

0.4

0.2

0.2

Bi = 10

Bi = 10

0.0 -3.0

log( ) -2.0

-1.0

0.0

1.0

log( )

0.0 -3.0

-2.0

-1.0

1.0

0.0

Fig. 4. Solid cylinder subjected to convection and radiation.

y

V2

h

r

φ

V1

h

x a

L−a

Fig. 5. Emery’s mode II problem.

ments, e.g. [13]. The increase in the number of degrees-of-freedom is marginal and fully compensates the gains in speed of convergence. In the context of hybrid elements, the main disadvantage is the increased difficulty of enforcing the inter-element thermal continuity condition (15) in the affected elements, as their bases now include temperature approximation functions with high gradients and heat flow fields that tend to infinity. In order to overcome this difficulty, the alternative technique reported here consists in implementing one set of functions (50) per each crack tip, independently of the number of elements that contain the singularity sources. To this effect, the support of those functions is defined as the domain in analysis, V. A subdomain defined by the combination of the elements of the mesh relevant to the analysis can also be used. For instance, the plate shown in Fig. 5 is solved here with the two-element mesh there identified. While the polynomial temperature approximation (18) is implemented in each element, with supports V 1 and V 2 , the domain of the plate is defined as the support of the local enrichment basis defined by Eq. (50), V s ¼ V ¼ V 1 [ V 2 . Moreover, the singularity condition is extended into form 12 6 n < 1 to assess the effect on convergence of lowerorder singular terms, characteristic of nonlinear applications. The implementation of this technique is illustrated with the problem defined in Fig. 5 assuming L ¼ h ¼ 2a, as taken from [5,14]. The initial temperature is null and adiabatic conditions are applied on edges x ¼ L, y ¼ 0 and y ¼ 2h. Dirichlet conditions are applied on edge x ¼ 0, where T ¼ T for y < h and T ¼ þT for y > h to induce a mode II singularity.

This problem is solved in [5] with a (biased) mesh of 36 conform (8-node) elements and 8 singularity-enriched hybrid elements, each with 16 temperature degrees-of-freedom, besides those associated with the heat flux approximation. The 8 (enriched) hybrid elements form a small disk centred on the crack-tip, with radius r  0.025L. The results presented in Figs. 6–8, where h ¼ T=T and the time scaling (49) is used, are obtained with polynomial approximations of degree five in each element and on each Dirichlet boundary (N ¼ 60, with 42 degrees-of-freedom for temperature and 18 for heat flux). The basis is enriched by the pairs of temperature approximation modes (50) associated with each value of the singularity order, n, being tested, that is: 12, 23 and 56  1. The temperature distributions on the edges of the crack presented in Figs. 6 and 7(a) illustrate the improvement in convergence that is attained simply by adding the singular solution: they show the solution obtained with the two-element mesh using only the (fifth-degree) polynomial approximation and the solution obtained enriching the same polynomial basis with the stronger singular mode, n ¼ 12. The contribution of the lower-order singular terms is weaker. Besides the improvement of the temperature estimate on the edges of the crack shown in Fig. 6, they contribute mainly in the earlier stages of the analysis to improve the modelling of the adiabatic boundary conditions on the edges of the crack, as shown in Fig. 7 (b). It is noted that these conditions are not explicitly enforced in the solving system (40). The temperature distributions presented in Fig. 8 recover the patterns and the values reported in [5]. They are obtained with the three sets of singular heat flow modes being

1.0

Time : τ = 1

θ

0.6

Edge y = h+ 0.2 -0.2

Edge y = h− -0.6 -1.0 0.0

0.1

0.2

0.3

0.4

x /L 0.5

Polynomial basis

n = 1/2

n = 1/2; 2/3

n = 1/2; 2/3; 5/6

Fig. 6. Temperature on crack edges at time s = 1.

21

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

0.8 0.6

0.8

0.4

1.0

1.0

0.4

0.2

1.0

0

0.2

0.2

0.8

0.4

1.0

0.6

0.8

1.0

0

1.0

0.2 1 2

Mode n

0.8 )

(a) Higher singular mode (

0

1.0

0.6

With mode n

Polynomial basis

1.0

0.2

0

0.4

0.2

0.2

0.6

0.2 1 2

With all modes

0.1 )

(b) Lower singular modes (

Fig. 7. Effect of singular modes in temperature distribution.

0.4

0.8 0.6

0.8

1.0

0.4

0

0

0

0.2

0.2

0.2

1.0

0.4

0.2

0.4

1.0

0.4

0.8 0.6

0.2

1.0

0.2

1.0

0.4

0.4

1.0

0.2

Time:

0.6

0.2

0

1.0

0.8

0.4

0.2

1.0

0.6

0.8

0.4

0.6

0.6

0.8

0.6

1.0

Fig. 8. Temperature distributions for Emery’s mode II problem.

tested, n1 ¼ 12, n2 ¼ 23 and n3 ¼ 56. The effect of the latter singular solution is marginal, as expected. 7.6. Application to cement hydration The simulation of cement hydration in a concrete block cast on a rock foundation is presented in [8] assuming that the moisture necessary to feed the chemical reactions is always available throughout the hydration process. This problem is reanalysed here assuming that moisture varies in space and in time, depending on the level of hydration and the effect of the boundary conditions of the problem (mainly, surface drying). As shown in [10], simulation of this hygro-thermo-chemical problem couples the heat transfer problem (4)–(7) with a moisture transport problem with a similar state equation,

$T cm km $h þ Q_ m ¼ cm h_ in V e

ð51Þ

where now h represents the relative humidity, cm is the moisture capacity, km the diffusivity matrix and Q_ the rate of moisture consumption. They depend on the chemical reaction laws described in terms of degrees of hydration, silica fume reaction and silicate polymerization, e.g. [15], and are defined in [10]. Boundary and initial conditions (5)–(7) can be similarly set. Solving systems similar to Eqs. (40) and (41), defined in detail in [10], are obtained for the moisture transport problem, as approximations (18) and (24) are extended to the relative humidity and moisture flux fields using Legendre bases (19) and (25), with dimensions still defined by Eq. (20). For practical temperature field simulations in concrete structures at early ages, conductivity and specific heat are often assumed as constant and isotropic during the cement hydration

process. Moreover, variation of the ratio between the concrete and air temperatures is relatively small, meaning that condition (12) is usually taken as linear, with a convection-radiation coefficient that varies with the wind speed. Consequently, the heat source is the only nonlinear and coupled term in the heat transfer problem [8]. The levels of nonlinearity and coupling are substantially stronger in the moisture transport problem. Permeability is a nonlinear function of temperature and relative humidity, the convective coefficients in the Neumann conditions depend on moisture capacity, defined as the slope of the sorption/desorption isotherms, which depend on the degrees of chemical reaction and relative humidity. Geometry and dimensions of the concrete block (subject to a hygro-thermo-chemical analysis) and of the rock foundation (subject to a thermal analysis) are defined in Fig. 9, under the assumption that symmetry conditions hold on surfaces x ¼ 0 and y ¼ 0. The cement hydration function used in the analysis is presented in Fig. 10, and the additional supporting data, taken from [10,15,16], is summarized in the Appendix. The boundary conditions of the heat transfer and moisture transport problems are defined in Fig. 11. A convective coefficient hcr ¼ 5 W/m2 K is used to simulate the formwork placed (throughout the test) on the lateral surfaces of concrete block, while a convection-radiation condition is applied to the top surfaces of both concrete and rock bodies using a combined coefficient hcr ¼ 9:53 W/m2 K, depending on the range of air temperature, wind velocity (1.0 m/s) and surface emissivity (0.9) [8]. Adiabatic conditions are assumed on the symmetry planes and heat transfer between the concrete block and the rock foundation is modelled with conductance hcd ¼ 50 W/m2 K.

22

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

with the rock foundation. The initial values for temperature, relative humidity and degree of cement hydration in the concrete block are 15°C, 100% and 103, respectively. The initial temperature in the rock foundation is 9.5°C. The variation in time of the air temperature in the exposed surfaces of the concrete block and rock foundation is presented in Fig. 12. The problem is analysed using a single element for the concrete block and four elements to discretize the rock foundation, as shown in Fig. 9. The temperature and relative humidity fields are approximated with bases of degree d ¼ 7, sufficient to capture the boundary layer effect associated with surface drying (convergence has been checked using ninth and eleventh degree approximations). The degree of approximation of the heat flux is d ¼ 4. As the Dirichlet boundary of the moisture transport problem is empty, the moist flux field is not independently approximated. The dimension o system (41) is N ¼ 660 for the heat transfer problem and N ¼ 120 for the moisture transport problem. The problem is analysed using a single element for the concrete block and four elements to discretize the rock foundation, as shown in Fig. 9. The temperature and relative humidity fields are approximated with bases of degree d ¼ 7, sufficient to capture the boundary layer effect associated with surface drying (convergence has been checked using ninth and eleventh degree approximations). The degree of approximation of the heat flux is d ¼ 4. As the Dirichlet boundary of the moisture transport problem is empty, the moist flux field is not independently approximated. The dimension of system (41) is N ¼ 660 for the heat transfer problem and N ¼ 120 for the moisture transport problem.

z

2 Concrete 4 y Rock

1 4

x

4m

1

Fig. 9. Concrete block cast on a rock foundation.

fc 1.0 0.8 0.6 0.4

Ta (oC)

0.2

15

αc

0.0

0.0

0.2

0.4

0.6

0.8

12

1.0

9

Fig. 10. Cement hydration function.

6 The concrete surfaces are exposed to a constant air relative humidity (0.5) and the moisture transport coefficient used to sim-

3

ulate the drying process is hcr ¼ 3:5  109 cm . The same value has been assumed for the surfaces enclosed by formwork because no values could be found in the literature. This option seemed more realistic than assuming sealed conditions. They are applied on the symmetry planes and in the interface of the concrete block

0 0

6

18

z Convection-radiation

Moisture convection

Sealed

Convection

z=5 Convection-radiation Adiabatic

Adiabatic

z=4 Conductance

Sealed

z=0 Adiabatic

(a) Heat transfer problem

24

Fig. 12. Variation of air temperature.

z z=6

12

(b) Moisture transport problem

Fig. 11. Boundary conditions in the cement hydration test.

Day 30

23

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

50

T (oC)

0.8

z=5

30

z=4

20

z=6

z=5

z=6 z=5 z=4: in concrete z=4: in rock z=0 Ta

40

αc z=4

0.6

z=6

0.4

z=0 10

0.2

Day

0 0

6

12

18

24

30

0.0 0

6

12

18

24

Fig. 13. Temperature at control points. Fig. 16. Hydration degree at control points.

12h

24h

60h

120h

Fig. 14. Evolution of the temperature field (back view) [45 °C

24h

60h

5 °C].

120h

Fig. 15. Evolution of degree of cement hydration (back view) [0.75

672h 0.05].

30

Day

24

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

1.00

θ

0.99

z=4

0.98

z=5

0.97 0.96 0.95

z=6

Day

0.94 0

6

12

18

24

30

Fig. 17. Relative humidity at control points.

Variation of the temperature at the control points indicated in Fig. 11(a) is shown in Fig. 13. As expected, the peak temperature of 45.4 °C (at 47 h) develops at the centre of the concrete bock (z ¼ 5 m) where most of heat induced by cement hydration is released. Variation of the temperature at the top of the hydrating concrete block (z ¼ 6 m) is strongly affected by the air temperature, while temperature at the base contacting the rock foundation remains almost unchanged (about 9.5 °C). The effect of cement hydration in the temperature of the rock foundation is visible (z ¼ 4 m). As time evolves, temperature at the top surface of the concrete block converges to the air temperature, and the temperature jump in the rock-concrete interface disappears. These effects

Back view

and the gradual cooling of the concrete block are shown in Fig. 14. This 3D solution recovers the patterns reported in [17,18] for similar tests. Heat transfer between concrete and rock modelled by the conductance condition is well captured. The evolution of the cement hydration degree presented in Figs. 15 and 16 shows that most of the heat of hydration is generated in the first six days, a6day  0:70, about 85% of its limit value, c a1 c  0:817, and confirms that hydration develops faster in the core of the concrete block. Air cooling delays significantly the hydration process, which is yet to be uniformly completed after the standard 28-day (672h) period. The results confirm that the thermo-chemical component of cement hydration is smooth in behaviour and not much demanding in terms of numerical modelling. The coupled moisture transport is distinct, as it is typically associated with edge effects caused by drying. The solutions obtained at the test control points and on the surfaces (of the symmetric simplification) of the concrete block are shown in Figs. 17 and 18. The development of relative humidity at the control points in Fig. 17 reflect different mechanisms of loss of moisture. At centre point of the concrete block (z ¼ 5 m) the drop in relative humidity during early-ages, say 6 days, is mainly caused by cement hydration. This effect is weaker at the rock interface point (z ¼ 4 m) due to the delay in hydration, and because no moisture exchange is assumed between the concrete block and the foundation. Conversely, at the point on the top surface (z ¼ 6 m), directly exposed to the air relative humidity, the moisture loss is mainly caused by drying. The results presented in Fig. 17 also confirm that the drop in relative humidity caused by cement hydration is relatively small [19] in normal strength concrete, being thus negligible in practice.

10 days

28 days

60 days

90 days Front view

Fig. 18. Evolution of relative humidity [1.0

Back view

Front view 0.78].

J.A. Teixeira de Freitas et al. / Computers and Structures 182 (2017) 14–25

The results presented in Fig. 18 show that the boundary layer caused by drying is captured using a single finite element. 8. Closure A hybrid formulation of the finite element method is applied to the analysis of nonlinear, transient problems. The approximation of geometry is independent of the one for the state variables. Consequently, the formulation can be implemented on meshes defined by alternative methods for domain decomposition and geometry mapping. As the number of elements is basically dictated by an adequate description of the geometry of the domain, thermal analyses can be performed using coarse meshes of high-degree elements. To enhance numerical stability under high-order p-refinement, naturally hierarchical, orthogonal functions are used to set up the domain and boundary approximations. Modelling is flexible because these approximations are independent and, for each field, different degrees can be used in each spatial direction. Moreover, it leads to solving systems that are sparse and well suited to adaptive refinement and parallel processing. Numerical validation shows that two main properties of hybrid elements extensively illustrated in Solid Mechanics applications are preserved in transient heat transfer problems, namely high rates of convergence, in particular under p-refinement, and a low level of sensitivity to gross mesh distortion. They also illustrate the modelling flexibility offered by hybrid formulations of the finite element method in two basic aspects, namely the implementation of approximation bases with overlapping supports and their enrichment to model local effects that may hinder convergence, typically the presence of high-gradient, singular or quasi-singular fields. It is shown in [8–10] that, in terms of relative performance, the hybrid element compares well and recovers the levels of accuracy reported for the conform finite element. Acknowledgements This research work has been supported by the Portuguese Foundation for Science and Technology (FCT) through Research Units CERIS (IST-ID) and CONSTRUCT (University of Porto). The authors also gratefully acknowledge the advice and support given by colleague C. Tiago. Appendix A The following data is used in the analysis of the hydration of the concrete block (hygro-thermo-chemical behaviour) cast on a rock foundation (thermal behaviour): (a) Cement hydration data: activation energy, Eac ¼ 50 kJ/mol; latent heat of hydration, Q 1 c ¼ 727 kJ/kg; ratio of the

25

maximum value of the heat production rate to the latent heat of hydration, Ac ¼ 0:2991  104 s1; cement content, cc ¼ 220 kg/m3; water-to-cement ratio, w=c ¼ 0:739; asymptotic degree of cement hydration, a1 c ¼ 0:818. (b) Concrete heat transfer data: conductivity, k ¼ 2:6 W/m K; specific heat c ¼ 2:4 MJ/m3 K. (c) Rock heat transfer data: conductivity, k ¼ 2:79 W/m K; specific heat c ¼ 2:04 MJ/m3 K. (d) Moisture transport data: Arrhenius energy, Ead ¼ 22:448 kJ/mol; isothermal diffusivity coefficients [16], aH ¼ 0:05, hH ¼ 0:8, n ¼ 15, h0 ¼ 22:85 °C. The data on water content is defined as in [10,15], with c g c ¼ 1:50 and kv g ¼ 0:225. References [1] Zienkiewicz OC, Cheung YK. Finite elements in the solution of field problems. Engineer 1965;220:507–10. [2] Bathe KJ, Khoshgoftaar MR. Finite element formulation and solution of nonlinear heat transfer. Nucl Eng Des 1979;51:389–401. [3] Fraeijs de Veubeke BM, Hogge MA. Dual analysis for heat conduction problems by finite elements. Int J Numer Meth Eng 1972;5:65–82. [4] Jirousek J, Qin QH. Application of hybrid-Trefftz element approach to transient heat conduction analysis. Comput Struct 1996;58(1):195–201. [5] Chen W-H, Ting K. Hybrid finite element analysis of transient thermoelastic fracture problems subjected to general heat transfer conditions. Comput Mech 1989;4:1–10. [6] Chen HT, Chang SM. Application of the hybrid method to inverse heat conduction problems. Int J Heat Mass Transfer 1990;33(4):621–8. [7] Farhloul M, Mounim AS. A mixed-hybrid finite element method for convection-diffusion problems. Appl Math Comput 2005;171(2):1037–47. [8] Freitas JAT, Cuong PT, Faria R, Azenha M. Modelling of cement hydration in concrete structures with hybrid finite elements. Finite Elem Anal Des 2013;77:16–30. [9] Freitas JAT, López C, Cuong PT, Faria R. Hybrid finite element thermal modelling of fire protected structural elements strengthened with CFRP laminates. Compos Struct 2014;113:396–402. [10] Freitas JAT, Cuong PT, Faria R. Modelling of cement hydration in high performance concrete structures with hybrid finite elements. Int J Numer Meth Eng 2015;103:364–90. [11] Freitas JAT, Moldovan ID, Toma M. Mixed and hybrid stress elements for biphasic media. Comput Struct 2010;88(23–24):1286–99. [12] Sucec J, Kumar A. Transient cooling of a solid cylinder by combined convection and radiation at its surface. Int J Numer Meth Eng 1973;6:297–303. [13] Freitas JAT, Ji Z-Y. Hybrid-Trefftz equilibrium model for crack problems. Int J Numer Meth Eng 1996;39(4):569–84. [14] Emery AF, Neighbors PK, Kobayashi AS, Love WJ. Stress intensity factors in edge-cracked plates subjected to transient thermal singularities. J Pressure Vessel Technol 1977;5:100–5. [15] Di Luzio G, Cusatis G. Hygro-thermo-chemical modeling of high performance concrete. I: Theory. II: Numerical implementation, calibration, and validation. Cement Concr Compos 2009;31:301–24. [16] CEB-FIP Model Code 1990. London: Thomas Telford; 1993. [17] Lee Y, Yi ST, Kim MS, Kim JK. Evaluation of a basic creep model with respect to autogenous shrinkage. Cem Concr Res 2006;36(7):1268–78. [18] Lee Y, Kim JK. Numerical analysis of the early age behavior of concrete structures with a hydration based microplane model. Comput Struct 2009;87 (17–18):1085–101. [19] Bazˇant ZP, Najjar LJ. Nonlinear water diffusion in nonsaturated concrete. Matériaux et Construction 1972;5(1):3–20.