{'fraction_non_alphanumeric': 0.03714715216301988, 'fraction_numerical': 0.015734090529480542, 'mean_word_length': 4.239978995273937, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CFD results calibration from sparse sensor observations with a case study for indoor thermal map

Current CFD calibration work has mainly focused on the CFD model calibration. However no known work has considered the calibration of the CFD results. In this paper, we take inspiration from the image editing problem to develop a methodology to calibrate CFD simulation results based on sparse sensor observations. We formulate the calibration of CFD results as an optimization problem. The cost function consists of two terms. One term guarantees a good local adjustment of the simulation results based on the sparse sensor observations. The other term transmits the adjustment from local regions around sensing locations to the global domain. The proposed method can enhance the CFD simulation results while preserving the overall original profile. An experiment in an air-conditioned room was implemented to verify the effectiveness of the proposed method. In the experiment, four sensor observations were used to calibrate a simulated thermal map with 167x365 data points. The experimental results show that the proposed method is effective and practical.


I. INTRODUCTION
The estimation of the indoor physical field (e.g., thermal map, airflow pattern, and pressure field) is an important problem in air-conditioning system design and indoor architectural-design. In recent years, the computational fluid dynamics (CFD) simulation has been widely used in building design to predict the thermal map, pollution dispersion, and ventilation performance [1], [2]. The accuracy of the CFD simulation strongly depends on the setting of boundary conditions and CFD models [3]. In practice, the real boundary conditions of the rooms or buildings are unavailable and the idealized version is used instead. Even though CFD modelling has been extensively researched, the reliability of the results remains a key concern when performing CFD simulations.
Currently, in order to produce reliable simulation results, we need to properly set the boundary conditions and adjust the computational model (numerical and/or physical models) input parameters to amend the agreement between the simulated results and the corresponding experimental data [4]- [6]. The current CFD model calibration work has mainly focused on minimizing the mismatch between simulated and experimental results. Tens to hundreds of simulation runs are required for the calibration process of finding acceptable model parameters. However, even a well calibrated CFD model cannot remove the gap between simulated and experimental data, which implies that the results of calibrated CFD models still have room for further improvement.
CFD simulations can map the sampled spatial domain to the airflow-field data. It can provide the global information of airflow fields of interest, which is rather expensive or even impossible to obtain from real experimental observations. On the other hand, sensors (e.g., thermocouples) can only practically provide sparse and isolated observations, but the observations are usually much more accurate than the simulation results. Can we calibrate/correct the simulated physical fields with sparse sensor observations? In this paper, we develop a methodology to solve this problem. To the best of our knowledge, no reported publication has focused on this problem.
The recent work [4]- [6] and the references therein mainly focused on the calibration of a CFD model. In this paper, we directly calibrate the results of CFD simulations. The proposed methodology can be applied for general CFD simulations and be a valuable additional step of the standard CFD simulation procedure.
In this work, we propose a methodology to calibrate the physical fields obtained from CFD simulations using sparse sensor observations. Since sensor observations are much more accurate than simulation results, we ignore the difference between sensor observations and the real values and use them as the ground truth. In addition, we assume that the simulation results can provide the rough global profile of the real physical fields, i.e., the profile of the simulated physical field is similar to that of the real physical field. The authors have taken inspiration from the image editing technique [7], [8]. As shown in Fig. 1, an image can be globally edited with sparse control samples. The right image is fused by the original image and the sparse control samples. The result of image editing reserves the content of the original image but with the style of the sparse control samples. For the calibration problem of simulated physical field, we can view the sparse sensor observations as the control samples. The calibration problem can be viewed as a simulated physical field editing problem. We formulated the simulated physical field calibration as an optimization problem. The objective function is a sum of two parts. One part is responsible to narrow the difference between the sensor observations and the simulation results, while the other part guarantees the similarity between the simulated results and the calibrated results. The optimization problem is solved by a numerical method. The proposed method was tested in an air-conditioned room. With the help of four sensor observations, indoor thermal map obtained from CFD simulation was well calibrated.
We summarize the contributions of this work as: 1) we proposed a formulation for the calibration of simulated physical fields from sparse sensor observations; 2) we formulated the proposed calibration problem as an optimization problem; 3) we provided a computationally efficient method to solve the optimization problem, and 4) we verified the effectiveness of the proposed methodology in an airconditioned room.
The rest of this paper is organized as follows. In Section II, we formulate the CFD-result calibration problem. In Section III, we detailed the proposed methodology in which we first formulate the calibration problem as an optimization problem, then present a numerical method to solve it, and finally we provide a computationally efficient approximated solution. In Section IV, we calibrated a thermal map of an air-conditioned room to verify the effectiveness of the proposed method. The concluding remarks are given in Section VI.
II. PROBLEM STATEMENT We denote a real physical field (e.g., thermal map, airflow pattern, and pressure field) of interest by f (x), where x = (x, y, x) ∈ Ω with Ω being the domain of interest. The physical field obtained from CFD simulation is a numerical solution, which can be simply interpolated to be a continuous one. We denote a continuous physical field interpolated from a CFD solution by f c (x), which can be described as where v(x) represents the error of the simulated physical field. We assume that max x∈Ω v(x) ≤ v max , and v max is a small value compared with the magnitude of f (x). We denote one of the sensor observations of the physical field by where n k is the measurement noise of the k-th observation, Ω s is the set of all sensing locations, and m is the number of available sensor observations. Without loss of generality, we assume that |n k | |v(x k )| for all 1 ≤ k ≤ m.
In practice, f c (x), x ∈ Ω and s k , 1 ≤ k ≤ m are known. In what follows, we estimate v(x), x ∈ Ω based on f c (x), x ∈ Ω and s k , 1 ≤ k ≤ m. With the estimated v(x), we can easily correct the simulated physical field f c (x).

A. Formulation of the optimization problem
Since |n k | |v(x k )| for all 1 ≤ k ≤ m, we ignore the measurement noise. We denote the error of the simulated physical field at the sensing locations by e k , and from (1) and (2) we can obtain In the remainder of this paper, we ignore the measurement noise and view e k as the real error of f c (x k ).
To estimate the error of the simulated physical field f c (x), i.e., v(x), we minimize the following cost function where is a priori known weight, which can be directly calculated from the simulated physical field. Note that w(x, y) represents the affinity between v(x) and v(y). Both σ m and σ d are positive variances. It is obvious that w(x, y) ∈ [0, 1] and that w(x, y) is very small unless the locations x and y are near and their corresponding physical value f c (x) and f c (y) are similar. Solving the minimization problem (3), we can find an estimation v(x), which we denote byv(x). Witĥ v(x), we can calibrate the simulated physical field f c (x), and the calibration result iŝ The cost function (3) is the summation of two parts. The first part is responsible for the regions around sensing locations. In the region around x k , the value of v(x) is forced to be similar to e k , especially when f c (x) is similar to f c (x k ). Consequently, in the region around x k , the physical value is calibrated to be similar to the sensor observation s k , especially for the subregion in which f c (x) is similar to f c (x k ). In this way, in the region around x k ,f (x) preserves the profile of f c (x) but the exact values are adjusted to be similar to the sensor observation s k . The first term in (3) can guarantee a good local calibration for the simulated physical field, but it has insignificant influence on the regions far away from x k for all On the other hand, the second term has a global influence, which can be explained from the following two perspectives: 1) the second term can transmit the adjustment from the local regions around the sensing locations to the whole physical field; 2) the second term can guarantee that in any local region, the gradient ofv(x) is small if in that region the simulated values f c (x) are similar, which guarantees the profile of f (x) is similar to that of the simulated physical field f c (x).
The optimization problem is controlled by three parameters: λ, σ m , and σ d . λ, which we call the balance factor, can balance the contribution of the two terms in (3). A small λ leads to good calibration in the local regions around sensing locations, while a large λ can globally reduce the gradient of the adjustment at the cost of the calibration accuracy in the local regions around the sensing locations.
w(x, y) is the product of two Gaussian functions. σ m and σ d control the rate of decay of the two gaussian functions, respectively. σ m , which we call the magnitude variance, mainly depends on the dynamic range of the simulated physical field. If the local gradient of the physical field is large, a large σ m is proper. Otherwise the first Gaussian function in (4) will be very small which may dilute the influence of the second Gaussian function. On the other hand, if the thermal map has a low dynamic range, a relatively small σ m is proper.
σ d , which we call the distance variance, controls the scope of the local region around one of the sensing locations in which the physical data can be calibrated by minimizing the first term in (3). In other words, in the particular local region, the sensor observation have a significant influence on the physical field calibration.

B. Numerical method to solve the optimization problem
It is difficult for us to directly minimize J 0 in (3) to find the estimation of the error, i.e., v(x), x ∈ Ω. Numerical methods can make the optimization problem easily solvable.
We denote the discrete space of the domain Ω in the CFD simulation by Ω D = {y 1 , y 2 , ..., y N }, where N is the number of mesh points. We assume that Ω s ⊂ Ω D , which implies that each sensing location is constrained to be one of the mesh points. Hence, we can denote x k = y l k and Ω s = {y l 1 , y l 2 , ...y lm }, where l k is the mesh index of the k-th sensing location. Then, we can approximate the cost function J 0 in (3) as i.e., w ij is the entry of W located at the i-th row and j-th column. For a more elegant description,we use the following cost function J instead of J 1 where v = [v 1 , v 2 , ..., v N ] T , 1 ∈ R N is a vector with all 1 entries, W D j = diag{w j }, w j is the j-th column of W, and A j is a matrix with all 1 entries in j-th column and all 0 entries in others. Here, diag{·} and Tr(·) are diagonal matrix and trace operators, respectively. It is clear that To minimize the cost function J, we take the derivative of (7) with respective to v and set it as zero, i.e., dJ dv = 2(Hv − b) = 0 (9) From (9), we can easily obtainv Substitutingv into (5), we obtain the calibrated physical field For the matrix H ∈ R N ×N , N is the number of mesh points. Generally, N is very large and therefore, the computational cost of solving H −1 in (10) is extremely expensive. In practice, we cannot directly obtain v from (10). Thus, we provide a computationally efficient method to find the approximated solution of v.

C. Computational efficient approximated solution
To avoid directly solving the inverse of the matrix H, which has very large dimension, we provide an approximated solution of the inverse problem. A similar method has been previously applied for image editing [7].
Considering the definition of w ij , i.e., w ij = w(y i , y j ), we find that each column of W represents the affinities between the point ( whose mesh index is the column index) and the whole physical field. The neighbouring points have similar affinities with other points. Hence, the neighbouring columns have similar entries, especially for the entries far away from the diagonal, which implies that the neighbouring columns of W have a large correlation coefficient. Therefore, the matrix W has many near zero eigenvalues, and it can be approximated by a low rank matrix.
Since the matrix W s = 2W ∈ R N ×N , W s can also be approximated by a low rank matrix. We rewrite the symmetric matrix W s as We assume that the number of dominant eigenvalues of W s is n. Using the eigen decomposition and ignoring all the near zero eigenvalues, we can approximate W s as where U is an orthogonal matrix, Λ is a diagonal matrix with all eigenvalues of W s in the diagonal, Λ 1 = Λ(1 : n, 1 : n) consists of all dominant eigenvalues of W s in its diagonal, and V ∈ R N ×n has exactly rank n. We denote V by where X and Y have the same sizes as A and B in (12), respectively. Then, we can obtain that Comparing (13) and (12), as illustrated in Fig. 2, we can find Substituting (14) into (12) yields is the first n column of W s . Using this way, we can find a low rank approximation of W s , which can both resolve the storage and computational cost problems in solving H −1 .
Substituting the approximation of W s in (15) into (8) and using the Woodbury formula [9], we can obtain Then, substituting (16) into (10) yieldŝ In (17), we need to solve the inverse of two matrices. D is a diagonal matrix and its inverse can be easily solved. Because (−A + Z T D −1 Z) ∈ R n×n and n N , solving (−A + Z T D −1 Z) −1 is computationally much more efficient.
Substituting (17) into (11), we can obtain the calibrated physical field. Next, we show an experiment to support the proposed methodology.

A. Statement of the experiment
We tested the proposed methodology for a 7.28m×3.32m×2.5m room, which is shown in Fig. 3. In the room, an Active Chilled Beam system (ACBs) was installed for air conditioning, as shown on the ceiling of the top picture in Fig. 3. One air outlet is at the corner of the room. The walls and the door are made by thermal insulation materials, and the windows are double glazing. Hence, in the CFD simulation, we simply set the walls, the door, and the windows as insulated walls.
temperature sensors @ z=1.2m temperature sensor @ z=2.5m velocity sensor @ z=2.5m As shown in Fig. 3, the room contains four heaters (0.76m×1.2m×0.02m). As shown in the bottom picture of Fig. 3, the left two heaters are 300W, and the other two are 400W. The height above the ground of the top surfaces of the heaters is 1.02m. The bottom faces of the heaters were made by insulated materials. One sensor was used to measure the velocity of the inlet cooling air, as shown in the leftmost picture of Fig. 4. Another sensor was used to measure the inlet cooling air temperature, which is shown in the middle picture in Fig. 4. The other eight sensors were used to measure the temperature of a horizontal plane 1.2m above the ground. Four sensors were hung and suspended just above centers of the four heaters, respectively. The four sensors are represented by stars (s1-s4) in the bottom picture of Fig.  3. One of the four is shown in the right-most picture of Fig. 4. Three sensors (s5-s7) were placed at the middle plane of the room. The sensor s8 was placed 20cm from the wall. All the eight sensors (s1-s8) were 1.2m above the ground. Fig. 4: Sensors used for indoor parameters measurement. Left: velocity sensor used to measure the inlet air velocity. Middle: temperature sensor used to measure the inlet cooling air temperature. Right: temperature sensor used to measure the temperature at 1.2m height, which is hung above the heater.
The experiment was carried out on 26th May, 2016. The four heaters and the air-condition system were turned on at 9:30 am. We collected the velocity and temperature of the inlet cooling air (via the sensors s v and s t , respectively) from 11:17 to 13:10. The data is shown in Fig. 5. The sensors s1-s7 recorded the temperature of the 1.2m plane from 9:30 to 13:10, the data for which is shown in Fig. 6. The sensor s8 does not store the sensor readings. We have recorded that after 11:00, the readings of sensor s8 remained at 25.6 • C.

B. CFD simulation
The simulation was performed using the commercial CFD software Ansys FLUENT v.13.0. The proposed calibration methodology has the potential to be applied for dynamic CFD simulations if both the real-time simulation results and sensor observations are available. In this work, however, for simplicity we only consider the steady state simulation results.
The indoor air was assumed as the ideal gas with density of 1.225 kg/m 3 . The viscous model was chosen as the realizable k − ε model. Since the test-bed was built in a laboratory, we did not consider the solar radiation. The door, ground, ceiling, windows, and walls are assumed to be insulated walls. During the experiment, all lights were turned off. For the velocity and temperature of the inlet cooling air, we use the mean values of the sensor readings in the time period 12:10 to 13:10 in Fig. 5. The bottom faces of the heaters were made by insulated material, but they were not completely insulated. We assumed that 5% of energy was transferred via the bottom faces of the heaters. Corresponding to the four sensors s1-s4, we denote the four heaters by heater 1, heater 2, heater 3 and heater 4, respectively. The boundary conditions of the CFD simulation are detailed in Table I.  Ti− inlet cooling air temperature; Vin− inlet cooling air speed; θ− the angle of Vin with the ceiling; To− outlet air temperature; P relative − pressure relative to the reference pressure (P reference = 101325Pa); hc− convective heat flux.
We created the structured mesh using Gambit v2.4. Since all the walls are insulated, we did not consider the boundary layer of the walls and uniformly divided the horizontal plane (X − Y plane) into 2cm × 2cm squares. In the z-direction, we made the mesh dense near the ceiling (z = 2.5m) and heater top faces (z = 1.02m), and divided the entire 2.5m height into 137 intervals. The total mesh size is ∆ 1 = 167 × 365 × 137 − 4 × 37 × 59 × 1 = 8, 342, 103.
To verify the grid independency, we made the mesh dense and increased the mesh size to ∆ 2 = 18, 527, 205. We did CFD simulations for both mesh files. The results were found to be in good agreement. Therefore, the mesh file with 8,342,103 cells is proper and the number of cells is large enough for our experiment.
The convergence criteria were set as 10 −6 for energy and 10 −3 for all the other parameters. As shown in Fig. 7, the simulation converged with 442 iterations. We are interested in the thermal map of the room at 1.2m height. The thermal map obtained from the CFD simulation is shown in Fig. 8a.

C. The calibration results
We use the sensor observations obtained from sensors s1-s4 to calibrate the thermal map in Fig. 8a. The four sensors s1-s4 were just above the centers of the four heaters, respectively.
We set σ m = 1000 and σ d = 1. The thermal map in Fig. 8a consists of 60955(= 167 × 365) temperature values of all the mesh points at the plane z = 1.2m. The balance factor λ is related to the number of mesh points and the number of sensor observations. We set where 0 < α ≤ 1. The first part in (6) is the summation of m × N terms while the second part is the summation of N ×N terms. If we set α = 1, it means that we have the same trust on the sensor observations and the simulation results on all mesh points. Obviously, sensor observations are more accurate. We set α = 0.01 and the calibrated thermal map is shown in Fig. 8b. Here, α = 0.01 implies that the influence of each sensor observation in (6) is 100 times important than the simulated data on each mesh point. Comparing Fig. 8b with Fig. 8a, we can easily find that the calibrated thermal map preserves the profile of the CFD results. However, Fig.8b cannot show how the simulated thermal map is calibrated. We show the estimated error (i.e.,v) in Fig. 9a, which is the difference between Fig. 8a and Fig. 8b. Considering (11), we can obtain the simulated thermal map by subtracting the estimated error from the simulated thermal map.

D. The influence of balance factor
To check the effect of the balance factor, we set α as 0.01, 0.1, 0.5, and 1 respectively. The corresponding estimated errors of the thermal map in Fig. 8a is given in Fig 9. The observations of sensors s1-s4 are used to estimate the error. To show the effectiveness of the proposed method, we use the observations of sensor s5-s8 as the ground truth to test the estimation performance. The error of the CFD simulation results and the calibrated thermal map at the 8 sensing locations are shown in Fig. 10. Fig. 9 shows that with different balance factors, the estimated errors (i.e.,v) almost have the same profiles. However, the magnitude of the estimated errors decrease when the balance factor increases, which implies that the influence of sensor observations decreases and we can only provide a very small adjustment for the original simulated thermal map. In other words, the balance factor can control the magnitude of the estimated errors, i.e., the magnitude of the adjustment for the simulated thermal map. As shown in Fig. 9 and Fig. 10, with a smaller balance factor, the magnitude of the estimated error is larger and we can obtain a better performance.
The root mean square error (RMSE) of the thermal map data at sensing locations are given in Table  II. Both Fig. 10 and Table II show that we can enhance the thermal map obtained from CFD simulation with the four sensor observations (s1-s4). With the calibration, the RMSE of the data at sensing locations  s5-s8 can be improved by 33.38%. However, we can see from Fig. 10 that the errors of the regions around s6 and s7 are still large. The CFD simulation results around s6 and s7 have large errors, i.e., around 2 • C and 1.5 • C, respectively. The proposed method reduces the errors to around 1.5 • C and 1 • C, which are still large. To address this issue, we should consider the capacity of the proposed method. The largest magnitude to which the proposed method can adjust the CFD results almost equals the maximum error of the CFD results at all sensing locations where the sensor observations were used to calibrate the CFD results. It can be seen in Fig. 10 that the maximum error of the CFD results at sensing locations s1-s4 is around 1.3 • C. Therefore, the maximum value of the estimated error of the simulated thermal map is less than 1.3 • C. The balance  factor can control the magnitude of the maximum value of the estimated errors. In other words, the balance factor can control the magnitude of the adjustment for the simulated thermal map, but the upper bound of this magnitude is determined by the sensing locations.

E. The influence of magnitude variance
To check the influence of magnitude variance, we set σ m as 1000, 100, 10, and 1, respectively. With the four sensor observations (s1-s4), we estimated the error of the thermal map in Fig. 8a. The estimated errors are given in Fig.11. The error of the estimated thermal map at the 8 sensing locations are shown in Fig. 12, and the RMSE of the calibrated thermal map at the 8 sensing locations are given in Table III. Fig. 11 shows that with different magnitude variances, the maximum magnitude of the estimated errors (i.e.,v) are almost the same, which implies that different with the balance factor, the magnitude variance has no significant influence on the the maximum magnitude of the estimated errors.
It can be seen from Fig. 11c-11d that the region in the two red circle is nearer to sensor s4 than s3 but the estimated error is similar with the estimated error around s3. We can see form Fig. 8a that the temperature values of the thermal map in the region marked by the two circles (in Fig. 11c-11d) are similar to those in the region around sensor s3. The proposed thermal map calibration work is a process of fusing sensor information and the information of CFD results. A smaller magnitude variance can enhance the influence of the CFD results in this information fusion work.
In addition, Fig. 11d shows that if σ m = 1 in the regions far away the sensing locations s1-s4, the estimated errors are almost zero. In other words, the proposed method has no influence for these regions. As we mentioned before, a small magnitude variance leads to a small value of the first Gaussian function in (4), which can dilute the influence of the second Gaussian function. The smaller the magnitude variance σ m , the smaller the weight w(x k , y), especially when y is far away the sensing location x k . Therefore, if the magnitude variance is too small, the proposed method could not globally calibrated the thermal map.  Fig. 12 shows that with the decrease of the magnitude variance, the errors around sensing locations s1-s4 also decrease but those around sensing locations s5-s8 do not. From Table III we can find that the RMSE at sensing locations s5-s8 increases with respect to the decrease of the magnitude variance. If we reduce the magnitude variance to 1, the calibration work makes the accuracy of the simulated thermal map worse. Hence, we conclude that reducing the magnitude variance can improve the local calibration but has no significant help for global calibration. If enough sensor observations are available we can set a small magnitude variance unless a large one is proper.

F. The influence of distance variance
To check the influence of distance variance, we set σ d as 1, 0.5, 0.2, and 0.1, respectively. With the four sensor observations (s1-s4), we estimated the error of the thermal map in Fig. 8a. The estimated errors are given in Fig.13. The error of the estimated thermal map at the 8 sensing locations are shown in Fig. 14, and the RMSE of the calibrated thermal map at the 8 sensing locations are given in Table IV.  From Fig. 13 we can find that like the magnitude variance, the distance variance has no significant influence on the maximum value of the estimated errors. It is clear shown in Fig. 13 that with a smaller distance variance, the scope of non-zero region of the estimated errors dramatically decrease. Fig. 14 shows that when the distance variance σ d ≤ 0.5, the proposed method can provide almost perfect calibration at the four sensing locations (s1-s4). Table IV shows that with a smaller distance variance, the errors of the calibrated thermal map at the sensing locations s5-s8 increase. These agree with our previous analysis that the distance variance controls the scope of the region on which the sensor observation has significant influence. A small distance variance σ d can lead to very good local calibration while a large σ d can provide global influence for an isolated sensor observation. V. CONCLUDING REMARKS CFD simulation results calibration from sparse sensor observations is a new and interesting problem. We formulated this problem as an optimization problem. The proposed method can provide an effective solution for this problem. The experimental results show that with four sensor observations (s1-s4), the proposed method has 33.38% improvement on the accuracy of the simulation data around sensors s5-s8. If we consider the whole thermal map (including the region around sensors s1-s4), the improvement can be even larger.
We need to set three parameters for the proposed method: balance factor, magnitude variance, and distance variance. The balance factor controls the magnitude of the adjustment for the simulated thermal map, but the upper bound of this magnitude is determined by the sensing locations. The distance variance controls the scope of the regions on which an isolated sensor observation has significant influence. A small distance variance can lead to very good local calibration while a large one can provide global influence for an isolated sensor observation. The magnitude variance controls the influence of the CFD results and scales the influence of the distance variance. Reducing the magnitude variance can improve the local calibration but has no significant help for global calibration. If enough sensor observations are available we can set both small magnitude variance and small distance variance. If the number of available sensor observations is very limited, we should set relatively larger values for the two variances to achieve global performance.
In addition, the capacity and the performance of the proposed method is closely related to the number of sensor observations and the sensing locations. The issue of how many sensors are required and where to place them are very interesting topics for future research.