Abstract

Various simplification or optimization techniques are sought that reduce demands of computational modeling on time or computing power while keeping a sufficient level of accuracy. In this paper, determination of hygrothermal performance of a brick block is presented using two homogenization techniques based on different principles. While the computational homogenization technique uses a multiscale method realized on the master/slave computer system, the materials homogenization comes out from the effective media theory (EMT), and after the determination of effective material properties, the whole isotropic problem can be transformed to one dimension. Contrary to most applications of EMT, free parameters of mixing formulas are not determined based on an experimental measurement of a single material property but on a complex hygrothermal performance of the element where the distribution of moisture and temperature over a reference year is taken into account. The calculated results from both techniques are compared with results obtained by high-performance computing without any computational simplifications. For materials homogenization, the best results are achieved when k = 0.9 in Lichtenecker’s mixing rule is assumed, which corresponds to a nearly parallel arrangement of the block. The root mean square error (RMSE) of relative humidity (RH) and temperature distribution is only 0.992% and 0.566°C, respectively. This is even better than the results of computational homogenization (RMSE: 1.502% of RH and 0.629°C). Besides obtaining sufficiently precise results, a significant time-saving is achieved, accounting for more than 99%, while being solved on a single-processor computer.

1. Introduction

Computational modeling of hygrothermal performance represents nowadays a very powerful tool for the assessment of building materials, elements, or structures. As Medjelekh et al. [1] or Brischke and Thelandersson [2] reported in their reviews, the computational approach has become an indisputable part of building materials assessment during the past several decades. Since this approach is supposed to be time-saving and cost-effective at the same time, it is predisposed to be used even more in the nearest future. The possibility of predicting hygrothermal behavior of building materials enables a broad range of practical applications as heat or moisture transport is involved in many processes resulting in deformation or damage of building materials [36].

It cannot be stated that computational modeling provides absolutely precise results, in general. The accuracy of input parameters, material properties in particular, and the quality of applied physical and mathematical models can be considered as the most important conditions for achieving that goal. The application of advanced experimental techniques for acquisition of material characteristics in combination with the utilization of well-calibrated models being able to reach an excellent match between computational and experimental results presents the most effective tools in that respect. There are numerous examples of a successful use of such approaches. For instance, Hassani et al. [7] developed a rheological model for wood that is able to investigate multispecies structures when exposed to moisture or mechanical straining. Using this model, they were able to predict deformations of the glued-laminated wood sample after several drying steps very precisely. The numerical results they achieved were confirmed by experimental measurements. Using computational modeling, moisture-induced stresses and drying shrinkage in concrete were predicted by Azenha et al. [8], while the obtained results were satisfactorily validated by means of experimental measurements. Other computational models that have been successfully validated were reported, e.g., by Ullah et al. [9], Kacmarczyk et al. [10], or Bozorgzad et al. [11].

All the mentioned researches, among many others, utilized the finite element or finite volume method for numerical solution of used models to achieve results as precise as possible. On the contrary, such a high level of accuracy is often balanced by high time demands or by extensive requirements on computing power [9, 12] because very detailed finite element discretizations usually result in very high numbers of degrees of freedom. Handling such extensive problems may often be beyond the possibilities of many research departments. Therefore, alternative computational approaches are being sought that would bring a compromise between computing time, requirements on computing power, and accuracy of obtained results. Within this paper, several methods are described. Being demonstrated on modeling of hygrothermal performance of a selected example, the best method is proposed on the basis of the previously mentioned criteria.

2. Materials

A highly perforated ceramic block has been chosen as a representative of heterogeneous systems of which hygrothermal performance must be primarily solved at least as two-dimensional. These bricks, which may contain different cavity fillings, have become very popular in the building market because of their excellent thermal-insulating properties as well as very good mechanical parameters.

Within this research, a ceramic block filled with expanded polystyrene has been selected which was provided with a lightweight plaster containing perlite from both sides. Because the main heat and moisture transport processes in the block are supposed to be in the horizontal plane, the simulation was performed as two-dimensional with the emphasis put on preservation of the sophisticated integral structure of the block. The material properties were determined experimentally following the laboratory procedures proposed in [13]. The summary of data, which was taken from the study in [1416], is given in Table 1, where ρ is the bulk density (kg·m−3), ψ is the open porosity (%), λ is the thermal conductivity (W·m−1·K−1), c is the specific heat capacity (J·kg−1 K−1), μ is the water vapor diffusion resistance factor (—), κapp is the apparent moisture diffusivity (m2·s−1), and is the hygroscopic moisture content (m3·m−3). Having the cavity/body ratio equal to 0.558, the investigated block is captured in Figure 1.

3. Methods

The hygrothermal performance of the ceramic block can be predicted using the modified Künzel’s model [17] of which heat and moisture balance equations are formulated aswhere ρw (kg·m−3) is the water density, (m3·m−3) is the water content, pv (Pa) is the partial pressure of water vapor in the air, t (s) is the time, Dg is the global moisture transport function, H (J·m−3) is the enthalpy density, T (K) is the temperature, λ (W·m−1·K−1) is the thermal conductivity, Lv (J·kg−1) is the latent heat of evaporation of water, and δp (s) is the water vapor permeability.

The simulation of coupled heat and moisture transport was done for the time period of three years. The climatic boundary conditions were applied in the form of TRY (test reference year) which was provided by the Czech Hydrometeorological Institute according to ČSN EN ISO 15927-4 [18] on the basis of long-term average values of various weather parameters, such as temperature, relative humidity, precipitation, wind, and several kinds of sun radiation.

Two-dimensional modeling of the brick resulted in 85,158 nodes and 84,394 quadrilateral elements. In the case of more complicated geometry, the number of nodes and unknowns grows significantly and large systems of algebraic equations are obtained. Such systems cannot be solved on single-processor computers because of insufficient computer memory and extremely long computational time. Therefore, parallel computers have to be used. Detailed analysis of numerical algorithms usually reveals that the main bottleneck for parallelization of such tasks is the system of algebraic equations. One possible group of methods enabling parallelization of large systems is the group of domain decomposition methods. The original domain solved is decomposed into smaller subdomains, and the continuity on the subdomain interfaces has to be enforced. For the purposes of this paper, the Schur complement method was selected because it works for symmetric as well as nonsymmetric systems. Some material models of transport processes lead to nonsymmetric systems of equations. This method is further described in Section 3.1. Other alternatives are represented by applications of simplifying methods that reduce the computing time on the expense of accuracy. In these cases, it is crucial to set a compromise between required computing time and power on the one side and accuracy of obtained results on the other side. Two different techniques based on computational homogenization (Section 3.2) and materials homogenization (Section 3.3) are studied within this paper, and the results are compared and discussed.

3.1. Domain Decomposition Method

Parallel computing using the domain decomposition method presents one of the ways on how to deal with the excessive demands arising at the solution of complex problems with many unknowns. For example, the application of the Schur complement method for solving large nonsymmetric systems obtained from spatial discretization of partial differential equations describing hygrothermal phenomena in building materials was proposed by Kruis [19]. For the solution of heterogeneous systems, the multigrid method suggested by Miehe and Bayreuther [20] can be exploited in connection with parallel computing because it often cannot be executed on single-processor computers.

The Schur complement method is based on decomposition of a domain into subdomains, where the internal unknowns are eliminated from the system of equations and finally the continuity on the subdomain interfaces is enforced. Let a system of linear algebraic equations have the following special form:where is the vector of all internal unknowns defined in the jth subdomain, is the vector of all unknowns defined on subdomain interfaces, is the vector of the right-hand side of the jth subdomain, is the vector of the right-hand side connected with all nodes on subdomain interfaces, and and are the matrices connected with the jth subdomain. The special form of system of equation (3) can be obtained by a special ordering of unknowns. The unknowns in inner nodes are ordered first, and the unknowns in interface nodes are ordered last. This type of ordering assures nonsingularity of the diagonal blocks, and therefore, all internal unknowns defined on any subdomain can be eliminated. For example, the internal unknowns on the jth subdomain can be expressed in the following form:

Clearly, the internal unknowns for all subdomains can be eliminated simultaneously because there is no influence of the internal unknowns of other subdomains. Only the vector of all unknowns on interfaces is there. Substitution of all expressions from equation (4) into the last equation (3) leads to the reduced system:where only the block is unknown. When the reduced system (5) is solved, the vector is substituted into the expression in equation (4) and all unknowns can be solved.

The reduced system (5) can be solved by a finite method or by an iterative method. Typical example of the finite method is the LDLT or LU factorization, and the iterative methods can be represented by the conjugate gradient method or the biconjugate gradient method.

3.2. Computational Homogenization Technique Based on Multiscale Analysis and Processor Farming Method

One group of alternative approaches is based on the utilization of computational homogenization techniques. A multiscale approach based on computational homogenization belongs to this group. The computational methods for multiscale modeling of heterogeneous structures have emerged in recent decades. The original approach was proposed by Kaminski [21] who selected several composite materials for homogenization of transient heat transfer problems. The complex multiscale analysis was described by Ozdemir et al. [22] who aimed also at heat transfer parameters of solids. Analyzing masonry structures, Sykora et al. [23] adopted the previously mentioned techniques and used the homogenization procedure at the meso-level (a composition of stone blocks and mortar in masonry) which provided effective (macroscopic) transport parameters necessary for a detailed analysis of the whole structure at the macro-level. Contrary to Kaminski [21] or Ozdemir et al. [22], they involved moisture transport as well. Very similar principles were presented also by Larsson et al. [24] and Lu et al. [25] who simulated cracks in solids using the multiscale finite element method at the macroscale and mesoscale.

The multiscale method using multiscale analysis was one of the homogenization methods used for the hygrothermal simulation of the brick wall. The two-scale analysis comes out from the first-order homogenization method, where the function values and their gradients at the macro-level are used at the lower meso-level. The finite element method is used for the solution of governing equations of the coupled heat and moisture transport (equations (1) and (2)).

Overall-homogenized macroscale properties are found from the solution of the mesoscale steady-state problem carried out on the representative volume element (RVE). RVE is usually represented by a statistically equivalent periodic unit cell (SEPUC) which is constructed to match the real mesostructure in a statistical sense as close as possible. The effective capacity matrix given byand the conductivity matrix given byneeded at the macro-level are solved by the multiscale method which exploits processor farming based on the master/slave architecture of the computational cluster. In equations (6) and (7), are tensors of homogenized conductivity, CT and Cφ are capacity submatrices, contain homogenized capacity coefficients, the matrix N contains approximation functions for the temperature and the relative humidity, and BT and Bφ are matrices of gradients of approximation functions stored in the matrix N.

In this method, each macroscopic integration point or each macroscale finite element is connected with a certain mesoscopic problem represented by the RVE. The macro-problem is assigned to the master processor, while the homogenization at the meso-level is carried out on slave processors. In each time step, the actual temperature, moisture, and their gradients at a given macroscopic integration point are passed to the slave processor (imposed onto the associated RVE), which, upon completing the mesoscale analysis, sends the homogenized data (effective conductivity and capacity matrices ) back to the master processor.

Similarly, as in Domain Decomposition Method (Section 3.1), the speedup and load balancing are important indicators of effective computations. The ideal situation is obtained when the decomposition of the macro-problem reflects the meso-problem meshes. For instance, the ideal solution is to assign one meso-problem to one slave processor. If the meso-problems are relatively small, i.e., they contain a small number of finite elements, the corresponding analysis might be even shorter than the data transfer between the processors. In this case, the communication time associated with the data transfer between the master and slave processors may grow extremely. It is therefore reasonable to assign one type of meso-problem or several of them to single slave processor. The master processor then sends a larger package of data from many macroscopic integration points at the same time to slave processors. On the contrary, in the case of large mesoscale problems, the concept of aggregated elements on the macro-problem is recommended for the purpose to save computational time. In this approach, several neighboring elements of the same material types are aggregated, and the mesoscale problem can be solved only once for the whole aggregate. The aggregation has to fulfill several recommendations, as for the material interfaces and discontinuities, high gradients of unknowns, and others, which are discussed, e.g., in [26].

Within this research, the assumed block was decomposed into 29 subdomains—aggregates (Figure 2), being solved using 29 slaves and one master processor. It can be seen in Figure 2 that the macro-problem mesh is much coarser, being formed only by 4,330 nodes and 4,049 quadrilateral elements with linear approximation functions.

The different sizes and shapes of subdomains are caused by different sizes of finite elements in the mesh which is enforced by the complex shape of the brick. The meso-problem (Figure 3), which is assigned to aggregates of finite elements of the macro-problem, has 1,071 nodes and 1,004 quadrilateral elements only for the brick block and the cavities.

3.3. Materials Homogenization Technique Based on Effective Media Theory

The main contribution of materials homogenization lies in a fact that mixing rules based on effective media theory convert the multicomponent areas into single homogenized matter. Therefore, it allows for a transformation of a given isotropic problem to a less-dimensional one which is computationally much less demanding, regardless of the complexity of the solved detail. Even if this approach definitely affects the accuracy of obtained results, it significantly reduces the computing time, being less demanding on computing power.

There are lots of records in the scientific databases dealing with the application of mixing rules in the field of materials science. For example, Karinski et al. [27] used in their research both parallel (Voigt) and series (Reuss) models for prediction of mechanical properties of cement-based composites. However, before the application of the models, they had to assume a homogeneous stress field in the specimen. These models corresponded to the boundaries of the Wiener system [28]. The Wiener theory, in general, is used for the characterization of a system by a set of functions which are determined by analyzing the system response. It has been applied in several material engineering tasks as reported, e.g., by Suleiman [29] or Hossain et al. [30]. This theory was also adopted in Lichtenecker’s mixing rule [31, 32] that operates over the whole range between parallel and serial models by adjusting the k parameter between −1 and 1. Lichtenecker’s formulas applied to a variable F of a mixture of N phases, with Fi and fi as the variables related to the ith phase and volume fraction, respectively, result in the “k model” given byand the “mean geometrical” model given bywhere F is the effective homogenized parameter, Fi is the parameter of the ith element of the system, and fi is the volume fraction of the ith element of the system. When considering porous building materials, one is then able to describe different spatial arrangements. Lichtenecker’s mixing rule is considered to be rather empirical. It was successfully applied to express physical properties of multiphase mixtures [33]. Probably, the weakest point is a proper determination of its parameters that ensure the correct results. The parameter k may be considered as describing a transition from the anisotropy at k = −1 to another anisotropy at k = 1. However, Lichtenecker’s equation may also be applied to isotropic composites. In this case, the parameter k must be determined based on the literature research, expert estimations, or experimental testing. In this paper, the global parameter k is treated. It means a single k parameter is determined for all the transport and storage parameters of the block based on the long-term hygrothermal performance, trying to reach the highest match with the results obtained using the domain decomposition method.

In order to select the optimal value of the k parameter from the point of view of the global performance of the block, all the variations of k = {−1.0; −0.9; −0.8; … 0.8; 0.9; 1.0} were investigated and the complete sets of effective material parameters were determined and used in the simulations. Such a global k value was then selected that provided an identical response to the domain decomposition method presented in Section 3.1. The effective bulk density and open porosity were calculated on the basis of the known geometry of the block.

4. Results and Discussion

The main objective of the computations was to determine the hygrothermal performance (temperature and moisture content with respect to time) of the studied ceramic block with the complex geometry. To allow the comparison, the results were related to a selected profile (Figure 4) because one of the applied methods is based on the materials homogenization and thus converts the isotropic problem to one dimension.

The domain decomposition method was performed on a PC cluster of fourteen Dell Precision T5600 computers equipped with two Intel Xeon E5 processors (2 × 8 cores) with 3.3 GHz and 1600 MHz DDR3 16 GB RAM (4 × 4 GB). The master/slave system of computers connected within the cluster enabled a decomposition of the computational mesh (Figure 4) into a number of smaller subdomains that were solved in parallel. According to Madera et al. [12], whose principle was followed within this technique, the decomposition leads to superlinear speedup which is caused by smaller bandwidth of subdomain matrices in comparison with the bandwidth of the original matrix. According to the performed computational time test [12], the highest relative speedup is achieved when 4 subdomain matrices are used. However, the further increase of subdomains means additional time reduction, even if it is not as efficient as in the case of smaller numbers. Within this paper, the mesh was decomposed into 30 subdomains. This setup reduced computational time by about 86% when compared to classical single-processor solution, but the computations of three-year performance took 982,020 s (11 d, 8 h, and 47 m) anyway. The accuracy of obtained results depends, of course, on the accuracy of input data, as well as on the accuracy of the used mathematical model. Even if such an assessment was not the objective of this paper, one can refer to the model verification performed in [17] which confirmed a very good match of results of a laboratory critical experiment and its computational representation (Figure 5).

The accuracy of the model is very high, having R2 between 0.9870 and 0.9888 and RMSE not higher than 0.53°C and 1.36% of RH [17]. Therefore, the results obtained using the domain decomposition method can be considered to be very precise, being balanced by the extreme requirements on computing time. Additionally, the identical input data together with the same model and boundary conditions were used for all techniques applied in this paper so that one could focus only on the application of the simplifying techniques, or more specifically on their influence on the results from the point of view of the block performance.

The computational homogenization based on the multiscale approach using the processor farming method was carried out on the same cluster as the domain decomposition method. The three-year performance results were obtained in 200,398 s (2 d, 7 h, and 40 m). It represents approximately 79.6% of time-saving when compared to the domain decomposition method.

The global k parameter in Lichtenecker’s mixing rule was selected based on a robust analysis which comprised a comparison of temperature and relative humidity values in each node of the selected profile (Figure 4) with the values obtained using the domain decomposition method after every time step in the third year of simulations. Since the time step of calculations was one hour and there were 94 nodes in the reference profile, it meant that 823,440 results were compared and evaluated. This approach ensured that every deviation from the reference performance was reflected and taken into account in the final evaluation. The lowest RMSE (root mean square error) then indicated which value of the k parameter was the most precise. The results summary is given Table 2.

It could be observed in Table 2 that RMSE decreased with the increasing k value. The most precise results were obtained for the k value being equal to 0.9. In this case, the root mean square error was 0.992% of RH and 0.566°C. From the point of view of thermal performance, even better results were obtained assuming k = 0.8 (RMSE = 0.553152°C), but the RMSE of relative humidity was ∼3.7% which was much higher than that in the previous case. Therefore, k = 0.9 could be considered as the best solution with respect to the overall hygrothermal performance of the block. Lichtenecker’s mixing rule cannot account for the geometry of the components in any rigorous way, in general. Therefore, there is not any exact relationship between the k value and the structure. However, the two limits of k = −1 and k = 1 correspond to Wiener’s bounds; that is, they express a finely layered system oriented perpendicular or parallel, respectively, to the direction of heat and moisture flux. Regarding the spatial arrangement of the analyzed block, of which geometry strongly evokes layering, the optimal value of k = 0.9 seems to be logical.

According to the same evaluation principles, the multiscale method provided the average difference of temperature 0.629°C and relative humidity 1.502%. These numbers can also be considered to be very small, but they are higher than the best solution provided by the materials homogenization technique. The differences between the results obtained by the domain decomposition method and the particular homogenization techniques are presented in a more detailed way in Figures 6 and 7 for selected temperature and relative humidity profiles covering the four seasons of the year. Here, CH stands for computational homogenization and MH for materials homogenization.

The distance 0.525–0.550 m corresponds to the exterior plaster, which explains the step changes of the relative humidity profiles in Figure 7 as the plaster exhibits the fastest response to the changing boundary conditions on the exterior side of the studied wall. The correspondence of results of the multiscale method with those of the domain decomposition method was found mainly in warmer and drier periods of the year, when the ambient air had lower relative humidity. The highest differences between relative humidity values could be observed on the material interface between the brick and the mortar near the external surface, where the RH values often reached nearly 100%. Since in the two-level analysis, averaging techniques were used in computation of effective material conductivity and capacity matrices, nonlinear effects of material properties and also nonlinearity in the distribution of relative humidity took place. Moreover, the aggregation of elements at the macro-level resulted in averaged values of unknowns and their gradients which were sent to the meso-level. It caused the values of relative humidity obtained by the multiscale method to be smoother and more averaged. Such an effect could not be totally avoided, but the decomposition to smaller aggregates could eliminate such problems in detached zones of the macro model. However, despite the above-mentioned difficulties or inaccuracies, the multiscale method with two-level analysis seems to be an effective tool for modeling the hygrothermal performance of building structures. A great potential can be seen in the solution of real-world masonry and other heterogeneous structures that cannot be transformed to a simple, one-dimensional problem. Additionally, the materials homogenization technique does not have a general validity, and the k parameter must be found for a geometric configuration of each problem solved. This could be considered as the biggest disadvantage of this technique as opposed to the multiscale method which could be applied universally.

The summary of computational times is presented in Table 3. Since the preparation work of materials homogenization regarding the determination of all the effective transport and storage parameters of the ceramic block for k = {−1.0; −0.9; −0.8; … 0.8; 0.9; 1.0} is also very time-consuming, the time requirements presented in the table can be misleading. However, all the works have been done within this comparative research and the optimal k parameter has been determined, so the advantage of the materials homogenization technique and the significantly reduced time requirements can be fully exploited when any further analyses of this very brick will be performed in the future, the long-term simulations of hygrothermal performance in particular. For instance, the computational simulation of 20-year performance using parallel computing would take 75 days, while after transformation and one-processor solution, it would be reduced only to less than 40 minutes. Such a time reduction makes possible to solve relatively easily a broad range of practical problems, e.g., the optimization of service life of brick block-based building envelopes by means of computer-aided design of compatible materials (plasters and thermal insulations) or the long-term assessment of brick blocks performance under different climatic conditions.

5. Conclusions

The presented research was aimed at the assessment of hygrothermal performance of ceramic brick blocks, which can be considered as a typical example of a complex computing task in the field of building physics. Since this assessment is computationally very expensive and time-demanding at the same time, as the complex geometry of the block has to be preserved, various simplification techniques were applied and their contributions were analyzed and discussed. The evaluation of these techniques was done by means of comparison of results with the performance obtained by the domain decomposition method.

Being supposed to reduce demands on time and computing power, two different principles of simplifications were assumed, namely, the materials homogenization represented by Lichtenecker’s mixing rule and computational homogenization by means of the multiscale method. While the computational homogenization principles are based on the solution of the problem on meso-level and macro-level using the master/slave computing system, the materials homogenization utilizes the effective media theory for the estimation of effective (homogenized) material parameters of the block. This enables a problem transformation to one-dimension, which is much less demanding on computing power.

The best results were achieved using the materials homogenization, assuming the k parameter of Lichtenecker’s mixing rule being equal to 0.9. In this case, the RMSE of relative humidity and temperature distribution over the reference year was only 0.992% RH and 0.566°C, respectively, when compared to results of the domain decomposition method. It proved that a sufficient level of accuracy was preserved. These results were even better than those provided by computational homogenization (1.502% RH and 0.629°C) which is still computationally more expensive. Besides the substantial time-saving (more than 99%), the materials homogenization allows also the solution on single-processor computers; thus, it is available for a major part of research departments. The latter does not although mean that the computational homogenization technique presented in this paper provides wrong results. The somewhat higher demands on computing time only make its application effective rather for the solution of more complex problems, where the 1D transformation is impossible. Critical construction details, real masonry, or construction parts may be stated in that respect. On the contrary, materials homogenization based on the effective media theory seems to be more advantageous for the analyses of heterogeneous building elements in the case of their development, innovation, or long-term assessment. Additionally, the presented technique is suitable to be applied on any building material with the periodical nonhomogeneous structure that can be homogenized. The application on more unstructured materials should be considered individually, but it probably would be questionable unless a periodical unit would be identified.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the Czech Science Foundation, under project no. 17-01365S.