Sensitivity Analysis of the Effects of Building Morphology on Indoor–Outdoor Wind–Light–Thermal Coupling Processes

2.1 Building morphology parameterization and morphological non-uniformity calculation

2.1.1 Selection and definition of morphological parameters

The accurate selection and quantification of morphological parameters are the basis for achieving an accurate relationship between building morphology and the thermodynamic irreversibility of indoor–outdoor wind–light–thermal coupling processes. This study breaks through the limitations of traditional studies that use macroscopic and coarse parameters such as shape coefficient and surface area, and selects five types of core morphological parameters that can be directly regulated in building design and have significant effects on the coupling process. All parameters are described by continuous variables, taking into account both the engineering practicality and optimization feasibility of the parameters, thereby achieving refined quantification of building morphology. The aspect ratio, as the core parameter describing the overall geometric proportion of a building, is defined as the ratio of building width W, depth D, and height H, and the quantitative expression is λ = W:D:H. The value range is set as 1:0.8:1–1:1.2:1. This range covers the common proportions of high-density urban buildings, and can accurately reflect the influence of overall geometric characteristics on airflow bypass patterns and solar radiation reception efficiency, avoiding the defect that traditional single proportion parameters cannot comprehensively characterize the three-dimensional geometric characteristics of buildings. The concave–convex depth of the south-facing exterior wall δ is defined as the maximum protrusion of the concave–convex part, with a value range of 0–0.5 m and an interval of 0.05 m. The interval setting is derived based on the critical size of near-wall flow field separation, ensuring that the refined effects of different concave–convex degrees on near-wall vortex generation, airflow separation, and radiation reception area can be captured, breaking through the limitation that traditional macroscopic parameters cannot characterize local geometric mutations. The roof inclination angle θ is the angle between the roof and the horizontal plane, with a value range of 0–30° and an interval of 3°. This range covers the commonly used roof inclination angle range in buildings. Its quantified values match the seasonal variation of solar altitude angle, directly determining the reception efficiency of direct solar radiation and the roof heat conduction path, providing a basis for the accurate calculation of radiation entropy generation rate in subsequent analysis. The facade opening ratio η is defined as the ratio of the facade window area to the total facade area, with a value range of 0–0.3 and an interval of 0.05. It is determined in combination with building energy-saving design specifications, and can accurately characterize the regulatory effect of window distribution on indoor–outdoor ventilation efficiency, radiation heat gain, and heat exchange processes, avoiding the quantification deviation of heat exchange caused by the coarse values of traditional opening ratio. The eave overhang length L is the length by which the eave protrudes from the exterior wall of the building, with a value range of 0–1 m and an interval of 0.1 m. It is set based on the shading requirements of summer solar altitude angle and the disturbance law of near-wall airflow, which can not only reflect its shading effect on summer solar radiation to reduce indoor heat gain, but also capture its disturbance effect on near-wall airflow distribution, providing support for the refined analysis of viscous dissipation entropy generation rate in the flow field. The continuous values and reasonable interval design of all parameters provide standardized and refined geometric inputs for subsequent multiphysics weak coupling simulation, while also laying the foundation for the study of continuous parameter variation and interaction effects in global sensitivity analysis, thereby achieving an accurate relationship between morphological parameters and thermodynamic irreversibility. As shown in Figure 1, this study establishes a parametric description framework for building core morphology and realizes the generation of morphological non-uniformity potential function.

Figure 1. Schematic diagram of building core morphology parameterization and morphological non-uniformity potential function generation

2.1.2 Construction and calculation of the morphological non-uniformity index (ξ)

To break through the limitation that traditional macroscopic morphological parameters cannot capture the refined effects of local geometric mutations on wind–light–thermal coupling processes, this study constructs morphological non-uniformity ξ as a scalar evaluation index for characterizing building morphology complexity based on the concept of “friction potential” in non-equilibrium thermodynamics. Its core lies in quantifying the spatial variation characteristics of building morphology through the distance potential function, thereby achieving an accurate relationship between local morphological mutations and thermodynamic characteristics. Let the position vector of any spatial point in the computational domain be r, and define the distance potential function φ(r) as the shortest distance from this point to the nearest building surface. This function can intuitively reflect the spatial fluctuations and local mutation characteristics of building morphology, overcoming the defect that traditional macroscopic parameters can only describe the overall contour. The level set method is used to solve the spatial distribution of φ(r), and its governing equation is:

$\frac{\partial \varphi}{\partial t}+u \cdot \nabla \varphi=0$                          (1)

where, u is the virtual velocity field, and its value is set as 1/10 of the characteristic velocity of the building. This value has been verified through numerical stability analysis, which can effectively avoid numerical oscillation during the evolution process of the level set function and ensure solution convergence. The finite volume method is used to discretize the governing equation, and the grid size is set as 1/50 of the building characteristic dimension. This size is determined based on the characteristic scale of local geometric mutations, which can accurately capture the spatial variations of local details such as exterior wall concave–convex features and eave overhangs, significantly improving the calculation accuracy of φ(r), and laying the foundation for the accurate calculation of subsequent morphological non-uniformity.

Morphological non-uniformity ξ is obtained by integrating and normalizing the squared gradient of the distance potential function over the entire computational domain. Its core lies in quantifying the intensity of spatial morphological variation through the integral of the squared gradient. A larger gradient indicates a more significant local morphological mutation and a stronger potential influence on thermodynamic irreversible losses. The calculation formula is:

$\xi=\frac{1}{V} \iiint_V|\nabla \varphi|^2 d V$                     (2)

where, V is the total volume of the computational domain. The integral is solved by the Gaussian integration method. The computational domain is divided into several tetrahedral elements. Linear interpolation is used within each element to approximate the gradient distribution of φ(r), and the global integral result is obtained through element integration accumulation. Numerical verification shows that this method can control the integral error within 3%, ensuring the calculation accuracy of ξ. To support subsequent multiphysics simulation and global sensitivity analysis, the Latin hypercube sampling method is used to generate 64 independent morphological configurations within the reasonable value ranges of the five types of morphological parameters. Among them, 48 are used for multiphysics simulation and PCE surrogate model training, and 16 are used for model validation and method effectiveness verification. During the sampling process, parameter distribution uniformity verification is carried out to avoid distortion in subsequent analysis caused by sampling bias, while ensuring that the sampled configurations can comprehensively cover the variation range of morphological parameters and fully characterize the thermodynamic response characteristics under different morphology complexities, thereby achieving seamless connection between morphological non-uniformity and subsequent multiphysics simulation and sensitivity analysis.

2.2 Multiphysics weak coupling simulation framework lattice Boltzmann method–monte Carlo–finite element method

2.2.1 Simulation domain and boundary condition settings

The accurate setting of the simulation domain and boundary conditions is the basis for improving the simulation accuracy of wind–light–thermal coupling and ensuring the reliability of subsequent entropy generation rate calculations. This study breaks through the limitations of unreasonable simulation domain size and coarse grid discretization in traditional studies, and constructs a standardized and refined simulation domain system in combination with the physical characteristics of building wind–light–thermal coupling processes. A rectangular computational domain is adopted, covering the building block and the surrounding atmospheric boundary layer. The overall structure of the indoor-outdoor wind-light-thermal coupled model and the setting rules of multi-physical field boundaries are depicted in Figure 2. The size of the computational domain is set according to building height (H) as 10H (x direction), 8H (y direction), and 5H (z direction). Numerical verification shows that this size can ensure the full development of the flow field around the building, no boundary reflection interference in the light field radiation, and stable temperature distribution in the thermal field, effectively avoiding the distortion problem of boundary effects caused by insufficient computational domain size in traditional studies. The interior of the building is divided into rooms according to actual functions, and discretized using structured grids, with grid sizes selected as 0.1–0.3 m, taking into account both computational efficiency and accuracy. For local geometric mutation regions such as exterior wall concave–convex features and eave overhangs on the building surface, the grid size is refined to 0.05 m. This refinement size is determined based on the characteristic scale of local flow field separation and radiation energy exchange, which can accurately capture near-wall flow field details and differences in surface radiation distribution, providing high-quality grid support for the refined calculation of subsequent local entropy generation rate (especially viscous dissipation and radiation absorption entropy generation rate), and breaking through the limitation that traditional uniform grids cannot characterize local physical quantity variations.

Figure 2. Schematic diagram of indoor–outdoor wind–light–thermal coupled spatial physical model and multi-field boundary coordination

The key to boundary condition setting is to achieve coordinated matching and accurate quantification of the three-field boundaries. By avoiding simulation errors caused by the disconnection of traditional boundary conditions and empirical parameter estimation, a reliable boundary foundation is established for multiphysics weak coupling iteration. The flow field boundary adopts an inlet condition of logarithmic wind profile that conforms to the actual atmospheric boundary layer, and its expression is:

$u(z)=\frac{u_*}{\kappa} \ln \left(\frac{z}{z_0}\right)$                                (3)

where, u* = 0.4 m/s is the friction velocity, κ = 0.41 is the von Kármán constant, and z₀ = 0.1 m is the surface roughness. This parameter combination has been verified by wind tunnel experiments and can accurately characterize the airflow distribution characteristics of high-density urban blocks. The outlet adopts a standard atmospheric pressure boundary, the building surface and ground are set as no-slip boundaries, and the top of the computational domain adopts a symmetry boundary, forming a complete and realistic flow field boundary system. The light field boundary is based on the standard Commission Internationale de l’Éclairage (CIE) sky model, and a band-based setting strategy is innovatively adopted. Solar radiation is divided into two bands: visible light (0.38–0.76 μm) and near-infrared (0.76–2.5 μm), with band weights set as 0.45 and 0.55, respectively, conforming to the actual solar radiation spectral distribution and ensuring the accuracy of radiation energy calculation. The building surface and ground are both set as diffuse reflection surfaces, with reflectivity set according to the actual characteristics of building materials. The side surfaces of the computational domain adopt blackbody boundaries to avoid radiation reflection interference, solving the problem of radiation heat gain calculation deviation caused by traditional light field boundaries ignoring spectral characteristics. The thermal field boundary adopts composite wall parameters to conform to engineering practice. A constant heat load is set indoors to simulate heat dissipation from occupants and equipment. The outdoor ambient temperature adopts hourly data from a typical meteorological year (TMY). The core is that the convective heat transfer coefficient is directly output by LBM flow field simulation, replacing traditional empirical formula estimation, effectively reducing the influence of heat transfer coefficient errors on thermal field simulation accuracy, achieving coordinated linkage between flow field and thermal field boundaries, and providing reliable temperature gradient data for the accurate calculation of heat conduction entropy generation rate.

2.2.2 Solution methods of each physical field

The innovation in the solution of the flow field and light field focuses on efficiency improvement and calculation accuracy optimization, aiming to break through the performance bottlenecks of traditional solution methods in coupling scenarios, and provide high-quality basic data support for subsequent entropy generation rate calculation. The flow field is solved using the multiple-relaxation-time LBM, and the D3Q19 three-dimensional 19-velocity model is selected. Its evolution equation is:

$f_i\left(x+e_i \Delta t, t+\Delta t\right)=f_i(x, t)+M^{-1} \mathrm{~S}\left[M\left(f_i^{e q}-f_i(x, t)\right)\right]$                       (4)

where, $f_i$ is the particle distribution function, $e_i$ is the velocity vector, $\Delta t=1 \mathrm{~s}$ is the time step, $M$ is the transformation matrix, and $S$ is the relaxation matrix. The core lies in optimizing the configuration of relaxation coefficients. The viscous relaxation coefficient, momentum relaxation coefficient, and energy relaxation coefficient are set as $1.0,0.8$, and 0.9 , respectively. Numerical stability verification shows that this configuration can effectively suppress numerical oscillation in flow field simulation and improve the solution accuracy of complex flow fields (such as near-wall vortices). At the same time, by directly solving the non-equilibrium stress tensor:

$\tau_{i j}=\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$                               (5)

The viscous dissipation entropy generation rate is calculated without additionally solving the Navier–Stokes equations, greatly reducing computational complexity and improving solution efficiency. The light field adopts the MC ray tracing method, and innovatively adopts an independent band-based solution strategy. For the two bands of visible light and near-infrared, 10⁷ rays are generated separately. The ray emission direction and energy follow the actual solar radiation distribution, and the radiation characteristics are accurately matched through random sampling methods. The interaction between rays and surfaces adopts the diffuse reflection model, accurately calculating the absorbed and reflected energy of each ray, and finally outputting the absorbed radiation power density $Q_{\text {rad}}$ of the surface grid elements of the building, which is directly used as the internal heat source term for thermal field solution, achieving accurate linkage between the light field and the thermal field. A spatial partition acceleration algorithm is introduced. By dividing spatial regions, ineffective interactions between rays and grids are reduced, significantly decreasing computational cost and solving the problem of low efficiency in traditional ray tracing methods.

The key to thermal field solution is the introduction of the non-Fourier heat conduction equation, fully considering the non-equilibrium thermal relaxation effect of building envelopes under short-term thermal disturbances, and breaking through the limitation that traditional Fourier equations are difficult to describe unsteady heat conduction, with high-accuracy solution achieved by the FEM. The non-Fourier heat conduction equation adopts the Cattaneo–Vernotte form, including the heat flux equation and the energy conservation equation. The specific expressions are:

$\tau_q \frac{\partial q}{\partial t}+q=-k \nabla T$                   (6)

$\rho c_p \frac{\partial T}{\partial t}=-\nabla \cdot q+Q_{r a d}$                                   (7)

where, $\tau_q=10$ is the thermal relaxation time, $q$ is the heat flux density, $k$ is the thermal conductivity, $T$ is the temperature, $\rho$ and $c_p$ are the material density and specific heat capacity, respectively, and $Q_{\text {rad}}$ is the absorbed radiation power density output from the light field. To improve the calculation accuracy of the temperature field, tetrahedral elements are selected for finite element discretization, and quadratic interpolation functions are used to fit the temperature distribution. The discretized linear equation system is solved by the conjugate gradient method, and a residual less than $10^{-6}$ is set as the convergence criterion, ensuring that the calculation error of the temperature gradient is controlled within an acceptable range. This solution method can not only accurately capture the unsteady heat conduction process of building envelopes, but also achieve deep coupling between the light field and thermal field by directly introducing $Q_{\text {rad}}$ output from the light field, providing high-quality temperature gradient data for the accurate calculation of heat conduction entropy generation rate. At the same time, it connects with the convective heat transfer coefficient output from the flow field, constructing a complete system for coordinated solution of the three fields.

2.2.3 Weak coupling strategy and program implementation

This study innovatively adopts a weak coupling strategy of “interactive transfer at each time step,” breaking through the limitations of low computational efficiency in traditional strong coupling and data disconnection in commercial software assembly, and constructing a closed-loop collaborative iterative mechanism of the flow field, light field, and thermal field, thereby achieving accurate characterization and efficient simulation of non-equilibrium thermodynamic processes. Figure 3 shows the logical framework of weak coupling, asynchronous collaborative iteration and data closed loop for the combined LBM, MC and FEM methods. The coupling logic takes the time step as the core link, and the time step is set as Δt = 1 s. This value takes into account both the evolution accuracy of physical processes and computational efficiency, and can accurately capture the dynamic variation characteristics of wind–light–thermal coupling processes. After initializing the initial parameters of the flow field, light field, and thermal field, within each time step, the flow field is first solved by Multiple-Relaxation-Time (MRT)–LBM, and the near-wall temperature T_wall and convective heat transfer coefficient h are output and directly transferred to the FEM thermal field solution module. At the same time, the light field is solved by MC ray tracing, and the absorbed radiation power density Q_rad on the building surface is output and directly injected into the energy conservation equation of FEM as the internal heat source term, achieving accurate linkage between the light field and the thermal field. Subsequently, the thermal field is solved by FEM, and the building surface and indoor–outdoor temperature field T are updated. The updated T_wall is then fed back to LBM as the wall boundary condition for the next time step, forming a closed loop of interactive iteration among the three fields. This coupling strategy does not require complex coupling interface design. It not only avoids the sharp increase in computational cost caused by synchronous solution of each physical field in strong coupling methods, but also solves the problems of data reading and writing delay and insufficient coupling accuracy in commercial software assembly. Finally, core parameters of each grid, including velocity u, temperature T, and radiation absorptivity α_λ, are output, providing high-quality data support for the accurate calculation of subsequent local entropy generation rate.

Figure 3. Logical diagram of Lattice Boltzmann Method (LBM)–Monte Carlo (MC)–Finite Element Method (FEM) weak coupling asynchronous collaborative iteration and data closed loop

The core of program implementation is multidimensional optimization design, seeking a balance between computational efficiency and simulation accuracy, thereby adapting to the large-scale simulation requirements of multiple morphological configurations. A core solver is written in C++ language, combined with Compute Unified Device Architecture (CUDA) parallel computing technology. Parallel optimization is mainly carried out for the LBM flow field solution and MC ray tracing modules with large computational cost. The number of parallel threads is set as 1024. This thread number has been verified by GPU computing power matching, which can maximize the utilization of hardware resources and control the 24-hour simulation time of a single morphological configuration within 4–6 hours, improving efficiency by more than 60% compared with traditional serial programs. An in-memory shared data interaction module is innovatively designed to realize direct memory data transfer among the three solution modules of LBM, MC, and FEM, without the need for data interaction through disk reading and writing, greatly reducing data delay and further improving coupling efficiency. At the same time, an adaptive grid refinement algorithm is embedded to identify key regions in real time, such as local geometric mutation regions on the building surface and flow field vortex regions, and automatically refine the grid size to 0.05 m, while maintaining the original grid size in non-critical regions. Under the premise of ensuring the calculation accuracy of local physical quantities, this avoids the waste of computational cost caused by global grid refinement, achieving a balance between accuracy and efficiency, and providing technical support for the efficient simulation of 64 morphological configurations in subsequent global sensitivity analysis.

2.3 Irreversibility evaluation system based on local entropy generation rate

2.3.1 Thermodynamic basis of entropy generation rate

Based on the second law of non-equilibrium thermodynamics, the core of this study is to break through the limitation that traditional entropy generation rate calculations only focus on the global total amount and cannot distinguish the contributions of different irreversible mechanisms, and to construct a local entropy generation rate evaluation basis suitable for the wind–light–thermal coupling system, clarifying the composition mechanism and physical meaning of the total entropy generation rate, and providing theoretical support for the accurate quantification and decomposition of subsequent irreversible losses. The spatial distribution characteristics of multiphysical irreversible effects and the decomposition mode of 3D local entropy generation rate are demonstrated in Figure 4. The local entropy generation rate $\dot{s}_{\text {gen}}(r, t)$, as the core physical quantity describing the intensity of irreversible processes, is essentially the entropy generated by the system per unit volume and per unit time, strictly following the Clausius inequality. Its core characteristic is that it is only caused by irreversible processes and is always non-negative. According to the physical characteristics of the wind–light–thermal coupling system, the total local entropy generation rate is innovatively decomposed into three independent components corresponding to the flow field, thermal field, and light field, respectively, achieving accurate separation of different irreversible mechanisms. Its expression is:

$\dot{S}_{\text {gen}}(r, t)=\dot{S}_{v i s c}+\dot{S}_{\text {cond}}+\dot{S}_{\text {rad}}$                  (8)

where, $\dot{s}_{\text {visc}}$ corresponds to the irreversible loss caused by viscous dissipation in the flow field, $\dot{s}_{\text {cond}}$ corresponds to the irreversible loss caused by heat conduction under finite temperature difference in the thermal field, and $\dot{s}_{\text {rad}}$ corresponds to the irreversible loss caused by solar radiation energy conversion in the light field. All three components satisfy the non-negative constraint. The innovative value of this decomposition equation lies in breaking through the limitation that traditional single entropy generation rate cannot quantify the contributions of different physical processes. It can accurately match the flow field, thermal field, and light field parameters output from the multiphysics weak coupling simulation described above, achieving separate quantification and comprehensive evaluation of losses in each irreversible process. It provides a unified and accurate response quantity for identifying the influence laws of different morphological parameters on each irreversible mechanism in subsequent global sensitivity analysis, ensuring the rigor and pertinence of thermodynamic analysis.

Figure 4. Spatial mapping diagram of multiphysics irreversible mechanisms and three-dimensional local entropy generation rate decomposition

2.3.2 Calculation methods of entropy generation rate components

For different irreversible mechanisms in the wind–light–thermal coupling system, this study innovatively designs accurate calculation methods for each entropy generation rate component, closely connecting with the output parameters of the multiphysics weak coupling simulation described above, breaking through the limitations of insufficient accuracy and poor adaptability of traditional calculation methods, and achieving separate quantification of each irreversible loss. The viscous dissipation entropy generation rate is caused by fluid viscous friction and vortex motion in the flow field, and is the core irreversible loss of mechanical energy conversion into thermal energy. The innovation in its calculation lies in directly using the non-equilibrium stress tensor output by MRT–LBM, without additionally solving the Navier–Stokes equations, greatly improving calculation efficiency and accuracy. The calculation formula is:

$\dot{s}_{v i s c}=\frac{\mu}{T}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) \frac{\partial u_i}{\partial x_j}$                    (9)

where, μ is the aerodynamic viscosity of air, with a value of 1.8 × 10⁻⁵ Pa·s, T is the local temperature, u_i and u_j are the velocity components in different directions, respectively, and the velocity gradient ∂u_i/∂x_j is directly solved from the second moment of the LBM distribution function. Its calculation accuracy is significantly better than that of the traditional Navier–Stokes equation derivation method, and it can accurately capture the viscous dissipation characteristics in near-wall vortex regions. The entropy generation rate of heat conduction under finite temperature difference originates from the non-uniform temperature distribution of the system. The innovation in calculation lies in distinguishing the thermal conductivity of building materials and air, while avoiding temperature singularity problems. The calculation formula is:

$\dot{s}_{\text {cond}}=\frac{k}{T^2}|\nabla T|^2$                                (10)

where, k takes values according to the medium type, with 0.8 W/(m·K) for building materials and 0.026 W/(m·K) for air. The temperature gradient ∇T is calculated from the temperature field solved by FEM using the central difference method, ensuring gradient accuracy. T ≥ 273 K is set, effectively avoiding calculation singularity when T approaches 0, while conforming to the actual temperature range of building thermal environments, thereby achieving accurate quantification of irreversible heat conduction losses and forming seamless connection with the thermal field solution method described above.

The radiation absorption entropy generation rate is the irreversible loss of solar radiation energy conversion into low-grade thermal energy. Its core lies in considering the spectral characteristics of solar radiation, and using band-based integral calculation, breaking through the errors caused by traditional single-band calculations. The calculation formula is:

$\dot{s}_{\mathrm{rad}}=\int_\lambda \frac{1}{T}\left(1-\frac{T}{T_{\mathrm{sun}}}\right) \alpha_\lambda I_\lambda d \lambda$                          (11)

where, T_sun is taken as 5778 K, which is the equivalent blackbody temperature of the sun, α_λ is the spectral absorptivity of the building surface, with a value of 0.6 in the visible light band and 0.8 in the near-infrared band, conforming to the actual spectral absorption characteristics of building materials, and I_λ is the incident radiation intensity directly output by the MC ray tracing method, ensuring the accuracy of radiation parameters. The integration range is set as 0.38–2.5 μm, covering the main bands of solar radiation. The Simpson integration method is adopted, and the integration step is set as 0.01 μm. Numerical verification shows that the integration error can be controlled within 2%. This calculation method accurately matches the band-based solution strategy of the light field described above, achieving refined quantification of irreversible losses of radiation absorption. At the same time, it works collaboratively with the calculation methods of viscous dissipation and heat conduction entropy generation rate, constructing a complete calculation system for each entropy generation rate component, and providing accurate data support for the integration of total entropy generation rate and subsequent sensitivity analysis.

2.3.3 Calculation of total entropy generation rate and evaluation indicators

To achieve comprehensive quantification and comprehensive evaluation of the thermodynamic irreversibility of the wind–light–thermal coupling system, this study innovatively integrates the calculation results of each entropy generation rate component and designs a hierarchical evaluation indicator system, breaking through the limitation of traditional studies that only focus on the total entropy generation and ignore spatial distribution characteristics, and providing multidimensional response quantities for subsequent global sensitivity analysis. The total local entropy generation rate is obtained by directly summing the three components of viscous dissipation, heat conduction, and radiation absorption, accurately reflecting the intensity of irreversible processes in each grid element and achieving spatially refined characterization of irreversible losses. The daily total entropy generation, as the core evaluation indicator of global irreversibility, is obtained by performing a 24-hour spatiotemporal integration of the local entropy generation rate over the entire computational domain. The calculation formula is:

$S_{\text {gen}}^{\text {total}}=\iiint_V \int_0^{24} \dot{s}_{\text {gen}}(r, t) d t d V$                                (12)

Its core lies in using the trapezoidal integration method to solve the time integration and the Gaussian integration method to solve the spatial integration. The trapezoidal integration can accurately capture the dynamic variation characteristics of entropy generation rate over time, while Gaussian integration can improve the calculation accuracy of spatial integration. Numerical verification shows that this combined integration method can control the overall integration error within 3%, which is significantly better than traditional single integration methods. To supplement the evaluation of the thermodynamic uniformity of the indoor–outdoor thermal environment, the coefficient of spatial variation of entropy generation rate is innovatively introduced. The calculation formula is:

$C V=\frac{\sigma}{\mu}$                            (13)

where, σ is the standard deviation of entropy generation rate in the computational domain, and μ is the mean value of entropy generation rate. This indicator can quantify the degree of dispersion of the spatial distribution of entropy generation rate. A smaller CV indicates a more uniform distribution of entropy generation rate and higher thermodynamic efficiency of the indoor–outdoor thermal environment. This evaluation indicator system not only reflects the total characteristics of irreversible losses in the system through daily total entropy generation, but also supplements the spatial distribution characteristics through the coefficient of spatial variation, achieving dual evaluation of the total amount and spatial distribution of irreversibility. It accurately matches the refined output of the multiphysics simulation described above, and provides a comprehensive and reliable evaluation criterion for quantifying the influence of morphological parameters on thermodynamic characteristics in subsequent global sensitivity analysis.

2.4 Global sensitivity analysis method based on polynomial chaos expansion

2.4.1 Construction of polynomial chaos expansion surrogate model

In view of the limitation that multiphysics coupling simulation has high computational cost and cannot directly support large-scale parameter sampling, this study innovatively constructs a PCE surrogate model adapted to the wind–light–thermal coupling system. The core breakthrough lies in accurately matching morphological parameters with thermodynamic evaluation indicators, achieving high-accuracy and high-efficiency response prediction, and providing reliable support for subsequent global sensitivity analysis. The innovative design of model input and output vectors is defined as follows. The input vector is a multidimensional vector composed of five types of core morphological parameters: $x=\left(x_1, x_2, x_3, x_4, x_5\right)$, which respectively correspond to aspect ratio, concave–convex depth, roof inclination angle, facade opening ratio, and eave overhang length, fully covering the controllable morphological parameters defined above. The output vector is: $Y=\left(S_{\text {gen}}^{\text {total}}, C V\right)$, which simultaneously includes daily total entropy generation and entropy generation rate spatial variation coefficient, achieving dual prediction of the total amount and spatial distribution characteristics of thermodynamic irreversibility, breaking through the limitation of traditional surrogate models that only focus on a single output. Based on the uniform distribution characteristics of morphological parameters, Legendre orthogonal polynomials are innovatively selected as PCE basis functions. Their orthogonality can effectively reduce correlation among basis functions and improve model fitting accuracy. The PCE expansion is:

$Y=\sum_{\alpha \in \mathrm{N}^d} c_\alpha \Psi_\alpha(X)$                    (14)

where, $d=5$ is the input parameter dimension, $\alpha$ is a multiindex, $c_\alpha$ is the polynomial coefficient, and $\Psi_\alpha(X)$ is the ddimensional Legendre orthogonal polynomial basis function. The polynomial order is innovatively set to 4 . Numerical tests show that this order can avoid overfitting while fully capturing the nonlinear relationship between morphological parameters and entropy generation rate, balancing model accuracy and computational efficiency, and solving the problem of inaccurate results or low efficiency caused by arbitrary order selection.

The solution and validation of model coefficients rely on an adaptive optimization strategy to ensure high accuracy and reliability of the surrogate model for subsequent global sensitivity analysis. The least angle regression (LARS) algorithm is used to solve the polynomial coefficients $c_\alpha$. This algorithm can adaptively select important polynomial basis functions and parameter interaction terms, automatically remove redundant terms, significantly reduce model complexity, improve computational efficiency, and avoid fitting distortion caused by retaining redundant terms in traditional methods. Using the 64 morphological configurations generated above, 48 are selected as training samples, and the $S_{\text {gen}}^{\text {total}}$ and CV data output from multiphysics simulation are used to complete coefficient estimation. The remaining 16 configurations are used as validation samples. The coefficient of determination R² between predicted and simulated values is used to verify model accuracy, strictly controlling R² > 0.95 and relative error < 5%. This validation standard, verified by thermodynamic analysis, ensures that the surrogate model can accurately reproduce the nonlinear response relationship between morphological parameters and entropy generation rate. Its prediction accuracy fully meets the requirements of global sensitivity analysis. At the same time, the response computation time of a single morphological configuration is reduced from 4–6 hours to millisecond level, greatly improving the efficiency of subsequent global sensitivity analysis, achieving dual optimization of accuracy and efficiency, and laying a solid foundation for Sobol index calculation and response surface construction.

2.4.2 Calculation of sobol global sensitivity index

Based on the high-precision PCE surrogate model constructed above, this study innovatively adopts the MC integration method to calculate Sobol global sensitivity indices. The core breakthrough lies in achieving accurate quantification of the influence of morphological parameters on thermodynamic response variables, while simultaneously identifying independent effects and interaction effects, breaking through the limitation of traditional local sensitivity analysis that cannot capture global correlations and parameter interactions, and providing a reliable basis for key morphological parameter screening. Daily total entropy generation and entropy generation rate spatial variation coefficient are selected as output response variables. Through the combined calculation of first-order and total-order Sobol indices, the influence mechanism of each morphological parameter is comprehensively analyzed. The first-order sensitivity index represents the independent contribution of a single morphological parameter to the variance of the output response:

$S_i=\frac{\operatorname{Var}\left[E\left(Y \mid X_i\right)\right]}{\operatorname{Var}(Y)}$                          (15)

where, $\operatorname{Var}(Y)$ is the total variance of the output response, and $\operatorname{Var}\left[E\left(Y \mid X_i\right)\right]$ is the expectation of conditional variance of the output response after fixing the $i$-th morphological parameter. This index can accurately identify the single parameter that plays a dominant role in thermodynamic irreversibility. The total-order sensitivity index represents the sum of the interaction effects between the $i$-th morphological parameter and all other parameters, as well as its own independent effect:

$S_{T i}=1-\frac{\operatorname{Var}\left[E\left(Y \mid X_{-i}\right)\right]}{\operatorname{Var}(Y)}$                            (16)

where, $X_{-i}$ is the vector of all morphological parameters except the i-th parameter. Its innovation lies in compensating for the inability of the first-order index to capture parameter interaction effects, enabling accurate quantification of the synergistic influence of parameters on entropy generation rate and solving the problem of biased conclusions caused by ignoring parameter interactions in traditional sensitivity analysis.

To ensure the accuracy of sensitivity index calculation, this study performs targeted optimization of computational settings and forms a standardized calculation procedure. The MC sampling size is innovatively set to 10⁶, which is verified to control the calculation error of indices within 2%.

Latin hypercube sampling is adopted, consistent with the uniform distribution characteristics of morphological parameters, effectively avoiding sampling repetition and uneven distribution, and improving sampling efficiency and representativeness. Variance is calculated using an unbiased estimation method, correcting sample variance bias to further reduce computational error and ensure accuracy of Var(Y) and conditional variance calculations. This method is closely linked with the PCE surrogate model described above. By leveraging the efficient predictive capability of the surrogate model, the computational efficiency of global sensitivity analysis is improved by two orders of magnitude. At the same time, through index decomposition and accuracy optimization, fine-grained analysis of morphological parameter influence mechanisms is achieved, providing high-quality sensitivity data support for key parameter identification, response surface construction, and thermodynamically optimal morphology localization.

2.4.3 Morphology–elasticity coefficient and construction of thermodynamic response surface

To accurately quantify the sensitivity of morphological parameter variations on thermodynamic irreversibility, and to break through the limitation that traditional Sobol indices can only represent relative importance of parameters and cannot quantify response sensitivity, this study innovatively defines the morphology–elasticity coefficient to realize a quantitative analysis of the response relationship between morphological parameters and daily total entropy generation. The calculation formula of the morphology–elasticity coefficient is:

$\varepsilon_i=\frac{\partial \ln S_{\text {gen}}^{\text {total }}}{\partial \ln x_i}$                             (17)

Its core lies in using a logarithmic form to construct the relationship between relative changes of parameters and relative changes of entropy generation rate. The physical meaning is clear, and it can directly reflect that when the i-th morphological parameter changes by 1% relatively, the relative percentage change of daily total entropy generation. The sign and absolute value of this coefficient have clear physical meaning. A larger absolute value indicates that the morphological parameter has higher sensitivity on thermodynamic irreversibility. A positive sign indicates that an increase of the parameter will lead to an increase of entropy generation rate, i.e., intensifying system irreversible loss. A negative sign indicates that an increase of the parameter can reduce entropy generation rate and optimize thermodynamic performance. The introduction of the morphology–elasticity coefficient compensates for the limitation of Sobol indices which can only qualitatively judge parameter importance, and realizes quantitative measurement of morphological parameter sensitivity. It provides an accurate basis for prioritizing control of key parameters, and is closely connected with the prediction results of the PCE surrogate model described above, ensuring high accuracy and reliability of coefficient calculation.

The key of the thermodynamic response surface is to focus on key morphological parameters, realize visualization of the relationship between morphology and entropy generation rate, and accurately locate optimal morphology, breaking through the limitation that traditional parameter analysis cannot intuitively present nonlinear relationships and optimal intervals. Based on the results of the Sobol global sensitivity analysis above, the two most important morphological parameters with the largest first-order Sobol index are selected: concave–convex depth and roof inclination angle. Other morphological parameters are fixed at baseline values, i.e., aspect ratio 1:1:1, opening ratio 0.15, and eave overhang length 0.5 m, ensuring the pertinence and simplicity of the analysis. Within the reasonable value ranges of the two key parameters, continuous fine variation is implemented. The concave–convex depth ranges from 0–0.5 m with a step of 0.01 m, and the roof inclination angle ranges from 0–30° with a step of 0.5°. This interval has been verified to be sufficient to capture the subtle influence of parameter changes on entropy generation rate and avoid missing optimal intervals caused by too large step sizes. Through the high-precision PCE surrogate model constructed above, the daily total entropy generation corresponding to each parameter combination is predicted, and the Kriging interpolation method is used to fit the nonlinear relationship between parameters and entropy generation rate. A three-dimensional response surface and two-dimensional contour map of concave–convex depth–roof inclination angle–daily total entropy generation are plotted. By analyzing the surface shape and contour distribution, the local minimum region of daily total entropy generation is accurately identified. This region is the thermodynamically optimal building morphology interval, achieving quantitative and visual identification of optimal morphology. This provides a clear quantitative target for subsequent experimental verification and low-entropy building morphology design, and establishes an intuitive relationship between morphological parameters and thermodynamic performance, completing the full system of global sensitivity analysis.

This paper conducts a systematic study on the sensitivity of indoor–outdoor wind–light–heat coupled processes with respect to building form variation. Taking entropy production rate as a unified thermodynamic evaluation index, a complete methodological system is constructed, namely “form parameterization – multi-physics weak coupling simulation LBM-MC-FEM – entropy production decomposition evaluation – PCE-based global sensitivity analysis”. This system effectively fills the application gap of non-equilibrium thermodynamics in the field of building form optimization and provides a standardized technical pathway for accurate thermodynamic optimization of building forms. Experimental studies confirm that the newly proposed shape non-uniformity index ξ shows significantly higher correlation with total entropy production rate than traditional shape parameters, and can accurately quantify the intrinsic relationship between geometric complexity and thermodynamic irreversibility of buildings; the self-developed LBM-MC-FEM weak coupling framework satisfies engineering and academic accuracy requirements, with all errors less than 8%, and computational efficiency significantly higher than both strong coupling methods and commercial software integration approaches, greatly improving the engineering applicability of multi-physics simulation; concave–convex depth is consistently the most critical geometric parameter affecting entropy production in the wind–light–heat coupling process, with its Sobol first-order index ranging from 0.40 to 0.48. Meteorological conditions significantly affect the sensitivity of each shape parameter, with influence weights of wind speed, solar radiation, and temperature being 42.3%, 37.6%, and 20.1%, respectively; there is a significant interaction effect between concave–convex depth and roof inclination angle, and their interaction contributes 12% of the variance of total entropy production rate. Based on this, the thermodynamically optimal building form range is identified as concave–convex depth 0.18–0.22 m and roof inclination angle 13–17°. Within this range, the thermodynamic irreversibility of the system under a single-day condition is minimized, and the spatial variation coefficient of entropy production rate decreases from 0.24 in the baseline configuration to 0.18, indicating optimal indoor–outdoor thermal uniformity. The overall accuracy of the proposed full-process method is high, with relative errors between simulation results and experimental measurements all below 8%, which can effectively guide low-entropy building form design and achieve coordinated improvement of building thermodynamic efficiency and energy consumption optimization.

The method proposed in this paper solves the core problems in traditional building form studies, including inconsistent evaluation indicators, coarse coupling simulation, and sensitivity analysis limited to local effects. It enriches the application scenarios of non-equilibrium thermodynamics in the field of building physics and provides important theoretical support and technical reference for building energy consumption optimization under global carbon neutrality goals. Annual hourly simulation results show that the representative optimal-range configuration reduces total annual entropy production by 24.4% compared with the baseline configuration, and by 42.5% compared with the maximum shape non-uniformity configuration; the annual average thermal environment variation coefficient is 0.19, and building energy consumption is reduced by 18.7% compared with the baseline configuration and by 32.1% compared with the maximum shape non-uniformity configuration. Different geometric configurations can effectively reduce thermodynamic irreversibility in multi-physics coupling processes, and the optimization effect is stable and reliable. At the same time, this study does not consider the influence of thermodynamic properties of building materials and dynamic shape variation. Future work will further couple material thermodynamic parameters and construct a dynamic form optimization model to expand the engineering application scope of the method and promote the implementation of low-entropy building design concepts in engineering practice.