Numerical analysis of coupled heat and moisture transport in masonry

Computers and Mathematics with Applications 74 (2017) 229–248

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Numerical analysis of coupled heat and moisture transport in masonry Tomáš Krejčí ∗ , Jaroslav Kruis, Michal Šejnoha, Tomáš Koudelka Faculty of Civil Engineering, Czech Technical University in Prague, Czech Republic

article

info

Article history: Available online 5 April 2017 Keywords: Heat and moisture transfer Masonry Finite element method Parallel computing Homogenization Processor farm

abstract Coupled analysis of heat and moisture transport in real world masonry structures deserves a special attention because the spatial discretization by the finite element method leads usually to large number of degrees of freedom. Thin mortar layers and large bricks or stones have very different material properties and the finite element mesh has to be able to describe correct temperature and moisture fields in mortar and in its vicinity in the blocks. This paper describes two possible solutions of such problems. The first solution is based on the domain decomposition method executed on parallel computers, where the Schur complement method is used with respect to non-symmetric systems of algebraic equations. The second alternative method is the application of a multi-scale approach in connection with a processor farming method, where the whole structure is described by a reasonably coarse finite element mesh, called the macro-scale problem, and the material parameters are obtained from the lower-level problems, called the meso-scale problem, by a homogenization procedure. In this procedure, the macro-problem is assigned to the master processor and the meso-scale problems belong to the slave processors in the processor farm. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Historical masonry structures such as, e.g., masonry arch bridges are seriously influenced by climatic loading. In many cases, it is responsible for the nucleation and further development of cracks and damage. Restoration and reconstruction of such damaged structures need complex analyses which are able to combine experimental work and numerical simulation [1]. A thermo-hygro-mechanical simulation of the response of masonry structures to climatic loadings can be shown as a suitable example of such complex numerical analysis. Numerical modelling of coupled hygro-thermal problems in masonry by the finite element method is hardly solvable on single processor computers. The huge number of degrees of freedom is usually obtained in case of detailed finite element discretization of thin mortar layers between bricks and stones. Miehe and Bayreuther [2] suggest to solve such systems by multigrid method but very large systems cannot be executed on a single processor computers. One of the possible solutions is an application of a domain decomposition method [3–5] in connection with parallel computers. There are two groups of domain decomposition methods: non-overlapping and overlapping methods. In the case of overlapping method, suitable overlap is prescribed between two adjacent subdomains. Overlapping methods are represented by the Schwarz method. In non-overlapping methods, subdomains share only nodes, edges or surfaces. The Schur complement method is very robust tool for solution of symmetric as well as non-symmetric systems. Another example of non-overlapping method is FETI (finite element tearing and interconnecting) method. Since 1991 when the method

∗

Corresponding author. E-mail address: [email protected] (T. Krejčí).

http://dx.doi.org/10.1016/j.camwa.2017.03.022 0898-1221/© 2017 Elsevier Ltd. All rights reserved.

230

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

was introduced [6], several variants were formulated. The original method is denoted FETI-1. In 1998, so called FETI-2 was introduced in [7,8] and FETI-DP was published in [9]. FETI method with all floating subdomains was described in [10]. The other possible solution of the given difficulties of masonry and heterogeneous materials and structures can be seen in an application of parallel computing and a multi-scale approach. When dealing with these problems, the application of homogenization techniques is inevitable [11]. Solving a set of problem equations on a meso-scale (a composition of stone blocks and mortar) provides with up-scaled macroscopic equations. They include a number of effective (macroscopic) transport parameters, which are necessary for a detailed analysis of the state of a structure as a whole [12]. Such a multi-scale approach to coupled heat and moisture transfer is based on a clear description of transport phenomena. An extensive review of this topic can be found in [13]. Averaging theories (a micromechanics-based approach), formulated, e. g., in [14,15], are considered as a counterpart to phenomenological models (a macromechanics-based approach), see [16,17]. When calculating the heat and moisture transfer in building materials, phenomenological models are still preferred to averaging ones. Both approaches are explained in detail in [18]. A phenomenological model proposed by Künzel and Kiessl [19] was chosen for studying transport processes in masonry structures due to its ability to describe all substantial phenomena of heat and moisture transport. While models for transport processes have been developed during several decades, the computational methods for multiscale modelling of these processes in masonry on meso and macro scales have emerged only recently. An original approach to homogenization of transient heat transfer for some composite materials is proposed in [20]. The complex multi-scale analysis for pure heat transfer in heterogeneous solids is offered in [11], where the authors established a macro to micro transition in terms of the applied boundary conditions and likewise a micro to macro transition formulated in the form of consistent averaging relations. The similar study with applications to textile composites can be found, e.g. in [21]. Homogenization strategies for coupled heat and mass transfer used in this paper are described in detail in [22,12,23]. Therein, governing equations of the coupled heat and moisture transport are derived in the framework of coupled two-scale analysis (the firstorder homogenization approach) in conjunction with the finite element method. The homogenized macro-scale fields are found from the solution of a certain sub-scale (meso-scale) problem performed on a representative volume element (RVE). This paper presents two computational methods based on parallel computing for effective simulation of coupled heat and moisture transfer in masonry. The first approach is the domain decomposition method, where the original domain is split into several smaller subdomains which are assigned to processors of a parallel computer. The second approach is the processor farming method in connection with a multi-scale analysis. In this method, each macro-scopic integration point or each finite element is connected with a certain meso-scopic problem represented by an appropriate representative volume element or a periodic unit cell (PUC). The solution of a meso-scale problem then provides effective parameters needed on the macro-scale. Such an analysis is suitable for parallel computing because the meso-scale problems can be distributed among many processors. In this regard, the master–slave strategy can be efficiently exploited. The processor farming method based on multi-scale analysis differs from classical parallel computing methods which come out from the domain decomposition. The macro-problem is assigned to the master processor while the solution at the meso-level is carried out on slave processors. At each time step the current temperature and moisture together with the increments of their gradients at a given macro-scopic integration point are passed to the slave processor (imposed onto the representative volume element), which, upon completing the small scale analysis, sends the homogenized data (effective conductivities, averaged storage terms and fluxes) back to the master processor. It means that the conductivity and capacity matrix of an element of the macro-scale problem are assembled as a result of the homogenization on meso-scale. The presented parallel methods have been implemented into a parallel version of SIFEL (Simple Finite Elements) code [24] with the distributed memory scheme and the MPI communication library. The code uses the master and slaves concept, where the master processor manages communication among all processors and it also controls the computation. The paper consists of seven sections. After the first Introduction, the basic approach based on Künzel and Kiessl coupled heat and moisture transfer model is summarized in Section 2. The basic principles of domain decomposition together with Schur complement method are discussed in Section 3. The first order homogenization procedure is described in Section 4. Section 5 deals with the processor farming method and results of two numerical simulations of heat and moisture transfer are illustrated in Section 6. The first benchmark, the 2D analysis of a brick wall, shows the application of both methods presented. The second example is the analysis of Charles bridge in Prague which illustrates the application of processor farming method to a real world masonry structure. Finally, the benefits of used methods and their comparison are discussed in Section 7. Customary matrix notation is used throughout the text. Matrices are denoted by uppercase boldface italic letters, e.g., D, P, etc. Conversely, lowercase boldface italic letters stand for vectors g , j, etc. ∇· is the divergence operator, while ∇ stands for the gradient operator. 2. Description of coupled heat and moisture transfer and numerical solution For the description of transport phenomena, phenomenological modelling has still been preferred to micromechanicsbased (averaging) description outlined, e.g., in [18]. In this paper, we focus just on one particular model implemented in the SIFEL computer code by which all the results discussed in this paper were obtained. It is the diffusion model proposed by Künzel and Kiessl [19]. It is much simpler than micromechanics-based models and it describes all substantial phenomena very well. This phenomenological model introduces two unknowns, relative humidity ϕ [-] and temperature T [K] in a material point. The practical usage of the main advantages of this model can be seen, e.g., in Refs. [25,26], where building

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

231

structures are simulated under common climatic conditions with the help of material properties measured experimentally. The model supposes the moisture transport mechanisms relevant to numerical analysis in the field of building physics are just water vapour diffusion and liquid transport. Vapour diffusion is most important for large pores, whereas liquid transport takes place on pore surfaces and in small capillaries. Vapour diffusion in porous media is described in the model by the Fick’s diffusion and effusion in the form jv = −δp ∇ p = −

δ ∇ p, µ

(1)

where jv [kg m2 /s] is the water vapour flux, δp [kg m s−1 Pa−1 ] is the vapour permeability of the porous material, p denotes vapour pressure [Pa], µ is the vapour diffusion resistance number and δ [kg m s−1 Pa−1 ] is the vapour diffusion coefficient in air. The liquid transport mechanism includes liquid flow in the absorbed layer (surface diffusion) and in the water filled capillaries (capillary transport). The driving potential in both cases is capillary pressure (suction stress) or relative humidity ϕ . The flux of liquid water is described in the form jw = −Dϕ ∇ϕ,

(2)

where jw [kg m /s] is the liquid water flux, the liquid conductivity Dϕ = 2

Dw ddw ϕ −3

−1

[kg m s ] is the product of the liquid

diffusivity Dw [m2 s−1 ] and the derivative of water retention function. w [kg m ] is the water content of the material. The heat flux is proportional to the thermal conductivity of the moist porous material and the temperature gradient (Fourier’s law) q = −λ∇ T ,

(3)

where q [W/m2 ] is the heat flux and λ [W m−1 K−1 ] denotes the thermal conductivity of the moist material. The heat and moisture balance equations are closely coupled because the moisture content depends on the total enthalpy and thermal conductivity while the temperature depends on the moisture flow. The resulting set of differential equations for the description of simultaneous heat and moisture transfer, expressed in terms of temperature, T , and relative humidity, ϕ , have the form of partial differential equations defined on a domain Ω dw ∂ϕ



ρC +

dϕ ∂ t ∂ Hw  ∂ T

∂T

∂t

  = ∇ · Dϕ ∇ϕ + δp ∇(ϕ psat ) ,

(4)

    = ∇ · λ ∇ T + hv ∇ · δp ∇(ϕ psat ) ,

(5)

where ddw [kg/m3 ] is the derivative of moisture storage function, Hw [J m−3 ] denotes the enthalpy of the material moisture, ϕ

hv [J kg−1 ] is the evaporation enthalpy of the water, psat [Pa] is the water vapour saturation pressure, ρ [kg m−3 ] is the material density, C [J kg−1 K−1 ] is the specific heat capacity and t [s] denotes time. Boundary of the domain Ω is split into parts ΓT , Γϕ , ΓqpT , Γjpϕ , ΓqcT and Γjc ϕ which are disjoint and their union is the whole boundary Γ . The system of equations (4) and (5) is accompanied with three types of boundary conditions:

• Dirichlet boundary conditions T (x, t ) = T (x, t ), x ∈ ΓT ϕ(x, t ) = ϕ(x, t ), x ∈ Γϕ

(6) (7)

• Neumann boundary conditions dT

= q(x, t ) = q(x, t ),

x ∈ ΓqpT ,

(8)

= j (x, t ) = j (x, t ), dn • Cauchy boundary conditions

x ∈ Γjpϕ ,

(9)

−λ

dn dϕ

−D ϕ

q(x, t ) = α(T (x, t ) − T∞ (x, t )),

x ∈ ΓqcT ,

(10)

j (x, t ) = β(p(x, t ) − p∞ (x, t )),

x ∈ Γjc ϕ ,

(11)

where T (x, t ) is the prescribed temperature, ϕ(x, t ) is the prescribed relative humidity, q(x, t ) is the prescribed heat flux, j (x, t ) is the prescribed moisture flux, α [W m−2 K−1 ] and β [kg s−1 Pa−1 ] are the heat and mass transfer coefficients, T∞ is the ambient temperature and p∞ is the ambient water vapour pressure. The finite element method is used for spatial discretization of the partial differential equations (4) and (5). The weighted residual statement is applied to the mass balance equation assuming δ T = 0 on ΓT and δϕ = 0 on Γϕ

 Ω

δϕ

 dw ∂ϕ dϕ ∂ t

  − ∇ · Dϕ ∇ϕ + δp ∇(ϕ psat ) dΩ = 0

(12)

232

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

and also to the energy balance equation

 Ω

δT



∂ Hw ∂T

ρC +

   ∂T − ∇ · (λ ∇ T ) − hv ∇ · δp ∇(ϕ psat ) dΩ = 0. ∂t



(13)

Applying Green’s theorem the weak formulation for the mass transfer yields

 Ω

δϕ

 dw ∂ϕ  dϕ ∂ t





− Γjpϕ

δϕ Dw



dΩ + dw dϕ

dw



Ω

∇δϕ Dw

+ δp psat

 ∂ϕ ∂n

dϕ

   dpsat  ∇δϕ δp ϕ ∇ T dΩ + δp psat ∇ϕ dΩ +

dΓ −

Ω



dT

dpsat  ∂ T δϕ δp ϕ dΓ = 0 dT ∂n



ΓqpT

(14)

and the weak formulation for the heat transfer



   ∂ Hw  ∂ T dpsat  ∇ T dΩ δT ρ C + dΩ + ∇δ T λ + hv δp ϕ ∂T ∂t dT Ω Ω       ∂ϕ ∇δ T hv δp psat ∇ϕ dΩ − δ T hv δp psat + dΓ ∂n Ω Γjpϕ   dpsat  ∂ T dΓ = 0. − δ T λ + hv δp ϕ dT ∂n ΓqpT

(15)

In the finite element method, the temperature T and relative humidity ϕ are approximated in the form T = N (x)rT ,

ϕ = N (x)rϕ ,

(16)

and the gradients of temperature and relative humidity are also needed in the form

∇ T = B(x)rT ,

∇ϕ = B(x)rϕ ,

(17)

where N (x) denotes the matrix of approximation functions, B(x) is the matrix of their derivatives, rT denotes the vector of nodal temperatures and rϕ denotes the vector of nodal relative humidities. Using approximation Eqs. (16) and (17) in Eqs. (14) and (15), a set of the first order differential equations is obtained in the matrix form



Kϕϕ KT ϕ

Kϕ T KTT

  rϕ rT

 +

Cϕϕ CT ϕ

    r˙ϕ j = ϕ . r˙T qT

Cϕ T CTT

(18)

The matrices Kϕϕ , Kϕ T , KT ϕ and KTT create the conductivity matrix of the problem and they have the form

 Kϕϕ =

Ω

 KT ϕ =



BT Dϕϕ BdΩ ,

Kϕ T =



BT DT ϕ BdΩ ,

Ω

BT Dϕ T BdΩ ,

(19)

BT DTT BdΩ .

(20)

Ω

KTT =

Ω

The material conductivity matrices Dϕϕ , Dϕ T , DT ϕ and DTT are diagonal matrices and the diagonal entries are equal to following conductivities dϕϕ = Dw

dw dϕ

+ δp psat ,

dT ϕ = hv δp psat ,

dϕ T = δp ϕ

dpsat dT

,

dTT = λ + hv δp ϕ

dpsat dT

(21)

.

The matrix elements Cϕϕ , Cϕ T , CT ϕ and CTT in the capacity matrix of the problem have the form

 Cϕϕ =

Ω

 CT ϕ =

Ω

N T Hϕϕ N dΩ , N T H T ϕ N dΩ ,

 Cϕ T =

Ω

 CTT =

Ω

N T Hϕ T N dΩ ,

(22)

N T HTT N dΩ ,

(23)

with diagonal material capacity matrices Hϕϕ , Hϕ T , HT ϕ and HTT corresponding to the capacities hϕϕ =

dw dϕ

hT ϕ = 0,

,

hϕ T = 0 , hTT = ρ C +

∂ Hw . ∂T

(24)

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

233

The vectors jϕ and qT contain prescribed nodal fluxes and have the form

 jϕ =

Γj



N T jϕ dΓ ,

qT =

Γq

N T qT dΓ .

(25)

jϕ denotes the mass boundary flux and qT stands for the heat boundary flux. The system of differential equations (18) can be written more compactly in the form K (r )r + C (r )˙r = f .

(26)

The dependency of the conductivity and capacity matrices on the attained values of variables is explicitly denoted. Letters r and r˙ denote the nodal variables and their time derivatives, respectively. In the SIFEL code, the system of equations (26) is usually solved by the v -form of the generalized trapezoidal method [27] defined by the relationships rn+1 = rn + 1tvn+α ,

(27)

vn+α = (1 − α)vn + α vn+1 ,

(28)

where v denotes the first derivatives of nodal values with respect to time. The subscript n stands for the time step and it serves also as an index in the incremental method, called the outer iteration loop. It is assumed that all variables are known at the time tn and variables at the time tn+1 are searched. Substitution of expressions defined in Eqs. (27) and (28) to the system of differential equations (26) leads to relationship

(Cn + 1t α Kn ) vn+1 = fn+1 − Kn (rn + 1t (1 − α)vn ) ,

(29)

with the capacity Cn and the conductivity Kn matrices evaluated with the help of values rn . The system of algebraic equations (29) is generally non-linear and the Newton–Raphson method [28,29] has to be used at each time step. 3. Parallel computing It was already discussed in the Introduction that the sequential computation of the coupled heat and moisture transfer in masonry is very time and memory consuming. Even for relatively small segments of real geometry, the single processor computation is suffered by large memory requirements and the computational time is unacceptably long. Therefore, parallelization of the problem based on domain decomposition methods can be successfully used. The domain decomposition methods split the original domain solved into several smaller subdomains which are assigned to processors of a parallel computer. The continuity conditions on subdomain interfaces are prescribed. The domain decomposition method speeds up not only solution of the system of linear algebraic equations but also evaluation of constitutive equations and assembling of the system matrices. The last mentioned speed up of evaluation of constitutive equations is distinctly visible in the coupled problems, where the equations are quite complex. The presented parallel algorithm is based on the distributed memory scheme and MPI communication library. It uses the master and slaves concept where the master processor manages communication among all processors and it also controls the computation. 3.1. Schur complement method There are several approaches for the parallel solution of the algebraic equation systems (Schur complement method, FETI, FETI-DP, TFETI). The Schur complement method is convenient for coupled problems due to the non-linearity of the non-symmetric system of equations (29). The method is a non-overlapping domain decomposition which uses the original unknowns during the whole computation. It can be formulated as a special version of the Gaussian elimination method and can be built on the LDL T factorization. It means, nonlinear systems can be solved by modified Newton–Raphson method, which uses the initial Jacobi matrix, and factorization of the matrix is performed only once. In the domain decomposition method, the original domain is split into m non-overlapping subdomains. A special ordering of unknown results in the system of equations in the form

          

[ii]

[ib]   [i] 

0

K1

K1

[ii] [ii]

K3

..

. Km[ii]

K1

[bi]

K2

[i] 

f1

 [i]    d2   f 2      [i]  [ib] [i]   K3   d3      f3   .  =  . , ..      .   ..   ..      [i]  f  Km[ib]   d [i]  K2

.. . [bi]



[ib]   [i] 

K2 0

d1

[bi]

K3

...

[bi]

Km

K

[bb]

m

m

d [b]

f [b]

(30)

234

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

[i]

where K with the subscripts and superscripts denotes the submatrix, dj are the vectors of unknowns defined inside the jth subdomain, d [b] is the vector of unknowns on the subdomain interfaces, fj

[i]

are the vectors of right hand sides defined

inside the jth subdomain and f [b] is the vector of right hand side on the subdomain interfaces. The superscript i denotes an internal quantity while b stands for an interface one. The interface quantity can be understood as a boundary quantity on a subdomain. Clearly, ib and bi denotes coupling between internal and interface quantities. The number of equations collected [i] in the jth block is denoted by nj . The number of equations collected in the last block is denoted by n[b] . If the number of all rows in the system Eq. (30) is n, the following simple relationship is valid n = n[b] +

m 

nj . [i]

(31)

j =1

[ii]

The diagonal blocks can be written as Kj [i ] nj ×n[b]

[ib]

[i ]

∈ Rnj

[ii]

[i]

×nj

[i]

, while the vectors dj and fj

[i]

are from R

[i] nj

and the off-diagonal blocks

[ii]

Kj are from R . If all diagonal blocks K1 to Km are nonsingular, particular blocks of the vector of unknowns d can be expressed in the form [i]



[ii]

dj = Kj

−1 

fj

[i]

 − Kj[ib] d [b] .

(32)

After repeating this process for each row of Eq. (30) except the last one, the reduced system of equations becomes

 K

[bb]

−

m 

[bi]

Kj



[ii]

Kj

−1

j =1

 [ib]

Kj

d [b] = f [b] −

m 

[bi]

Kj



[ii]

Kj

−1

fj , [i]

(33)

j =1

[b]

[b ]

[b ]

where d [b] ∈ Rn is the vector of unknowns. The matrix of the system (33) is from the set Rn ×n . The original system of equations contains n unknowns while the reduced system (33) contains only n[b] unknowns. This reduction of unknowns is an important feature of the Schur complement method. When the vector d [b] is computed, all vectors of internal unknowns [i] dj are established from (32). More details about the method can be found in Ref. [5]. 3.2. Nonlinear system of equations Computer codes using finite element method are usually equipped with a module for solution of systems of linear algebraic equations. On parallel computers with distributed memory, the generalization for systems of nonlinear algebraic equations solved by the Newton–Raphson method requires an implementation of additional subroutines. Namely, the dot product of two vectors which are allocated on different processors and subroutine for interruption of the inner loop in the Newton–Raphson method running on different processors. The numerical tests showed that the master processor has to collect the contributions of the residual vector from slaves and it has to send a message about the interruption of the inner loop, see, e.g., [5]. The whole parallel algorithm in SIFEL computer code is described in Tables 1 and 2. 4. Homogenization and solution of multi-scale problem 4.1. First order homogenization of non-stationary coupled heat and moisture transport An alternative approach to domain decomposition method is a two-scale analysis with a processor farming method suggested by the authors. The two-scale approach with a macro-scale problem dealing with the structure and a mesoscale problem describing the type of masonry has been introduced in [22,12]. This formulation corresponds to the first order homogenization method, where only the function value and the gradient of the macro-level function are used at the meso-level. Governing equations of the coupled heat and moisture transport are derived in the framework of coupled two-scale analysis of finite element type. It is presumed that the homogenized macro-scale fields are found from the solution of a certain sub-scale (meso-scale) problem performed on a representative volume element which is represented by a periodic unit cell (PUC). Most often regular periodic unit cells have been assumed, e.g., in [30,23,21]. Irregular mesostructures have been approached either directly in the framework of stochastic finite elements [31] or indirectly by introducing the concept of statistically equivalent periodic unit cell (SEPUC) [32]. Such a unit cell is then constructed to match the real meso-structure, at least in a statistical sense, as close as possible. Examples of such SEPUCs for both regular as well as irregular masonry are plotted in Figs. 1 and 2. For non-mechanical fields, the path paved in [33] is followed. Therein, it has been advocated that for a finite size RVE the assumption of transient flow on both the macro- and the meso-scale introduces certain non-local, size dependent, terms in equations governing the macro-scopic response. To introduce this subject suppose that a local field a can be replaced by a spatially homogenized one ⟨a⟩ such that

 Ω

a dΩ ≈

 Ω

⟨a⟩ dΩ =

    1 a dΩ dΩ , Ω Ω Ω

(34)

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

235

Table 1 Parallel version of the Newton–Raphson iteration—outer loop.





Γ

    1 a dΓ dΓ , a dΓ ≈ ⟨a⟩ dΓ = Γ Γ Γ Γ

(35)

where Ω and Γ represent the internal and boundary parts of the SEPUC. Employing Eqs. (34) and (35) into weak formulations (14) and (15), the mass balance equation leads to the form

      dw   dw ∂ϕ  δϕ dΩ + ∇δϕ Dw + δp psat ∇ϕ dΩ  dϕ ∂ t  dϕ Ω Ω         ∂ϕ  dpsat  dw + ∇δϕ δp ϕ ∇ T dΩ − δϕ Dw + δp psat dΓ  dT dϕ ∂n  Ω Γjpϕ    dpsat  ∂ T  − δϕ δp ϕ dΓ = 0, dT ∂n  ΓqpT

(36)

236

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248 Table 2 Parallel version of the Newton–Raphson iteration—inner loop.

and subsequently the energy balance equation has the following form

     ∂ Hw  ∂ T  dpsat  δT ρ C + ∇ T dΩ dΩ + ∇δ T λ + hv δp ϕ  ∂T ∂t  dT Ω Ω          ∂ϕ  + dΓ ∇δ T hv δp psat ∇ϕ dΩ − δ T hv δp psat  ∂n  Ω Γjpϕ    dpsat  ∂ T  − δ T λ + hv δp ϕ dΓ = 0 . dT ∂n  ΓqpT

 

(37)

In the spirit of the first order homogenization it is assumed that the macro-scopic temperature and relative humidity vary only linearly over the SEPUC. This can be achieved by loading its boundary by the prescribed temperature T hom and relative humidity ϕ hom derived from the uniform macro-scopic temperature ∇ T and relative humidity ∇ϕ gradients, respectively. In such a case the local temperature and relative humidity inside the SEPUC admit the following decomposition T (x) = T (X 0 ) + ∇ T {x − X 0 } + θ ∗ (x) = T hom (x) + θ ∗ (x),

(38)

ϕ(x) = ϕ(X ) + ∇ϕ{x − X } + ϕ (x) = ϕ

(39)

0

0

∗

hom

(x) + ϕ (x), ∗

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

237

Fig. 1. 3D meso-structural mesh of regular bonding of masonry SEPUC1.

Fig. 2. 3D meso-structural mesh of irregular bonding of masonry SEPUC2.

where θ ∗ (x) and ϕ ∗ (x) are the fluctuations of local fields superimposed onto linearly varying quantities T hom (x) and ϕ hom (x). The temperature T (X0 ) and the moisture ϕ(X0 ) at the reference point X0 are introduced to link the local fields to their macro-scopic counterparts. For convenience the SEPUC is typically centred at X0 . Henceforth, the local fluctuations will be demanded to be periodic, i.e. the same values are enforced on the opposite sides of a rectangular SEPUC. Next, substituting Eqs. (38) and (39) into Eqs. (14) and (15) and collecting the terms corresponding to δ T hom , δϕ hom and δθ ∗ , δϕ ∗ split the weak forms of the original problems (14) and (15) into the homogenized (macro-scale) problems, for moisture transfer

     dw ∂ϕ   dw   + δp psat ∇ϕ dΩ δϕ hom dΩ + ∇δϕ hom Dw  dϕ ∂ t  dϕ Ω Ω        dw  ∂ϕ   dp sat + ∇δϕ hom δp ϕ ∇ T dΩ − δϕ hom Dw + δp psat dΓ  dT dϕ ∂n  Γjpϕ Ω    dpsat  ∂ T  − δϕ hom δp ϕ dΓ = 0 , dT ∂n  ΓqpT and for heat transfer

       dpsat  ∂ Hw  ∂ T  δ T hom ρ C + dΩ + ∇δ T hom λ + hv δp ϕ ∇ T dΩ  ∂T ∂t  dT Ω Ω

(40)

238

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

      + ∇δ T hom hv δp psat ∇ϕ dΩ − Ω

 −



 ∂ϕ    δ T hom hv δp psat dΓ ∂n  Γjpϕ

  dpsat  ∂ T  δ T hom λ + hv δp ϕ dΓ = 0, dT ∂n  ΓqpT

(41)

and the local sub-scale (meso-scale) problems, for moisture:

        dw dw ∂ϕ  δϕ ∗ dΩ + ∇δϕ ∗ Dw + δp psat ∇ϕ dΩ  dϕ ∂ t  dϕ Ω Ω          ∂ϕ  dp d w sat ∇δϕ ∗ δp ϕ δϕ ∗ Dw + ∇ T dΩ − dΓ + δp psat  dT dϕ ∂n  Ω Γjpϕ    dpsat  ∂ T  δϕ ∗ δp ϕ − dΓ dT ∂n  ΓqpT     dw 1 = + δp psat ∇ϕ d(∂ Ω )dΩ δϕ ∗ Dw j dϕ Ω |Ω | ∂ Ω    1 dpsat  + δϕ ∗ δp ϕ ∇ T d(∂ Ω )dΩ = 0, q dT Ω |Ω | ∂ Ω

(42)

and for heat:

     ∂ Hw  ∂ T  dpsat  dΩ + ∇ T dΩ δθ ∗ ρ C + ∇δθ ∗ λ + hv δp ϕ  ∂T ∂t  dT Ω Ω          ∂ϕ  + ∇δθ ∗ hv δp psat ∇ϕ dΩ − dΓ δθ ∗ hv δp psat  ∂n  Ω Γjpϕ    dpsat  ∂ T  − δθ ∗ λ + hv δp ϕ dΓ dT ∂n  ΓqpT     1 δθ ∗ hv δp psat ∇ϕ d(∂ Ω )dΩ = j Ω |Ω | ∂ Ω    1 dpsat  δθ ∗ λ + hv δp ϕ ∇ T d(∂ Ω )dΩ = 0, + q dT Ω |Ω | ∂ Ω

 

(43)

which are satisfied identically owing to the assumed periodic boundary conditions. Solving Eqs. (42) and (43) for the prescribed increments of ∇ T and ∇ϕ provides the instantaneous effective properties and storage terms that appear in the macro-scale equations (40) and (41), respectively. 4.2. Macroscopic conductivity and capacity matrices The evaluation of macro-scopic conductivity and capacity matrices in Eqs. (40) and (41) comes out from the solution of a simple steady state heat and moisture conduction problem using the finite element method. The stepping stone in the estimation of macro-scopic response is the Hill–Mandel lemma suggesting equality of the work of local fields averaged over the solution domain and the work of their macro-scopic counterparts [22]. Owing to the non-linear nature of underlying material models, the formulation must be written in an incremental form. It assumes that the heat and moisture balance is reached at the end of the ith step at the meso-structural level and the nodal non-equilibrated forces are equal to zero. The equilibrium state is then assumed at the end of the ith time step and a small increment of the macro-scopic temperature and moisture gradients d∇ T and d∇ϕ is considered resulting in an incremental change of local and macro-scopic fluxes

⟨−δ(∇ϕ)T dj ⟩ = −δ(∇ϕ)T dj M ,

(44)

⟨−δ(∇ T ) dq⟩ = −δ(∇ T ) dq ,

(45)

T

T

M

where the symbol ⟨⟩ represents the volume averaging and M denotes the macro-scopic level. The Dirichlet boundary conditions are admitted in the form of the prescribed gradients of macro-scopic temperature and moisture, and the terms on the right hand side would disappear (δ(∇ T ) = 0, δ(∇ϕ) = 0). For the calculation purposes these conditions are considered to determine the unknown fluctuation fields. The discretized forms of the local temperature and moisture gradients read

∇ T = ∇ T + BrT∗ ,

∇ϕ = ∇ϕ + Brϕ∗ ,

(46)

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

239

then

δ(∇ T ) = δ(∇ T ) + Bδ rT∗ ,

δ(∇ϕ) = δ(∇ϕ) + Bδ rϕ∗ .

(47)

Substitution of Eqs. (46) and (47) into (44) and (45) gives the macro-scopic fluxes for moisture and heat M dj M = −DM ϕϕ d(∇ϕ) − Dϕ T d(∇ T ), M

dq

DM Tϕd

=−

(∇ϕ) −

DM TT d

(48)

(∇ T ),

(49)

M M M where DM ϕϕ , Dϕ T , DT ϕ , DTT are macro-scopic material conductivity matrices. The first equality can be then written for moisture m M m ∗ ∗ (−DM ϕϕ + Dϕϕ )d(∇ϕ) + (−Dϕ T + Dϕ T )d(∇ T ) + Lϕϕ drϕ + Lϕ T drT = 0,

(50)

and for heat m M m ∗ ∗ (−DM T ϕ + DT ϕ )d(∇ϕ) + (−DTT + DTT )d(∇ T ) + LT ϕ drϕ + LTT drT = 0,

(51)

with meso-scopic material conductivity matrices dw

 + δp psat ,  dϕ  = hv δp psat , 

Dm ϕϕ = Dw Dm Tϕ

 dp sat

Dm ϕ T = δp Dm TT

 ϕ(X 0 ) ,

  dT dp sat ϕ(X 0 ) , = λ + hv δp

(52)

dT where the superscript m denotes the meso-structural level. The second equality is then provided for moisture T Lϕϕ d(∇ϕ) + LϕT T d(∇ T ) + Dϕϕ drϕ∗ + Dϕ T drT∗ = 0,

(53)

and for heat T LTTϕ d(∇ϕ) + LTT d(∇ T ) + DT drϕ∗ + DTT drT∗ = 0.

(54)

In previous equations, matrices are defined by: dw

 Lϕϕ = LT ϕ =

Dw

  + δp psat B ,

dϕ

   dpsat δp ϕ(X 0 ) B , dT    dpsat 0 LTT = λ + hv δp ϕ(X ) B , dT    dpsat LϕT T = BT δp ϕ(X 0 ) , dT    dpsat 0 T T ϕ(X ) . LTT = B λ + hv δp Lϕ T =

hv δp psat B ,



 





T Lϕϕ = BT Dw

LTTϕ = B



dw dϕ

hv δp psat

 T

+ δp psat 



,

,

(55)

dT

Assuming that the heat and moisture balance is reached at the end of the ith step at the meso-structural level and the nodal unbalanced forces are equal to zero, the following system of Equations is obtained



drϕ∗



drT∗



D = − ϕϕ DT ϕ

Dϕ T DTT

−1 

T Lϕϕ

LϕT T

LTTϕ

T LTT



d(∇ϕ) , d(∇ T )



(56)

where Dϕϕ , Dϕ T , DT ϕ and DTT are material conductivity matrices at the meso-structural level (see Eqs. (19) and (20)). The system of equations (50) and (51) with the help of Eqs. (56) can be rewritten as



DM ϕϕ

DM ϕT

DM Tϕ

DM TT



d(∇ϕ) d(∇ T )



 =

Dm ϕϕ

Dm ϕT

Dm Tϕ

Dm TT





L − ϕϕ LT ϕ

Lϕ T LTT

Dm ϕϕ

Dm ϕT

Dm Tϕ

Dm TT



 −1 

T Lϕϕ

LϕT T

LTTϕ

T LTT

 

d(∇ϕ) d(∇ T )

 (57)

to be satisfied for any changes of the macro-scopic heat and moisture gradients. Then the macro-scopic material conductivity matrix becomes



DM ϕϕ

DM ϕT

DM Tϕ

DM TT



 =

Dm ϕϕ

Dm ϕT

Dm Tϕ

Dm TT





L − ϕϕ LT ϕ

Lϕ T LTT

Dm ϕϕ

Dm ϕT

Dm Tϕ

Dm TT



−1 

T Lϕϕ

LϕT T

LTTϕ

T LTT

 .

(58)

For the determination of the macro-scopic heat and moisture storage capacities which appear in the macro-scopic balance equations (40) and (41), the following consistency rules based on the volume averaging are adopted, see [34,11]:



dw dϕ

M =

1

|Ω |ϕ(X 0 )

  Ω

dw dϕ

m

ϕ dΩ ,

(59)

240

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

Fig. 3. Scheme of a processor farm used on a masonry arch bridge.



ρC +

∂ Hw ∂T

M =

1

|Ω |T (X 0 )

  Ω

ρC +

∂ Hw ∂T

m

T dΩ ,

(60)

where |Ω | is the area of the periodic unit cell, T (X 0 ) is the macro-scopic temperature, ϕ(X 0 ) is the macro-scopic moisture. 5. Processor farming method The processor farming method differs from the domain decomposition method presented in Section 3. The processor farm architecture assumes each macro-scopic integration point or each macro-scale finite element connected with a certain meso-scopic problem represented by an appropriate representative volume element (periodic unit cell). The solution of a meso-scale problem provides effective properties needed on the macro-scale. The macro-problem is assigned to the master processor while the homogenization at the meso-level is carried out on slave processors (Fig. 3). At each time step the current temperature, moisture and their gradients at a given macro-scopic integration point are passed to the slave processor (imposed onto the associated periodic cell), which, upon completing the meso-scale analysis, sends the homogenized data (effective conductivities, averaged storage terms and fluxes) back to the master processor. It means that the material conductivity and capacity matrices of an element or integration point of the macro-scale problem (Eqs. (40) and (41)) are assembled. The integration is performed numerically and the matrix product BT DB is evaluated at each integration point. Therefore, each matrix D is obtained as a result of the homogenization on meso-scale (Eqs. (58)–(60)). The meso-scale problems are solved n times, where n denotes the number of all integration points or elements. The algorithm of this approach is described in Table 3. If the meso-scale problems are large enough, the ideal solution is to assign one meso-problem to one slave processor. For very small macro-problems with a few thousands of finite elements, the hardware requirements would be in such a case excessive. On the other hand, if the meso-problems are relatively small, e.g., they contain a small number of finite elements, the corresponding analysis might be even shorter than the data transfer between the processors. Then, the communication time associated with the data transfer between the master processor and many slave processors may grow excessively. It is worth mentioning that this time consists of two contributions. The first one represents the latency time (the processors make connection) which is independent of the amount of transferred data while the second contribution clearly depends on the amount of data being transferred. For small meso-problems, it is therefore reasonable to assign one type of meso-problems or several of them to a single slave processor. The master processor then sends a larger package of data from many macroscopic integration points at the same time to each slave processor so that the latency time does not play a crucial role. The ideal speed-up and load balancing are obtained when the decomposition of the macro-problem reflects the meso-problem meshes. However, this is considerably more difficult when compared to the classical mesh decomposition. From the physical point of view, it is luxury to solve a meso-scale problem for each integration point of the macro-scale problem. In case of large meso-scale problems, it is appropriate to use the concept of aggregated elements, where several neighbouring elements of the same material types are aggregated. The meso-scale problem can be solved for the whole aggregate thus considerably saving time. On the other hand, the aggregation has to take into account material interfaces and other sources of possible discontinuities. This approach is suitable if the temperature and relative humidity in all elements of the aggregate are comparable. The application of aggregates is shown in a hygro-thermal analysis of Charles bridge in Prague presented in Section 6.2. 6. Results of computation 6.1. Heat and moisture transfer analysis of a brick wall The benchmark study for the methods presented is a 2D hygro-thermal analysis of a brick wall consisting of ceramic blocks with the cavities filled by thermal insulation material and mortar layers on both sides (Fig. 4). This numerical example

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

241

Table 3 Algorithm of parallel processor farming method for two-scale coupled heat and moisture analysis.

Fig. 4. Scheme of brick wall.

follows the former studies published in [35,36] and compares the efficiency of the proposed processor farming method against the classical domain decomposition method. For the first, the 2D numerical model was discretized by the finite element mesh with triangular elements with linear approximation functions. The mesh is respecting the cavities filled with thermal insulation material and the inner shaped walls of the brick block and has 216,806 elements and 109,082 nodes. The domain decomposition was tested for different numbers of subdomains and computational times of one time step for various mesh decompositions are summarized in Table 4, where the ideal time is the total computational time without time of data files reading and printing, respectively. The decomposition of the brick wall is illustrated by Fig. 5. The domain is decomposed into subdomains which should contain the same number of nodes and elements in ideal case. The mesh was decomposed by METIS and the number of nodes and elements are nearly identical in all subdomains. The reason of different size and shape of subdomains is caused by different size of finite elements in the original mesh which is enforced by complicated shape of the brick. Moreover, the decomposition has to take into account material interfaces in the original domain. There are plasters and the brick itself and a single subdomain cannot contain two different materials. The analysis of the brick wall was performed on the DELL computer Precision T5600 equipped with two INTEL XEON processors E5 (2 × 8 cores) 2.4 GHz and RAM 16 GB (4 × 4 GB) 1600 MHz DDR3. The operating system is the Debian Linux 64-bit architecture. The computational time of determination of material properties and one time step took 3 min and 30 s and the whole computational time took about 4 days. Subsequently for the processor farming method, the 2D model was discretized by a reasonable coarse mesh of the macroproblem with 948 quadrilateral elements and 9227 nodes (Fig. 6 left) and the meso-problem, which is assigned to each finite element of the macro-problem (only for the brick block and the cavities), has 2597 elements and 2872 nodes. Moreover, the

242

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248 Table 4 Computational time for various mesh decompositions. Number of subdomains N Time (s) Ideal time (s)

1 211 210

4 52 43

10 29 21

20 17 11

Table 5 Comparison of computational times for three different methods. Solution type

One time step

Total time

Classical FE solution (216,806 elements and 109,082 nodes) Domain decomposition method with 20 subdomains Processor farm method with 20 subdomains

3 min and 30 s 17 s 28 s

4 days 10 h 15 h

Fig. 5. Domain decomposition of finite element mesh of brick wall—20 subdomains.

mesh of macro-problem was decomposed into 20 subdomains (aggregates), where the appropriate meso-problems were assigned to 20 slave processors. If the meso-scale problem is too large, the basic assigning of one meso-problem to each integration point of the macro-level problem is not the ideal solution. For this case, it is better to decompose the macroscale problem in a similar way as a classical domain decomposition method, where one subdomain (aggregate) is assigned to one meso-scale problem (and therefore also to one slave processor) respecting material types and the water vapour pressure (moisture) and temperature fields distribution. One meso-scale problem is then computed only once in one time step for the whole aggregate with averaged water vapour pressure, temperature and their gradients as input parameters entering the meso-problem. Subsequently, effective material properties are redistributed to each element of the aggregate. This approach can reduce significantly the computational time and can be applied if the temperature and relative humidity in all points of the aggregate are comparable. On the other hand, the aggregation has to take into account material interfaces and other sources of possible discontinuities in the response. The computational time of determination of material properties at meso-structural levels and one time step at macrolevel for the problem took 28 s and the whole computational time for 1800 time steps took about 15 h. The reduction of the computational time is evident in comparison with one subdomain in the decomposition method (one time step took 3 min. and 30 s). The speed up of processor farming method with two-level homogenization is not too significant in comparison with domain decomposition method but it has to be mentioned that the finite element mesh for the meso problem is much more finer than original mesh. For clarity, the computational times are compared in Table 5. Results of the computation are illustrated by the temperature and water vapour pressure history on the boundary between the brick block and the mortar near the internal surface (Figs. 7 and 8). The temperature and water vapour profiles in different times are depicted in Figs. 9 and 10. The good coincidence between numerical simulations using the classical FEM approach with domain decomposition method and the homogenization with parallel computing can be observed. The small differences in temperature and water vapour pressure history and profiles are influenced by the decomposition of the macro level into aggregates and by the size of finite element mesh at the macro-scale level. 6.2. Heat and moisture transfer analysis of Charles bridge in Prague In this section, the practical application of three-dimensional two-scale homogenization with processor farming method is presented. The 3D analysis is a part of a complex hygro-thermo-mechanical analysis of Charles Bridge in Prague (Fig. 11) which has been developed to determine the influence of climatic loading on the current state of the bridge [1]. The analysis established that climatic loading should be considered the most serious for the bridge because it is responsible for the nucleation and further development of cracks in the bridge. This fact is supported by results of former complex analysis of stone bridges and other masonry structures, see, e.g., [37–39] to mention a few references. Furthermore, it was shown that such a suitable thermo-mechanical (or hygro-thermo-mechanical) analysis thus becomes the first thing when

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

243

Fig. 6. Macro-level finite element mesh with representative volume element at meso-level of brick wall.

Fig. 7. History of temperature on the boundary between the brick block and the mortar on the right-hand side of the wall.

Fig. 8. History of pore water pressure on the boundary between the brick block and the mortar on the right-hand side of the wall.

contemplating any actions such as, e.g., designing further stages of the bridge repair. In case of Charles Bridge, it is the retrofit of the bridge’s sandstone cladding scheduled for the near future.

244

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

Fig. 9. Temperature distribution in the wall.

Fig. 10. Pore water pressure distribution in the wall.

The three-dimensional numerical model simulates the coupled heat and moisture transfer of one half of arch No. 3. The finite element mesh of the macro-problem was created using tetrahedron elements with linear approximation functions. It has 73,749 nodes and 387,773 elements (Fig. 12). Meso-problem is assigned to each finite element of the macro-problem. The mesh of macro-problem was decomposed into 15 subdomains according to material types and the appropriate mesoproblems were aggregated and assigned to 15 slave processors. The numbers of elements in subdomains which are equal to the numbers of meso-problems in groups are summarized in Table 6. Regular masonry (arch vault, pier and breast walls) is represented by SEPUC No. 1 (Fig. 1) which is formed from 1950 nodes and 1512 hexahedron elements with linear approximation functions, while irregular masonry (infill) is represented by SEPUC No. 2 (Fig. 2) with 2454 nodes and 12 269 tetrahedron elements with linear approximation functions. Other parts of the bridge (concrete slab and layers of pavement) are considered as homogeneous and isotropic. The decomposition of the macro-problem into subdomains is not ideal. In comparison with domain decomposition methods, the macro-problem has to be split with respect to heterogeneity of the material and therefore the number of elements in subdomains varies from 11,563 to 37,677. The ideal speedup and load balancing are obtained when the decomposition of the macro-problem into subdomains takes into account the meshes used in meso-problems. It should be noted that it is much more difficult task than the classical mesh decomposition. The analysis of Charles bridge was performed on two powerful computers each equipped with Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40 GHz with 32 cores and RAM 128 GB. The computers are based on Debian Linux and 64-bit architecture. Despite of such powerful configuration of these two computers, the computational time of determination of material properties at meso-structural levels and one time step at the macro-level for coupled heat and moisture transfer analysis took 10 h. It is evident, if the meso-scale problem is too large, the assigning of thousands of elements of one meso-problem to one slave processor is not the ideal solution. As already mentioned, for relatively small macro-problems with a few

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

245

Fig. 11. Detail of Charles bridge in Prague, Czech Republic.

Fig. 12. 3D finite element mesh model of one half of the arch No. 3 (left) and material decomposition of macro-problem (right).

thousands of finite elements, the hardware requirements would be in such a case excessive. For this case, it is better to decompose the macro-scale problem in a similar way as a classical domain decomposition method, where one subdomain (aggregate) is assigned to one slave processor. One meso-scale problem is then computed only once in one time step for the whole aggregate with averaged moisture, temperature and their gradients as input parameters entering the mesoproblem. Subsequently, effective material properties are redistributed to each element of the aggregate. Fig. 13 shows the decomposition of the finite element model to 63 aggregates respecting material types and the moisture and temperature fields distribution. This decomposition reduced significantly the length of one time step to 7 min. It should be mentioned, the first step of the hygro-thermal analysis was the verification and validation of heat and moisture transfer model and its material parameters. The boundary conditions for this analysis were specified with respect to climatic data, including the effect of sun radiation, wind, rain, heat conditions, and the structure’s orientation. Subsequently, the 3D non-stationary evolution of temperature and moisture fields in Charles Bridge was analysed in conjunction with its last repair, when a measuring system was installed in the bridge and the temperature and moisture content at selected gauge points have been continuously monitored for a period of two years. Results of computation are illustrated by the temperature distribution in the bridge arch No. 3 in the spring season depicted in Fig. 14. The comparison of temperature computed with in situ measurements in selected point is shown in Fig. 15, where the temperature profile is in good agreement with in situ measurements. More details about the hygrothermo-mechanical analysis of Charles bridge can be found in Refs. [1,40,41].

246

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

Table 6 Decomposition of the Charles bridge macro-problem into 15 subdomains. Aggregate No.

Bridge element

# elements

PUC

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Concrete slab Arch vault Pier Infill Breast wall Breast wall Infill Paving 2nd concrete layer 1st concrete layer Infill Pier Pier Pier Pier

11 198 31 962 37 677 19 351 24 267 23 745 26 812 11 563 11 589 12 088 26 813 37 677 37 677 37 677 37 677

– SEPUC1 SEPUC1 SEPUC2 SEPUC1 SEPUC1 SEPUC2 – – – SEPUC2 SEPUC1 SEPUC1 SEPUC1 SEPUC1

Fig. 13. Decomposition of arch No. 3 macro-problem into 63 subdomains (aggregates).

Fig. 14. Temperature profile [Kelvins] in bridge body.

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

247

Fig. 15. Distribution of temperature at selected gauge point [degree centigrade].

7. Conclusion The well-tried constitutive model by Künzel and Kiessl was addressed in conjunction with two parallel computing methods which were introduced for numerical analyses of coupled heat and moisture transfer in masonry structures. Both the first, the domain decomposition method, and the second processor farming method with two-scale approach seem to be effective tools for large problems with huge number of degrees of freedom. Both methods presented also exploit the master/slave processor architecture and MPI communication format. Whereas the domain decomposition method splits the original domain into non-overlapping subdomains and uses the Schur complement method for the solution of non-linear and non-symmetric system of equations, the processor farming method is based on the first order homogenization approach and assigns macro-scale problems to the master processor and meso-scale problems to slave processors. The speed-up of both methods used is evident in comparison with the original finite element discretization. The results of the 2D numerical simulation of coupled heat and moisture transfer in a brick wall show good agreement in both parallel methods. The small visible differences are particularly influenced by the quality of finite element mesh used at the macro-level and by the size of aggregates of the processor farming method, were the averaging of non-linear distribution of moisture and temperatures and their gradients play the crucial role. The main advantage of both parallel methods is the significant reduction of the computational time. Comparing both methods, the processor farming method in connection with homogenization procedure seems to be more effective than domain decomposition, because it allows to bring randomness of material composition into the computation. Despite these facts, domain decomposition method and processor farming method enable to solve more complicated geometries of real world masonry and other heterogeneous structures. Acknowledgements Financial support for this work was provided by GAČR project No. 15-17615S. This financial support is gratefully acknowledged. References [1] T. Krejčí, J. Šejnoha, Evolution of temperature and moisture fields in Charles Bridge in Prague, Computational prediction and measurements, Int. J. Archit. Herit.: Cons. Anal. Restor. 9 (8) (2015) 973–985. [2] C. Miehe, C.G. Bayreuther, On multiscale FE analyses of heterogeneous structures: From homogenization to multigrid solvers, Internat. J. Numer. Methods Engrg. 71 (2007) 1135–1180. [3] A. Toselli, W. Olof, Domain Decomposition Methods–Algorithms and Theory, in: Springer Series in Computational Mathematics, vol. 34, Springer, Berlin, 2005. [4] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, in: Numerical Mathematics and Scientific Computation, Oxfor Science Publications, Clarendon Press, Oxford, 2005. [5] J. Kruis, Domain Decomposition Methods for Distributed Computing, Saxe-Coburg Publications, Kippen, Stirling, Scotland, UK, 2006. [6] C. Farhat, F.X. Roux, A method of finite element tearing and interconnecting and its parallel solution algorithm, Internat. J. Numer. Methods Engrg. 32 (1991) 1205–1227. [7] C. Farhat, J. Mandel, The two-level FETI method for static and dynamic plate problems - Part I: An optimal iterative solver for biharmonic systems, Comput. Methods Appl. Mech. Engrg. 155 (1998) 129–152. [8] C. Farhat, P.S. Chen, J. Mandel, F.X. Roux, The two-level FETI method - Part II: extension to shell problems, parallel implementation and performance results, Comput. Methods Appl. Mech. Engrg. 155 (1998) 153–180. [9] C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, D. Rixen, FETI-DP: A dual-primal unified FETI method–Part I: A faster alternative to the two-level FETI method, Internat. J. Numer. Methods Engrg. 50 (7) (2001) 1523–1544. [10] Z. Dostál, D. Horák, R. Kučera, Total FETI–an easier implementable variant of the FETI method for numerical solution of elliptic PDE, Int. J. Numer. Methods Biomed. Eng. 22 (12) (2006) 1155–1162. [11] I. Özdemir, W. Brekelmans, M. Geers, Computational homogenization for heat conduction in heterogeneous solids, Internat. J. Numer. Methods Engrg. 73 (2008) 185–204.

248

T. Krejčí et al. / Computers and Mathematics with Applications 74 (2017) 229–248

[12] J. Sýkora, J. Vorel, T. Krejčí, M. Šejnoha, J. Šejnoha, Analysis of coupled heat and moisture transfer in masonry structures, Mater. Struct. 42 (8) (2009) 1153–1167. [13] R. Černý, P. Rovnaníková, Transport Processes in Concrete, Spon, London, 2002. [14] M. Hassanizadech, W. Gray, General conservation equations for multiphase systems: 1, averaging procedure, Adv. Water Resour. 2 (1979) 131–144. [15] M. Hassanizadech, W. Gray, General conservation equations for multiphase systems: 2, mass, momenta, energy and entropy equations, Adv. Water Resour. 2 (1979) 191–203. [16] M. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155–164. [17] R. De Boer, Highlights in the historical development of the porous media theory: toward a consistent macro-scopic theory, Appl. Mech. Rev. 49 (1996) 201–262. [18] R.W. Lewis, B.A. Schrefler, The Finite Element Method in Static and Dynamic Deformation and Consolidation of Porous Media, John Wiley & Sons, Chiester-Toronto, 1998. [19] H.M. Künzel, K. Kiessl, Calculation of heat and moisture transfer in exposed building components, Int. J. Heat Mass Transfer 40 (1997) 159–167. [20] M. Kaminski, Homogenization of transient heat transfer problems for some composite materials, Internat. J. Engrg. Sci. 41 (2003) 1–29. [21] B. Tomková, M. Šejnoha, J. Novák, J. Zeman, Evaluation of effective thermal conductivities of porous textile composites, Int. J. Multiscale Comput. Eng. 6 (2) (2008) 153–168. [22] J. Sýkora, M. Šejnoha, J. Šejnoha, Homogenization of coupled heat and moisture transport in masonry structures including interfaces, J. Appl. Math. Comput. 219 (13) (2013) 7275–7285. [23] J. Sýkora, J. Vorel, M. Šejnoha, J. Šejnoha, Effective material parameters for transport processes in historical masonry structures, in: Proceedings of the Tenth International Conference on Civil, Structural and Environmental Engineering Computing. Stirling, UK, Paper 190, 2005. [24] J. Kruis, T. Koudelka, T. Krejčí, Multi-physics analyses of selected civil engineering concrete structures, Commun. Comput. Phys. 12 (2012) 885–918. [25] V. Kočí, J. Maděra, R. Černý, Computer aided design of interior thermal insulation system suitable for autoclaved aerated concrete structures, Appl. Therm. Eng. 58 (1) (2013) 165–172. [26] V. Kočí, J. Maděra, R. Černý, Exterior thermal insulation systems for AAC building envelopes: computational analysis aimed at increasing service life, Energy Build. 47 (1) (2012) 84–90. [27] T.J.R. Hughes, The Finite Element Method. Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1987. [28] Z. Bittnar, J. Šejnoha, Numerical Methods in Structural Mechanics, ASCE Press, New York, USA, 1996. [29] M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, John Wiley & Sons Ltd., Chichester, UK, 1991. [30] A. Anthoine, Derivation of the in-plane elastic characteristics of masonry through homogenization theory, Internat. J. Solids Structures 34 (11) (1995) 137–163. [31] M. Lombardo, J. Zeman, M. Sejnoha, G. Falsone, Stochastic modeling of chaotic masonry via meso-structural characterization, Int. J. Multiscale Comput. Eng. 7 (2) (2009) 171–185. [32] J. Zeman, M. Šejnoha, From random microstructures to representative volume elements, Modelling Simul. Mater. Sci. Eng. 15 (4) (2007) 325–335. [33] F. Larsson, K. Runesson, F. Su, Variationally consistent computational homogenization of transient heat flow, Internat. J. Numer. Methods Engrg. 81 (2010) 1659–1686. [34] J. Sýkora, Multiscale modeling of transport processes in masonry structures (Ph.D. thesis), Czech Technical University in Prague, 2010. [35] V. Kočí, J. Maděra, J. Jerman, M. Trník, R. Černý, Determination of the equivalent thermal conductivity of complex material systems with large-scale heterogeneities, Int. J. Therm. Sci. 86 (2014) 365–373. [36] J. Maděra, J. Kočí, V. Kočí, J. Kruis, Parallel modeling of hygrothermal performance of external wall made of highly perforated bricks, Advances in Engineering Software, in print, http://dx.doi.org/10.1016/j.advengsoft.2016.08.010. [37] P.B. Lourenço, Computations of historical masonry constructions, Prog. Struct. Eng. Mater. 4 (3) (2009) 301–319. [38] P.B. Lourenço, J. Krakowiak, F.M. Fernandes, L.F. Ramos, Failure analysis of monastery of jerónimos, Lisbon: How to learn from sophisticated numerical models? Eng. Fail. Anal. 14 (2) (2007) 280–300. [39] A.H. Akhaveissy, G. Milani, A numerical model for the analysis of masonry walls in-plane loaded and strengthened with steel bars, Int. J. Mech. Sci. 72 (2013) 13–27. [40] M. Toesca, Thermo-Mechanical Analysis of Charles Bridge. Faculty of Civil Engineering (Master thesis), Czech Technical University in Prague, 2014. [41] T. Krejčí, J. Šejnoha, M. Toesca, in: J. Kruis, Y. Tsompanakis, B.H.V. Topping (Eds.), Coupled Hygro-Thermo-Mechanical Analysis of Charles Bridge, Prague, Computational Techniques for Civil and Structural Engineering, in: Computationas Science, Engineering and Technology Series, vol. 38, Saxe-Coburg Publication, Stirlingshire, Scotland, 2015.