Image-based microstructural simulation of thermal conductivity for highly porous wood fiber insulation boards

The thermal conductivity of wood fiber insulation boards is significantly influenced by the microstructure of the fiber network and in general, the efficiency of wood fiber insulation boards increases with porosity. For higher raw densities, the raw density is a good predictor for the thermal conductivity. For lower raw densities however, this simple relation does not hold anymore. Here, structural information gained from 3D computed tomography images at several scales, modeling of the microstructure, and numerical simulation of the thermal conductivity are combined to get deeper insight into which and how microstructural features influence the thermal conductivity. The model-based simulation as described here shows that the presence and orientation of wood fiber clusters impact the thermal conductivity significantly.


Introduction
The microstructure of highly porous insulation mats made of wood fibers is investigated in order to better understand the structure-property relationship for improving the thermal insulation properties. The final goal is to reduce heat transfer through the wood fiber mat. The European Commission (2020) states that buildings cause 40% of the final energy consumption and 36% of the greenhouse gas emissions in the European Union. 75% of the building stock lack energy efficiency. Thermal insulation is thus crucial to reduce energy consumption of new buildings as well as for energy retrofitting of the building stock.
Here, we report on thermal conductivity measurements, 3D imaging and subsequent processing and analysis of the image data, and on thermal conductivity simulations in virtual microstructures aiming at the microstructural optimization of wood-based and thus renewable insulation products. Micro-computed tomography ( μCT) using synchrotron radiation (SR) has been frequently used to spatially image and morphologically characterize the microstructure of cellulose fiber materials (Faessel et al. 2005;Lux et al. 2006;Peyrega et al. 2011). Lux et al. (2006) derive for very light wood-based thermal insulation material the size distributions of fibers and pores as well as the local orientation of the fibers using classical morphological tools. Thömen et al. (2008) use microstructural characteristics from SRμ CT images obtained by Walther (2006), Walther et al. (2006) and Walther and Thömen (2009) for simulation models for permeability and thermal transport.
Geometric models for cellulose materials are commonly based on sedimentation as described by Niskanen and Alava (1994) and Provatas et al. (2000). Peyrega et al. (2009) incorporate fiber-fiber interaction in 2D and elastic deformation of the fibers in 3D. The virtual compression suggested by Provatas and Uesaka (2003) yields realistically rough surfaces. Building on Niskanen and Alava (1994), Wang and Shaler (1998) and Wang (2000) develop a 3D geometry model for medium density fiber boards consisting of hollow, flattened fibers that bend depending on the structure they fall on, see also Walther (2006). All the sedimentation-based models are very complex. Thus, model fitting in the sense of deriving the model parameters from features measured based on images of the real structure, as described for example by Redenbach et al. (2014), seems out of reach.
Boolean models in the sense of Stoyan et al. (1995) or Schneider and Weil (2008) are much more accessible for fitting due to Miles' formulae relating the model parameters to microstructural features observable in 3D images (Ohser and Schladitz 2009, Sect 7.4.). Grains drawn from the distribution of the so-called typical grain are placed independently, uniformly randomly in the considered volume. The cellulose fibers thus do not interact in these models. Faessel et al. (2005) use Boolean models of curved fibers for modeling cellulose fiber materials with volume fractions between 25 and 75% and subsequently simulate thermal conductivity. Carvajal et al. (2019) give a comprehensive overview of modeling heat transport in fibrous insulating materials and compare especially glass fiber materials with relatively high densities. Predictions of several models for computing the effective thermal conductivity are compared to certified measured values for a standard glass fiber board. Carvajal et al. (2019) take into account coupled effects of radiation and conduction heat transfer although the paper focuses on the conductive thermal transport, only. The effective thermal conductivity can be computed precisely by numerical homogenization, where three periodic boundary value problems are solved in a representative volume element (RVE) to obtain the full second order thermal conductivity tensor as described for example by Monchiet and Bonnet (2013) and Höller et al. (2020). For the solution of the steady-state heat conduction equation there are several discretization methods-finite elements (FEM) as used for example by Höller et al. (2020), finite differences (FDM) as applied for example by Linden et al. (2015) and Wiegmann and Zemitis (2006), and many fast numerical solution methods, that can be used in the context of numerical homogenization. Methods based on the fast Fourier transform (FFT) on regular grids like those of Monchiet and Bonnet (2013) and Merkert et al. (2015) and methods that employ local mesh refinements like the one of Linden et al. (2015) are widely used in order to reduce the numerical effort.
Several software systems generate deterministic or stochastic structures and simulate macroscopic properties in realizations of these structure models. PAKKA by Niskanen et al. (1997) uses the sedimentation algorithm of Niskanen and Alava (1994) to generate systems of fibers with rectangular cross-section, isotropically oriented within the xy-plane, as described by Hellén et al. (2002). PaperSim by Thompson (1992) builds on Heyden (2000) and models curvatures and orientations, but no interaction. The fiber networks are usually fit to the real material based on a couple of simple features like the area weight. For the numerical examples in this work, the software GeoDict ® (2022) is used. This software provides two algorithms for the computation of the effective thermal conductivity-the explicit jump immersed interface method (EJ-HEAT) by Wiegmann and Zemitis (2006) and the LIR method by Linden et al. (2015).
In this paper, we study four wood fiber insulation boards representative for two different dry production processes, both with differing raw density ranges. The materials as well as methods applied to measure their thermal conductivities are described in "Materials and methods" section. For the two highly porous variations, 3D imaging and geometric characterization of their microstructures are detailed in "Computed tomography" and "Image processing and analysis" sections. "Microstructure modeling" and "Simulating thermal conductivity" sections yield methodology for model generation and thermal conductivity simulation. Results are presented in "Results" section and summarized in "Conclusion" section.

Wood-based high porosity insulation mats
Four commercially produced types of wood fiber insulation boards, abbreviated with B3, B4, B5, and B6, are investigated. The rigid boards B3 and B4 were produced from the same fiber assortment and bonded with isocyanate in a dry process. They differ however in raw density, with (B3) ≈ 170 kg/m 3 , whereas (B4) ≈ 105 kg/m 3 . Products B5 and B6 are flexible mats, bonded with a small fraction of thermoplastic binding fibers. The raw density of B5 is lower than the one of B6-(B5) ≈ 50 kg/ m 3 and (B6) ≈ 60 kg/m 3 . Moreover, B5 features coarser fibers than B6 as the geodesic lengths x LG and thicknesses x E listed in Table 1 show. More precisely, Table 1 contains the projection area weighted percentiles of geodesic length x LG and thickness x E , measured using an optical scanner with a resolution of 1200 dpi according to ISO 9276-6:2008ISO 9276-6: (2008. The percentiles are recalculated by IST AG's Fibre-Shape, version 5.2.2. Only fibers exceeding 100 μ m for both length and thickness were taken into account.
From these four products, specimens for thermal conductivity measurements and for CT scans were taken from six boards, abbreviated with B3-01, B4-01, B5-01, B5-85, B6-01, and B6-83. Note that B5-01 and B5-85 as well as B6-01and B6-83 are boards of the same product lines, respectively, but produced on different days. Microstructural differences are thus expected primarily due to the natural variability of the wood processed.
For thermal conductivity measurements, specimens of size 500 mm × 500 mm × board thickness were used. For CT imaging, smaller samples were cut from these plates. Additionally, unbound material-fibers and wood chips and chunks consisting of one or several layers of fibers, respectively-was provided for the SRμ CT experiments.

Measuring thermal conductivity
Thermal conductivity was measured using a guarded hot plate apparatus according to EN 12667 (2001). The relative measurement uncertainty is expected to be less than 1%. For each type of insulation, two specimens with dimensions 500 mm × 500 mm × thickness were used. The actual measurement area of the guarded hot plate is a square of side length 400 mm for B5-85and B6-83, and 300 mm in all other cases.
All measurements were performed on a stacked pair of specimens with the heating plate in between and cold plates on the outer sides. That way, the mean of the thermal conductivity for upward and downward heat flux directions is measured. This setup is favorable in terms of measurement uncertainty and representativeness, because any convection phenomena that can be influenced by gravity are excluded, even though convection is unlikely to occur in the observed case due to small temperature differences and high flow resistance. For B5-85 and B6-83, only one specimen was available. Their original thickness of 120 mm was therefore sliced into two slabs of thickness of 50 mm.
All measurements were performed at a temperature difference of about 15 K. The mean temperature m was varied from about 10 • C to 50 • C for samples B3-B6 and in the wider range from − 21 to 70 • C for B5-85 and B6-83. According to Kast (2008), the thermal conductivity can be predicted by a third-order polynomial ( m ) = a 3 3 m + a 2 2 m + a 1 m + a 0 for sufficiently small variation of m and temperatures below 100 • C. Before measurement, all materials were dried at 70 • C in a ventilated oven with pre-dried air supply until mass equilibrium is reached. Directly after drying, the specimens were sealed by a PE membrane and stored in a desiccator until the measurement is taken. To avoid re-adsorption of moisture, the specimens remain packaged during measurement. To exclude the additional thermal resistance of the membrane, the thermo-couples were led through the PE membrane and attached directly to the specimen's surface. Right after measurement, the specimens were weighed and re-dried, to determine the moisture content after the measurement. Very low moisture contents of 0.1-0.6 mass% with mean 0.3 mass% indicate a good sealing.
To interpret the measured thermal conductivities, the relevant heat transfer phenomena have to be identified. An influence of convection can be excluded according to ISO 10456:2007+ Cor. 1:2009: The most unfavorable material properties are represented by specimen B5-85, see Table 2. Thus, the flow resistance is measured for this specimen ( r = 4718 Pa s/m 2 ) resulting in the modified Rayleigh number Ra m (B5-85) = 0.52. This value is far below the critical values 15-30 for upward heat flux and 2.5 for horizontal heat flux according to ISO 10456:2007+ Cor. 1:2009. Treml et al. (2019) determined the fractions of heat flux due to conduction in the solid phase, conduction in the gaseous phase, and radiation experimentally for Table 2 Measured thermal conductivity , raw density , thickness t, temperature difference Δ , mean temperature m and coefficients of the polynomial from "Measuring thermal conductivity" section 1 3 B3 and a second material. In both materials, radiative heat transfer was found to contribute a very similar portion to the total heat flux-13% for B3 and 14% for the other material.

Computed tomography
3D images of the wood fiber insulation board's microstructure at lateral resolutions from sub-μ to about 25 μ m were taken to shed light on how specific microstructural features are related to thermal conductivity.

Sample preparation
For the CT imaging, cuboidal samples with a short edge length of 5 mm have to be prepared. This is extremely demanding due to porosities of 96% and higher of B5-xx and B6-xx and the fact that mechanical strength is achieved only by a low proportion of thermoplastic binding fibers.
The fibers have to be fixed to cut the specimen in such small dimensions without influencing the microstructure. Infiltration experiments (also under vacuum conditions) with typical resin systems used for fixation of microscopic specimens resulted in residual pores dominating the X-ray absorption contrast. Therefore, a dedicated sample preparation method was developed.
To avoid the problem of the low X-ray absorption contrast between the fibers and the fixation, the fixation has to be removed after cutting the specimen. Therefore specimens with larger dimensions were frozen in water. They were cooled in liquid nitrogen before being immersed in water to avoid swelling of the fibers. This procedure was repeated several times until a solid block of ice is formed. This block can then be cut on a band saw to the required dimensions of 5 mm short edge length. Directly after cutting, the frozen specimen were put into a plastic container to stabilize the fiber structure after the removal of the water. Finally, the ice-fixation was removed by sublimation in a vacuum chamber to avoid distortion of the microstructure by liquid water during melting.
Three specimens of each variation (B3-B6) were prepared this way for SRμ CT measurements. Additionally, four column-shaped samples of B5-85 and B6-83 of short edge length 15 mm and covering the complete board thickness (about 100 mm) were prepared the same way for laboratory CT imaging.

Sub-micro computed tomography using synchrotron radiation
The samples were imaged at ID19 of the European Synchrotron Facility (ESRF) in Grenoble, France. A photon energy of 19 keV (pink) was utilized with a singleharmonic undulator (u17.6 type) as a source. The resulting homogeneous wave front is excellently suited for propagation-based phase contrast imaging (camera-detector distance: 20 mm, detector pixel size: 0.65 μm). Paganin et al.'s (2002) singledistance phase-retrieval approach combined with an unsharp mask was applied to facilitate the segmentation of the solid material.

3
Wood Science and Technology (2023) 57:5-31 Micro computed tomography using a laboratory device ITWM's CT device featuring a Feinfocus FXE 225.51 X-ray tube with maximum acceleration voltage 225 kV, maximum power 20 W, and a Perkin Elmer flat bed detector XRD 1621 with 2048 × 2048 pixels was used. The tube voltage of 120 kV is rather low and the integration time of 1 s high to accommodate the low absorbing cellulose. Tomographic reconstructions were obtained from 1200 projections, each averaged over four, resulting in an overall integration time of 4 s. The voxel edge length of 12.8 μm is a trade-off between the wood fibers' wall thickness and the field of view that can be captured by one scan. Visualizations are shown in Fig. 1. The four column-shaped samples were scanned at three height steps (bottom, center, top) resulting in over all 12 scans per board examined.

Extraction of typical objects
The gray value (16bit) images were denoised by a 5 × 5 × 5 median filter. Binarization of the SRμ CT images is then easily achieved by a global gray value threshold.
Separation of objects is, however, not straightforward. The solid material forms essentially one connected component. Morphological transformations are not able to separate the objects as structuring elements have to be large enough to span over the large lumen areas. Using the lumen as a marker as suggested by Malmberg et al. (2011) does not work neither as the lumen is not completely separated from the pore space in many cases-due to very thin walls not being completely captured in the 3D image as well as due to fiber walls being ripped during production. Fiber separation approaches as the one of Viguié et al. (2013) do not work for the chips and chunks.
Here, a tailor suited method based on locally measured densities of the intrinsic volumes as described by Ohser and Schladitz (2009) was therefore used. More precisely, the specific surface area in cubic sub-volumes of an edge length slightly larger than the fiber thickness was considered. This value is higher in chunks, where the sub-volumes intersect with more than one layer of wood fibers compared to one layer thick chips and fibers. Analogously, but starting with smaller sub-volumes, chips and fibers were separated in a second step. See Fig. 2 for an example. The separation works fine as long as fibers and chips are not touching in a large region.

Segmentation of the solid component based on laboratory CT image data
Binarization of the laboratory CT images is much more difficult. In this case, nominal resolution is reduced significantly compared to the SRμ CT measurements in order to capture samples of a representative size. As a consequence, the partial volume effect affects the binarization strongly. Choosing a global gray value threshold such that the expected volume fractions are retained, results in a very rough, friable structure. Thus, we deliberately overestimated the solid volume fraction. That means, for the sake of preserving local connectivity, in particular of the fine individual fibers, a too high solid volume fraction is accepted, see Fig. 3 in Section Image-based structure characterization. In the absence of a ground truth, there is no Fig. 3 Solid volume fraction V V of samples B6-83 and B5-85, derived from the binarized CT images in sub-volumes of size 400 3 voxels corresponding to (5.12 mm) 3 . Sub-plots a-d correspond to the four column-shaped subsamples as described in "Computed tomography" section 1 3 means to validate the segmentation result. However, scrupulously treating all data sets in the same way ensures comparability of results. Note, moreover, that the error in estimating geometric characteristics other than the solid volume fraction is lower, see Table 3. Thus the oversegmentation seems acceptable.
Pre-processing shall ensure comparability of the gray value ranges of all scans: First, gray values exceeding the 99.99% percentile of the gray value distribution were set exactly to this value to remove bright imaging artifacts. Subsequently, gray values are locally transformed to remove gray value variations due to locally strongly varying material density. More precisely, in cubes of roughly 160 voxels edge length, the gray values are centered at mean zero and the gray value variance is fixed to one. As a result, voxels representing material, but also noise voxels clearly stick out.
Then hysteresis thresholding yields a binary mask for the solid phase: A high gray value threshold selects marker voxels surely belonging to the solid phase. A lower threshold yields all potential material voxel candidates, including bright noise voxels. This candidate set is slightly enlarged by a dilation with a 3 × 3 × 3 voxel cube as structuring element. The binary mask is then achieved by growing the marker set, restricted by the candidate set. Non-connected voxels are removed as noise.
The locally normalized gray value image is median filtered. Noise is removed by applying the just derived mask. A global threshold chosen by Otsu's method (Otsu 1979) finally yields the desired binarization of the solid component.

Geometric characterization of the solid component
Based on the binarized voxel images, volume V, surface area S, and related characteristics were deduced by Ohser's algorithm (Ohser et al. 1998). This algorithm yields generalized projection areas in 13 directions given by the image lattice as a by-product of the surface area measurement as described in Ohser and Schladitz (2009, Chapter 5). Here, the generalized projection areas in x, y, and z-directions were used as a measure of the degree to which the solid structure is oriented within the xy-plane. More precisely, the anisotropy factor a was calculated as the ratio of the projection areas in z-direction and the mean of those in x-and y-direction: For a perfectly isotropic structure, we would get a = 1 . If the wood fiber and fiber bundle structure is oriented transversally isotropic with respect to the z-direction and preference for the xy-plane, then p 2,x ≈ p 2,y and p 2,z > p 2,x and thus a > 1.
Surface area and volume divided by the total volume of the observation window yield the specific surface area S V and the solid volume fraction V V . Additionally, we use the surface area per volume of the solid component S/V, which is called specific surface area in biology and medicine when bone is investigated (Bouxsein et al. 2010). In general, we have S∕V = 1∕(4L) , where L denotes the mean intercept length (MIL), see Ohser and Schladitz (2009, Remark 5.6). In the special case of (1) a ∶= 2p 2,z ∕(p 2,x + p 2,y ).
the structure under consideration being a spatially homogeneous system of endless, non-overlapping tubes (hollow cylinders) of circular or rectangular cross-section with constant wall thickness t, then S∕V = 2∕t: Let Λ denote the length density of the cylinder system in the sense of Schneider and Weil (2008). That is, Λ is the expected total cylinder length per unit volume. Denote by R > 0 the outer and by r > 0 the inner radius of the circular cross-section, respectively. We have then as the wall thickness is t = R − r . In the case of a hollow square cross-section with outer edge length L > 0 and inner > 0 , the analogous consideration yields as the wall thickness is t = 1 2 (L − ).

Microstructure modeling
To derive a geometry model that reflects correctly essential features of the real structures, first, objects-fibers, chips, chunks-are extracted from the high-resolution SRμ CT images as described in "Image processing and analysis" section. These objects are modeled individually by hollow spherical cylinders (fibers) or arrangements of them using GeoDict ® (2022). Chips are formed by one layer of fibers, chunks by several layers, see Fig. 4 for the example of a chip. Diameters, volume, arrangement, and curvature of the fibers are chosen to resemble the real counterparts. Rescaling and reorienting these building blocks, finally placing them uniformly randomly, yields the desired model structures. More precisely, the model structures are realizations of a Boolean model with the typical grain being distributed according to a finite mixture of shapes. In "Simulated thermal conductivities" section, this strategy is used to generate virtual geometries resembling B6-83, starting from the solid volume fraction and altogether 27 grain shapes-7 fibers, 16 chips, and 4 chunks. Table 4 yields representative results from an extensive simulation study including the resulting anisotropies and thermal conductivities.
Clearly, model structures M17 and M18 fit the real structure best with respect to anisotropy. Model structures M10-M12 deviate and rather feature a higher proportion of surface normals in x-direction. The renderings in Fig. 5 favor model structure M7.

Simulating thermal conductivity
Heat transport in the porous structure is composed of thermal conduction, thermal radiation in the pore space, and thermal convection in the pore space. In this section, only the thermal conduction fraction is computed by numerical homogenization. As discussed in "Measuring thermal conductivity" section, convection can be neglected due to small temperature differences and high flow resistance. The radiative contribution has, however, to be taken into account for good agreement of simulated and measured thermal conductivities. Often, in simulations, the radiative component is approximated by an increased temperature-dependent thermal conductivity of air. This approach is not chosen here as we aim at temperature-independent effective quantities. Moreover, varying the microstructure as in "Simulated thermal conductivities" section would have no significant impact. We therefore added a radiative contribution estimated from experiments, see "Measuring thermal conductivity" section. Note that approximating the effective thermal conductivity following Treml et al. For simulating the thermal conductivity at the micro-scale, we solve the stationary heat equation using GeoDict ® (2022). The boundary conditions are set periodically in the lateral directions x and y and a temperature difference is applied in the through direction (z, board thickness). Then, GeoDict ® computes the temperature distribution for all computational cells (voxels). Employing standard homogenization theory, we finally obtain the effective thermal conductivity of the cellulose structures. The bulk thermal conductivites of air and the cellulose constituents as well as the settings for the GeoDict ® module ConductoDict are listed in Table 5. For settings not listed, the defaults are kept.
The ConductoDict module of GeoDict ® provides two solvers with different discretizations and solution methods. One method is FFT based and uses a finite difference method for the discretization, which captures coefficient jumps via an explicit jump implicit interface method as described by Wiegmann and Zemitis (2006). The second method is a standard finite volume method on a LIR tree space partition as mesh as devised by Linden et al. (2015). Both methods can be used to solve the problem considered here. The finite volume method was chosen because the contrast of the thermal conductivities is 380∕25.7 = 14.8 , i.e., not very high. Adaptive mesh refinement was applied to increase accuracy, and the iterative solution was accelerated using a multigrid scheme. This solution variant can be easily selected in ConductoDict via corresponding options. evacuated = radiation + solid = k rad T 3 + solid yields k radiation = 2.33 × 10 −7 mW∕(mK 4 ) and solid = 1.89 mW∕(mK). rad = 2.33 × 10 −7 ⋅ 293 3 mW∕(mK) = 5.86 mW∕(mK)

Experimental
The measured thermal conductivities are reported in Table 2. Thermal conductivity measured at a mean temperature of 10 • C ( 10 • C) as a function of the raw density of the material is plotted in Fig. 8. The observed difference of 10 • C between B3-01 and B4-01 ( Δ 10 • ≈ 4.7 mW/(mK)), is well explainable by the increasing proportion of heat transfer due to conduction in the solid structure caused by increasing raw density ( Δ ≈ 65 kg/m 3 ).
In comparison, Δ 10 • ≈ 3.9 mW/(mK) between the lowest (B6-01) and the highest value (B5-85) in the group of flexible mats is remarkably high, compared to the small difference in raw density Δ ≈ 9 kg/m 3 . This indicates that the thermal conductivity of the group of flexible mats is more strongly influenced by the microstructure of the fiber network.
For a sufficiently small variation of m , the thermal conductivity can be predicted by linear regression. For wider variation of m , the influence of temperature on the radiative heat transfer according to Stefan-Boltzmann's law is visible, causing a faster than linear growth of for increasing m , see Fig. 8. In this case, a third-order polynomial should be fit, see Kast (2008) in Section Measuring thermal conductivity. Measurements at further temperatures could increase the accuracy of the regression curves and are, however, beyond the scope of this work.
In general, the higher the raw density of the specimen, the lower is the temperature dependency of thermal conductivity. A closer look at the regression curves and the thermal conductivity gradients in Fig. 8 reveals deviation from this general rule of thumb. This indicates that the proportion of radiative heat transport increases with decreasing raw density. For the group of flexible mats, this effect is relevant to explain the increase in thermal conductivity with decreasing raw density.

Image-based structure characterization
Based on the binarized μ CT images, the densities of the intrinsic volumes and further characteristics derived from them can be efficiently measured by Ohser's algorithm. However, as explained in "Image processing and analysis" section, the μ CT images cannot be binarized such that the solid volume fraction is properly reflected and the SRμ CT scans do not cover representative volumes. Consequently, we cannot expect to derive realistic values of geometric characteristics like solid volume     . 4 Volume renderings of a chip (left) extracted as described in "Image processing and analysis" section and its counterpart generated in GeoDict ® combining hollow spherical cylinders fraction or specific surface area based on these images. The solid volume fraction derived from these images will clearly exceed the true value. Nevertheless, geometric characteristics measured based on the binarized μ CT images can be used to compare the microstructures and to reveal possible differences as all samples have very similar chemical compositions, samples are prepared and scanned exactly the same way, and finally the images are scrupulously treated identically, too. The analysis focuses on the counter-intuitive pair of samples B5-85 and B6-83, where B6-83 yields the lower thermal conductivity although its raw density is higher than the one of B5-85, see Table 2. As described above, four subsamples are taken from each of the two boards. These column like samples are addressed as B5-85 a, b, c, d and B6-83 a, b, c, d, respectively. Three scans of top, center, and bottom arise from each of them as described in "Computed tomography" section.  Table 4 1 3 The microstructure varies strongly. In order to decide, whether measured geometric features of B5-85 and B6-83 differ significantly, they are derived from subvolumes of size 400 3 voxels yielding 16 data points per subsample and height. These are shown as boxplots in Figs. 3, 6, and 7.
The anisotropies of the two samples quantified by a as defined above (1) are in the same range, but a( B6-83) is significantly higher than a( B5-85) , see Figs. 6 and 9. A Kolmogorov-Smirnov test confirms this visual impression given by the box and cumulative distribution function plots.

Surface area per solid volume
Analysis of the binarized images shows not only the difference in anisotropy but also a surplus of surface voxels in sample B6-83. This could indicate finer fibers but also be due to just slightly increased surface roughness. These two effects cannot be separated as a morphological closure even with a very small structuring element alters the structure considerably. In particular, a significant fraction of the lumen is filled afterward. A subsequent opening for morphological regularization would remove too much of the wall structure and thus whole fibers. Figure 7 yields S/V. The observed values of 50-75 mm −1 are within the rather wide range that can be expected 2.5-2000 mm −1 . However, the real structure surely does not fulfill the model assumptions from "Image processing and analysis" section, in particular w.r.t. surface roughness. The true surface is surely rougher, yielding a larger surface area. On the other hand, the preprocessing and binarization imply considerable smoothing and increased values for V V . Thus, again, the values reported here surely differentiate the two boards well. They are, however, no estimates of 2/t in the strict sense.

Representativeness of the analyzed sub-volumes
The sub-volume size is bounded from below by the pore size. On the other hand, we need enough sub-volumes to investigate the variation between samples B5-85 and B6-83 and within each mat, that is between locations a-d. Representativeness of the chosen 400 3 voxel subvolumes was checked following the idea of Kanit et al. (2003), Dirrenberger et al. (2014) as adapted to the 95% confidence interval by Easwaran (2017). That is, 144 sub-volumes of size 400 3 were analyzed, 24 of size 700 3 , and 12 of size 800 × 800 × 1500 voxels, and a power function was fit to the point variance as a function of sub-volume size. This finally relates the relative estimation error Fig. 6 Anisotropy of samples B6-83 and B5-85, measured by a from (1) for sub-volumes of size 400 3 voxels corresponding to (5.12 mm) 3 . Sub-plots a-d correspond to the four column-shaped subsamples as described in "Computed tomography" section. The degree of anisotropy varies with subsample position. The difference in anisotropy, varies, too. However, over all, B6-83 is clearly more anisotropic than B5-85 compared to the mean of the respective characteristic with the size of the representative volume.
See Table 3 for the RVE sizes when limiting the relative errors to 2% or 2.5%. For estimating the solid volume fraction of B5-85, the 400 3 sub-volumes are not big enough to ensure a relative error below 2% at the 95% confidence level. However, for a relative error below 2.5% they are.
Having estimated the relative error at the 95% confidence level, we can check whether the differences of the samples B5-85 and B6-83 are significant. Fig. 7 Fraction of surface voxels of samples B6-83 and B5-85, measured by S/V for sub-volumes of size 400 3 voxels corresponding to (5.12 mm) 3 . Sub-plots a-d correspond to the four column-shaped subsamples as described in "Computed tomography" section 1 3 Wood Science and Technology (2023) 57:5-31

Simulated thermal conductivities
A striking experimental result reported in Table 2 is that specimen B5-85 has a higher thermal conductivity than B6-83, even though the density is lower. In this section, simulations will be used to show that such an effect cannot be ruled out.
In order to shed more light on the difference in thermal conductivity observed for samples B5-85 and B6-83, the sensitivity of the thermal conductivity was examined with respect to geometric microstructural variations, based on geometric models starting close to sample B6-83. The starting point is the model structure found to correspond most closely to the real structure as observed in the μ CT images in "Microstructure modeling" section. Figure 5 shows a 160 3 voxel sub-volume of the insulation material with a voxel length of 12.8 μ m. The solid material comprises 3.68% of the volume including 0.7% chips and 0.3% chunks, very strongly oriented in the xy-plane. The remaining 2.68% of the volume being solid consists of wood fibers with Gaussian distributed diameters with mean 50 μ m and standard deviation 20 μ m. The fibers are 1 mm long and strongly oriented in the xy-plane, too.
For the sensitivity study, the following model parameters are varied: solid volume fraction, fraction of chips and chunks, orientation, diameter and length of the wood fibers. Note that the chips and chunks are built of fibers. Thus, orientation of these larger objects is just defined by the orientation of their building blocks. The wall thickness is fixed to 0.15 of the outer diameter. That is, the thicker the fibers, the thicker also their walls. All estimations are based on five simulated structures.

Fig. 9
Empirical cumulative distribution functions of the three geometric characteristics together with the 95% confidence intervals. B5-85 and B6-83 do differ significantly, not only with respect to solid volume fraction V V but also in structure fineness S/V and anisotropy a Table 4 reports the simulated thermal conductivities for the virtual structures generated according to "Microstructure modeling" section together with selected geometric characteristics. In contrast to the measured ones, the simulated effective thermal conductivities in Table 4 do not comprise thermal radiation as the focus here is not on comparing measurements and simulations. Here, we rather concentrate on the sensitivity of the geometric microstructure features in order to find ways to optimize the insulation material.
The features of B6-83c and the cut-out from it reported in Table 4 are estimated based on the CT image binarized as described in "Image processing and analysis" section. The solid volume fraction V V is therefore higher than specified as explained in "Image-based structure characterization". The structure fineness S/V is essentially in the range of the reciprocal voxel edge length as nearly every voxel is a surface voxel. Thus, S/V found in the CT images is lower than the respective value in the model realizations. Figure 10 yields the simulated thermal conductivities depending on the degree of variation of the respective parameters. It thus provides an overview over the sensitivity and fluctuation of the thermal conductivity. Keeping in mind that the variation in solid volume fraction investigated is very small compared to the feasible range while the anisotropy parameter runs through a considerable part of the reasonable range, and we see that as expected, variation in solid volume fraction (center left) has by far the strongest impact on the thermal conductivity. Note that the chips and chunks dominate this effect. This is apparent by the progressive curve shape (center right). Furthermore, fibers pointing in the main direction of the temperature gradient (bottom) increase the heat conduction considerably. Variations of the size of the individual wood fibers (top left and right) cause the least changes.
In the simulations, structure M16 has an approximately 54% higher solid volume fraction V V than structures M10 and M11. However, structure M16 has a significantly smaller thermal conductivity than M10 and a slightly smaller one than M11. Thus, the sensitivity analysis supports the experimental results regarding samples B85-5 and B86-3: Unfavorable orientation of wood fibers and chunks can overturn the effect of a lower solid volume fraction.
Visualizations of the temperature and the heat flow in Fig. 11 help to understand that fiber bundles and fiber chunks, which remain aligned in the direction of the heat flow (z-direction) during production (pressing), transport a lot of heat, and practically act as thermal bridges. Temperatures and heat fluxes in the pores are not shown in Fig. 11.
The simulated temperature fields for structure variants M12, M17, and M20 visualized in Fig. 11a, c, and e are similar. There are no fibers, chips or chunks in the structures that conduct heat so well that the temperature at any point in the central layer is as high or as low as those at the top and bottom boundary, respectively.
To understand why the thermal conductivities of the structures differ, we look at the local heat fluxes: Structure M17 has the smallest simulated thermal conductivity apart from variants M13 and M14, which have smaller solid volume fractions. The corresponding heat flux in Fig. 11d is much more homogeneous than the ones in M12 and M20. This is caused by the lower volume fraction of chips and chunks and the higher anisotropy in M17. Structure M12 is less anisotropic than M17. In 1 3 Fig. 10 Sensitivity of thermal conductivity with respect to variation of selected geometric features of the model structures. Fiber thickness, fiber length, and fraction of chunks as prescribed. Anisotropy factor a and solid volume fraction V V as measured in the images of the realizations Fig. 11a, more fibers and chunks can be seen, which are more vertically oriented and have a high heat flux mainly because of the vertical orientation. Structure M20 has a higher percentage of chips and chunks compared to the other model structures. Figure 11e shows that the heat flux in these larger objects is high when they are more vertically oriented.

Conclusion
In this paper, dedicated sample preparation, SR-sub-μ and μCT, and 3D image analysis were used to analyze the intricate microstructure of highly porous wood fiberbased boards for thermal insulation. It was proven that computed tomography on several scales yields microstructural information enabling detailed geometric characterization. The relation of microstructural features and measured thermal conductivities was established via quantitative image analysis and model-based numerical simulations.
The model-based simulation of the thermal conductivity allowed to study sensitivity with respect to microstructural features. Next to the solid volume fraction, the orientation of surface normals turned out to be most crucial. In the future, production parameters should therefore be optimized such that structural components are as strongly as possible oriented within the board plane while preserving the mechanical stability of the resulting mat. Rather low sensitivity was found with respect to fiber lengths or thickness. Hence, recycled fibers could be used, too. Finally, avoiding wood chunks and chips by further softening and refining them to separate fibers is desirable, however significantly less important than the orientation. As a consequence, the choice of the wood type is important only as far as it influences solid volume fraction and orientation.
Several challenges have to be solved before simulated and measured thermal conductivities can be compared quantitatively, too. New approaches for extracting walls thinner than one voxel from CT images are needed to extract representative microstructures with correct geometric features, in particular solid volume fractions. The geometry models need to be more realistic, for example, with respect to surface roughness on the one hand. On the other hand, they should allow for better control of the features and faster generation of more and larger representative volumes. Finally, heat transfer by radiation has to be incorporated into the simulation.
To summarize, the presented study sheds new light on how the microstructure of highly porous wood fiber insulation boards influences their effectivity and points a way to optimize them. The proposed workflow proved to be feasible but needs refinement of the methodology in image processing, modeling, and numerical simulation.
Wood Science and Technology (2023) 57:5-31 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.