Comparative full-field velocimetry of liquid flow within monolithic catalyst carriers via CFD simulations and MRV measurements

In reaction engineering, reactor performance can be improved in many cases by overcoming transport limitations. This requires detailed analyses of transport phenomena in the catalytic beds. Nuclear magnetic resonance (NMR) velocimetry measurements have been utilized for analyzing mass transport of gas flows within opaque monoliths. Comparisons to full-field computational fluid dynamics (CFD) simulations, however, show significant deviations. In this study, polyethylene glycol (PEG) and 3D-printed monoliths including one open-cell foam (OCF) and one honeycomb were used to demonstrate that both operating fluid and monolith morphology influence the achievable signal-to-noise ratio and resolution of NMR data. The velocity profiles measured by NMR in OCF agreed well with full-field CFD simulations with ± 5% deviation. In addition, the similarity between the simulated and experimental velocity fields was quantified by the similarity index, which is 1 for identical images. A mean value of 0.83 was determined for a 10 PPI OCF. Thus, using PEG as the operating fluid and a 10 PPI OCF allows to improve both spatial resolution by 34% and the quality of agreement by 13 percentage points compared to the published results of gas velocimetry within 20 PPI OCF. We further identified and quantified possible sources of deviation between CFD and MRV velocity fields. By limiting our analysis to velocities higher than 45% of the maximum velocity, we could achieve similarity indices of 0.95–0.99.


Introduction
Understanding transport phenomena within monolithic catalyst carriers requires numerical and experimental tools. Morphological properties of monolithic catalyst carriers such as open-cell foams (OCF) and honeycomb structures control the transport phenomena and have a considerable influence on the performance of the reactive systems (Sinn et al. 2021a, b). Generally, OCFs and honeycombs are classified based on the number of pores per inch (PPI) or number of channels per inch (CPI), respectively. The structure of OCFs is determined by their pore size (d p ), window size (d w ), porosity (ε), and specific surface area (S v ), while the morphology of honeycomb monoliths is determined by the channel size (d c ) and wall thickness between two channels (δ wall ).
Heterogeneous modeling using computational fluid dynamics (CFD) simulations and in situ measurement techniques have been widely used to investigate local velocity (Huang et al. 2017;Della Torre et al. 2014), 1 3 138 Page 2 of 12 temperature (Sinn et al. 2021a, b;Ridder et al. 2022), and concentration fields (Ulpts et al. 2017;Della Torre et al. 2016) within catalytic monoliths. There are two options for performing such direct comparison: The first option is to scan the monolith by micro-computed tomography (µCT) imaging, which allows to use its morphology in CFD simulations. The second option is to print an artificially generated monolith, which was used in CFD simulations, using 3D printers. Meinicke et al. (2017) investigated the flow of dimethyl sulfoxide solution within a transparent 20 PPI glass OCF by micro-particle image velocimetry (µPIV) and µCT-based CFD simulations. They qualitatively compared the time-averaged velocity fields from µPIV measurements at three different positions with velocity fields from CFD simulations. While they used the whole geometry of the glass OCF in their measurements, only a part of the OCF-a representative elementary volume (REV)-was considered in their simulations. This neglects the effect of neighboring pores on the velocity field and causes some of the reported discrepancies between numerical and experimental data.
In a study, Cooper et al. (2019) investigated the flow of SF 6 within a diesel particulate filter via 3D CFD simulations and magnetic resonance velocimetry (MRV) measurements. They reported a good agreement between the obtained CFD and MRV mean velocity profiles using a high spatial resolution of the obtained velocity fields (0.14 mm). However, they did not directly compare CFD and MRV velocity pattern. In our previous work (Sadeghi et al. 2020), we performed a full-field analysis of methane flows within 20 PPI OCFs by µCT-based CFD simulations and MRV measurements. The reported simulated and experimental velocity profiles showed an agreement with ± 10% deviation. The low spin density of methane gas led to rather low signal-to-noise ratio (SNR). Thus, using methane as operating fluid limited the spatial resolution of experimental velocity fields to 0.8 mm. In another work Mirdrikvand et al. (2021), we used commercial and 3D-printed honeycombs. We found a positive effect when using the 3D-printed honeycomb, because of its large channel size (2.4 mm which equals three MRV voxel edge lengths) and because it fits ideally into the reactor preventing flow bypasses often encountered in real samples. Due to their short measurement time, gradient-echo sequences with small flip angles and short repetition time T R are the common choice for in-vivo measurements. Also, they enable a low specific absorption rate (heating of tissue). However, disadvantage of this technique is the dependence on the effective transverse relaxation time T 2 * and thus the susceptibility to magnetic field inhomogeneities. The advantage of spin echo (SE) sequences is their independence of T 2 * and therefore better image quality is achieved in case of magnetic field inhomogeneities. Thus, SE sequences are a common choice for measuring flow velocity in porous media (Sederman et al. 1997;Ren et al. 2005;Huang et al. 2017;Sadeghi et al. 2020;Williamson et al. 2020;Clarke et al. 2021).
In this study, we investigate the influence of the operating fluid and the monolith morphology on the quality of agreement between full-field CFD simulations and MRV measurements. For this purpose, we perform velocimetry measurements in a 3D-printed 10 PPI OCF and in a 3D-printed honeycomb using polyethylene glycol (PEG) as operating fluid. A liquid-like PEG-should provide better SNR (and velocity-to-noise ratio, VNR) due to its higher spin density compared to methane gas. In addition, the relaxation time T 1 of PEG is much shorter than water and thus allows for shorter repetition time T R . Additionally, because of its higher viscosity, PEG is less sensitive to pulsatile flow than, e.g., water. A 10 PPI OCF with larger pore and window sizes than the 20 PPI OCF used in our previous work (Sadeghi et al. 2020) is utilized for the velocimetry measurements. The 10 PPI OCF was printed by a 3D printer to ensure good fitting of the sample and to avoid a fluid bypass. The structure of the 3D-printed OCF was obtained by geometrically modifying a commercial sample scanned by μCT imaging. We further identify potential sources of remaining deviations between CFD and MRV velocity fields and quantify their contributions.

CFD simulation
CFD simulations were carried out using the finite-volume CFD software OpenFOAM (version 4.2). The steady-state flow of PEG in this study was considered incompressible since the Mach number (Ma) is lower than 0.3. We used the solver simpleFoam within the OpenFOAM library to solve the governing equations of laminar and steady-state flows (OpenFOAM advanced tutorial 2016). SimpleFoam solves the pressure-velocity equation using the SIMPLEC (semi-implicit method for pressure-linked equationsconsistent) algorithm (Versteeg and Malalasekera 2007). Constant superficial velocities were considered at the inlet (u inlet = constant) and the outlet pressure was set zero (p = 0). A no-slip boundary condition (u = 0) was considered for the tube walls and the surface of monoliths. We considered appropriate entrance length distances (L e ) to ensure the flow has fully developed before reaching the monoliths: (1) L e = 0.056ReD in , with Here, D in is the tube diameter, Re the Reynolds number based on the tube diameter, u inlet the inlet velocity, ρ f the fluid density, and f the dynamic viscosity of the fluid. The Re number calculated based on the pipe diameter are 4, 5 and 5.6 for the measured flow rates. The pore-Re numbers (R p ) were calculated using the pore diameter instead of the tube diameter in Eq. 2 to be 1, 1.1 and 1.3. The flow regime within OCFs can be characterized based on the R p . The examined flows belong to the Darcy-Forchheimer regime which occurs in the range of 1-10 < Re p < 150. In this regime, flow is dominated by inertial forces. The fluid temperature was 20 °C during measurements. The structure of the 10 PPI OCF was reconstructed from micro-computed tomography (μCT) images by using the image processing software ImageJ (version 1.52a, https:// imagej. nih. gov/ ij/) and Meshlab (version 1.3.2, http:// www. meshl ab. net/) as described in our previous work (Sadeghi et al. 2020). The morphology of the 3D-printed OCF was generated modifying the original μCT images in ImageJ to have circular outer geometry and perfect fitting into the reactor which prevents flow bypasses. The geometry of the reactor and the utilized honeycombs were generated by a CAD software (Freecad, version 0.18.4, https:// www. freec adweb. org/). The computational mesh was generated via the commercial meshing software cfMesh (version 1.1.2, https:// cfmesh. com/). Grid dependency tests revealed computational networks with a minimum of 9.4 and 4.7 million cells were required for OCF and honeycomb structures, respectively. Grid independency was investigated by performing simulations with different cell numbers. In these simulations, the initial velocities were of 23.94 mm/s (705 mL/min) and 16.97 mm/s (0.5 mL/ min) for honeycomb and OCF, respectively. Afterward, the pressure loss caused by each monolith was calculated as a function of cell number to determine the number of cells at which pressure loss becomes independent from cell number. Table 1 presents the morphological properties of the OCF determined by analyzing the μCT images. To determine the open porosity, we considered a cylindrical volume with equal length and diameter as the representative volume element positioned in the center of the OCF. We increased the diameter of the cylindrical volume and calculated the open porosity for each diameter of the corresponding volume of the cylinder. The obtained open porosity reached a stable value at diameter of 11 mm after initial fluctuation for the small representative cylinder. Figure 1 represents the schematic illustration of the adopted method for calculating open porosity. Afterward, the distribution of both pore and window size were analyzed within the chosen diameter of the cylindrical volume. Subsequently, we could determine the averaged values of pore and window size and their range of variation (Kiewidt 2017). The utilized honeycomb has a channel size of d c = 2.4 mm and a wall thickness between two channels of δ wall = 0.8 mm.

MRV measurements
A spin echo (SE) and phase contrast (PC)-based (MRV) pulse sequence was implemented (Fig. 2). In the employed non-slice-selective 3D MRI sequence, spatial resolution was achieved by one read (along Z-axis) and two-phase encoding gradients (along X-and Y-axis). One-directional velocity encoding was performed using a 5-step scheme as follows. To improve VNR, a dual-velocity encoding (VENC) method was used (Lee et al. 1995;Nett et al. 2012;Schnell et al. 2017;Callahan et al. 2020). Aliasing in the low-VENC dataset was corrected by using a second dataset with a higher VENC value (see Table 2). For further VNR improvement, this high VENC dataset was generated similar to the technique presented in Johnson and Markl (2010). Two steps with velocity encoding gradients applied in opposite direction were done for both VENC values and an additional step was done without velocity encoding. More details about the encoding scheme are described in the supplementary information.
All measurements were performed on a horizontal 7 T scanner (BioSpec 70/20 USR, Bruker BioSpin MRI, Ettlingen, Germany), which was equipped with a 114-mm bore gradient system (BGA 12S2) with maximum gradient strength of 440 mT/m and a maximum slew rate of 3440 mT/(m•ms). ParaVision 5.1 was used as the user interface to perform the measurements. The post-processing of the data was performed using a self-developed MATLAB script (R2021a MathWorks, Natick, MA, USA). A horizontal 72-mm bore birdcage 1 H quadrature transceiver RF coil (Bruker BioSpin MRI, Ettlingen, Germany) was used in all measurements. The field of view (FOV) was 4.2 × 4.2 × 8.4 cm 3 (honeycomb) and 4.2 × 4.2 × 6.3 cm 3 (OCF) with a matrix size of 60 × 60 × 120 (honeycomb) and 80 × 80 × 120 (OCF) resulting in a spatial resolution of 0.700 mm/pixel and 0.525 mm/pixel. T R /T E was 50 ms/6 ms, FA (flip angle) of the rectangular excitation pulse α = 145° and the acquisition time for each echo was 2.4 ms. Flow encoding gradients had a duration of δ = 370 µs and a delay of Δ = 3.315 ms. The gradient amplitude was set depending on the desired VENC value. For each step, four accumulations were done and averaged to increase the SNR. Apodization was done by using a Hanning window to reduce Gibbs ringing. A binary mask was used based on Otsu thresholding (Otsu 1979) applied to magnitude images. This mask was applied to the velocity image to distinguish the fluid phase from the Fig. 2 Implemented SE PC MRV pulse sequence for 3D imaging and one-directional velocimetry (not to scale). Spatial encoding was done by gradients in read direction (along Z-axis) and two phase-directions (along X-and Y-axis). Velocity encoding gradients with duration δ and delay ∆ were applied. Rectangular excitation pulse was applied with angle α = 145°. Spins were refocused by a rectangular 180° pulse after duration τ. Spin echo was acquired after an additional τ. Crusher gradients at the end of the sequence were used to dephase and suppress unwanted transverse magnetization background (monolith structure). If the voxel size is much smaller than the pore size, fraction of partial volume effects are negligible, i.e., fluid phase and background are clearly distinguishable (Bruschewski et al. 2021). However, the OCF also had pores that were smaller than the voxel size. In this case, accurate segmentation was not possible.

Method of CFD and MRV comparison
The results of CFD simulation and MRV measurement were compared in terms of averaged axial velocity profiles and cross-sectional axial velocity fields. Figure 3 represents cross sections of the utilized monoliths at an arbitrary axial position which were obtained from μCT images of OCF or CAD file of the honeycomb, CFD simulations and MRV measurements. Parity plots including prediction bands were used for better quantification of agreement between the averaged values of velocity in each axial cross section. These bands were calculated by Here, m is the slope of the regression line (the intercept for all the regression lines is zero), T INV is the calculated so-called student's distribution for the datasets, SEM is the standard error, n is the of number the datapoints in the sample, X m is the averaged values of X and SSX is the sum of squared deviation of X values from the averaged value. In this equation, X and Y correspond to the MRV and CFD values, respectively.
To quantify the comparison between CFD and MRV velocity fields shown in Fig. 7, we compared both velocity fields voxel by voxel using a two-dimensional correlation coefficient also known as similarity index between two matrices with the same size. For that we used the corr2 function within the MATLAB image processing toolbox (R2021a MathWorks, Natick, MA, USA) as reported in Sadeghi et al. (2020). This approach comprised a voxelization procedure. CFD and MRV velocity fields were converted into the grayscale images with the same size and were considered as two matrices of A and B, respectively. Afterward, the similarity index was calculated by Here, A and B are the mean values of A and B, respectively.

Experimental setup
PEG 400 was purchased from Carl Roth GmbH + Co. KG (Karlsruhe, Germany). According to manufacturer specifications, PEG has a molar mass of M = 380-420 g/mol, mean chain length of n ≈ 8, and viscosity of 105-130 mPa·s. Flow was realized by using a peristaltic pump (Standard Digital Drive with Easy-Load II Pump Head; Masterflex, Vernon Hills, USA). The pump tubing was also from Masterflex and flexible PVC tubing was used for the connection to (4) Fig. 3 Cross sections from μCT images/CAD file, CFD simulations and MRV measurements of a OCF and b honeycomb, respectively the glass flow reactor (Fig. 4). The structure of the utilized honeycomb was printed from polylactide (PLA) by a fused deposition modeling (FDM) 3D printer (Ultimaker 3, Ultimaker BV, Utrecht, Netherlands) with the nozzle diameter of 0.4 mm and a resolution of 12.5 × 12.5 × 2.5 µm 3 (layer height in vertical direction: 200 µm). The OCF structure was printed by a Digital Light Processing printer (Elegoo Mars, https:// www. elegoo. com/, Shenzhen, China) with 2 K resolution (12.5 × 12.5 µm 2 ) in XY direction and layer height (Z) of 50 µm. Measurements for the OCF were done at (500 ± 6) mL/min, (625 ± 7) mL/min and (705 ± 8) mL/min. The honeycomb was measured at (530 ± 6) mL/min. The magnetic susceptibility of PEG and PLA were 0.36 and 0.21 ppm, respectively (Wapler et al. 2014;Souza et al. 2013).

Results and discussion
Full-field CFD simulations of PEG flow at three different flow rates of 500, 625, and 705 mL/min through a 3D-printed OCF were compared with results of MRV measurements in terms of velocity fields and averaged velocity profiles. To compare CFD and MRV data quantitatively, the averaged axial velocity profiles along the reactor were calculated for both simulations and experiments. For this purpose, 81 planes (slices) perpendicular to the flow direction at different axial positions were considered. The 81 planes included 43 slices within the OCF and 38 slices at axial positions before and after the OCF. The average of the axial component of the velocity (u z ) was calculated for each of the 81 planes. The velocity profiles from CFD and MRV agree well (Fig. 5). When comparing each value of the slice-averaged velocity between CFD and MRV, the deviation is within ± 5% at all flow rates (Fig. 6). The quality of agreement was improved in comparison with the flow of methane within 20 PPI OCF, which was ± 10% (Sadeghi et al. 2020). The shaded bands in Fig. 6, which were calculated as described in Sect. 2.3, present the predicted interval at the confidence level of 95%. The R-squared values of the regression line (R 2 ) and the calculated statistical values are shown in Table 3 for all flow rates. For the flow in the 3D-printed honeycomb, the obtained averaged velocity profiles agree well (see Fig.  S1 in supplementary information). The averaged relative error between CFD and MRV velocity values within honeycomb was 2.54%. This value is lower than the 4.17% reported in our previous publication (Mirdrikvand et al. 2021) in which we investigated flow of methane in the same 3D-printed honeycomb. The calculated averaged flow rates within OCF from MRV velocity data were 496, 617 and 693 mL/min while the corresponding value for flow within honeycomb was 537 mL/min. This compares well with the velocity measured by the flow meters during the experiment ((500 ± 6) mL/min, (625 ± 7) mL/min, (705 ± 8) mL/min, and (530 ± 6) mL/min)).  Fig. 5), we applied a threshold filter removing voxels with a velocity lower than a certain ratio of the inlet velocity. For this purpose, different thresholding values were examined to determine the best thresholding value by calculating their corresponding averaged relative deviation between MRV and CFD velocity values within the OCF (Fig. S4). The optimum percentage ratio of the thresholding value to the inlet velocities were 7.7, 7.7 and 7.5 for flow rates of 500, 625 and 705 mL/min, respectively. Afterward, the open porosity at each MRV slice was calculated by dividing the area of the remaining fluid region by the whole area of the slice. An axial profile of the resulting porosity from CFD and MRV is shown in Fig. 8. It should be noted that by applying such a thresholding filter, vortices and regions with negative velocity values are also removed from the velocity fields. Most of these vortices are located near the wall and at the entrance and terminal region of the OCF. This removal influences the quality of agreement between CFD and MRV. Furthermore, applying such filter causes the elimination of common voxels in the nearwall region (including tube wall and OCF surface) which  belong to both fluid and solid. This influence is shown by plotting velocity profile over an arbitrary line in several slice before and after applying the threshold filter (see Fig. S5).
To quantify the agreement between CFD and MRV velocity fields, CFD velocity fields were converted to the same size of the corresponding velocity fields determined by MRV. The obtained voxelized velocity fields are shown in Fig. 9 for the same slices as shown in Fig. 6 (flow rate of 705 mL/min). The voxelized velocity fields for flow rates of 625 and 500 mL/min are shown in Figs. S8 and S9 of the supplementary information. To perform a voxel-wise comparison, the resolution of the CFD velocity fields was decreased by convolution to the same resolution obtained by MRV. This method helped us to determine the angular deviation between CFD and MRV velocity fields which is a demanding task with the naked eyes due to complexity of the OCF structure. Therefore, the similarity index (result of the corr2 function) between MRV velocity fields and corresponding CFD velocity fields at each slice was calculated by shifting one data set in axial position and rotating it from 0° to 360°. By maximizing the similarity index as a function of angular rotation and axial deviation we could align the velocity fields and find the correct spatial displacement between both datasets (Fig. 10). As shown in Fig. 10a, the highest structural similarity was obtained at an optimum rotational angle at each slice. The optimum rotational angle per slice and the corresponding similarity index is shown in Fig. 10b (flow rate of 705 mL/min within 10 PPI OCF). The optimum rotational angle varies between 175° and 177°. With a step width of the rotation angle of 1°, it was 176° in most slices. These deviations between the highest similarity indices at the slices and the corresponding similarity indices at 176° were less than 1% and can be considered as a negligible measurement error. Therefore, 176° was used as the optimal angle of rotation. The mean of the obtained values of the similarity index, which is 1 for the identical images, was 0.83. The mean similarity index of PEG flow within the 3D-printed 10 PPI OCF in this work was much higher compared to reported mean similarity index for flow of methane within commercial 20 PPI OCFs, which was 0.7 (Sadeghi et al. 2020). This considerable improvement was achieved firstly by using a more suitable operating fluid (PEG), which enables a better SNR and allows us to achieve a better spatial resolution (0.525 mm) compared to methane Fig. 8 The obtained porosity within 10 PPI OCF from MRV measurements and CFD simulations for three flow rates of 500, 625, and 705 mL/ min Fig. 9 The voxelized CFD and MRV velocity fields corresponding to velocity fields shown in Fig. 7 (0.8 mm). Secondly, using 10 PPI OCF with larger pore and window size (d p = 5.76 mm, d w = 3.3 mm) compared to 20 PPI OCF (d p = 3.41 mm, d w = 1.90 mm) allowed to resolve the structure of OCF in MRV measurements. Thirdly, 3D printing of a 10 PPI OCF from PLA with perfectly circular outer geometry provided better fitting of the sample into the reactor and prevented the problems with flow bypass during measurements. Fourthly, using an improved 5-step multi-VENC encoding (see supplementary information) improved the VNR in comparison to single-VENC encoding, which was used in our previous publication (Sadeghi et al. 2020). The optimum rotational angle and the similarity index between CFD and MRV data for the flow rates of 625 and 500 mL/min are shown in Figs. S6 and S7 of the supplementary information.
Although the mean similarity index of CFD and MRV velocity fields showed significant improvement, there are obviously still several sources of deviation which limit the similarity index. For example, convolving CFD velocity fields (with pixel size of 200 microns) to the size of MRV velocity fields (with pixel size of 525 microns) influenced the original structure of the velocity fields and thus decreases the highest possible similarity index. To quantify the contribution of this convolution procedure to the similarity index, the size of utilized computational network for CFD velocity fields were convolved to the same voxel size as the MRV velocity fields. Afterward, the resolution of convolved CFD data was increased by linear interpolation to the same resolution as the original CFD data. Finally, the similarity index of the original CFD data and the interpolated CFD data was determined. The original CFD velocity fields and their corresponding convoluted, interpolated velocity fields and their similarity index are shown in Fig. S10. It was concluded that convolving CFD velocity fields reduces the similarity index by 2.0-2.2%.
In addition to the convolution effect, the MRV velocity fields included noise, which can be characterized by the standard deviation σ v . The contribution of MRV noise to the reduction of the similarity index was estimated by adding noise to the CFD axial slices using the same standard deviation determined in the corresponding slices of the MRV data. The similarity index between the noisy CFD data and original CFD data was determined (see Fig. S11). The experimental MRV noise reduced the highest achievable similarity index by 1.0-1.2%.
A potential axial shift and angular rotation between CFD and MRV data were considered in the alignment of their velocity profiles and velocity fields. Although the resolution of CFD velocity fields was much higher compared to those obtained by MRV, they were not continuous and the step size of CFD slices in axial direction was limited to their grid size. The similarity index calculated from two consecutive CFD slices is shown in Fig. S12. The value of 0.985 is below 1 which confirms that a slight mismatch in axial alignment of CFD and MRV can decrease the highest achievable similarity index by about 1.5%.
To quantify if regions with fast and slow velocities play a different role in determination of the similarity index, a threshold filter was applied to both CFD and MRV data. We applied a binary mask to both CFD and MRV datasets with a threshold value of 45% of the maximum velocity encountered in each respective CFD slice. Velocities below the threshold value were set to zero. We chose 45% as this gives us the highest similarity index, it thus appears to be the optimal value. Subsequently, the similarity index of the filtered CFD and MRV velocity fields was calculated. By using only voxels with fast velocities for the illustrated velocity fields in Fig. 5, the similarity index improved by 7-11 percentage points (see Table S2 and Fig. S13). Several reasons may be responsible for the lower similarity index for slow velocities (respectively higher similarity index when only considering fast velocities). In the experiments, it could not be guaranteed that PEG was completely free from air bubbles which influenced flow behavior. For fast velocities they 138 Page 10 of 12 may be easily removed, while for regions with slow velocities it was more likely that air bubbles remain at the surface of the OCF and thus affect the similarity index more. The effect of these air bubbles could not be quantified because their size was not large enough to be directly resolved in the MR images. Air bubbles will reduce the signal in the voxels, however, these minor signal losses could not be distinguished from the partial volume effect as air bubbles occur at the surface of the structures. The comparison between the obtained velocity fields from CFD simulation and map of the absolute signal from MRV measurements confirms the hardships in identifying air bubbles (Fig. S14). The effect of noise on the similarity index was already described before. It should be noted that this noise causes a bigger relative error for slow velocities. In addition, the effect of consecutive slices on similarity index was already described (Fig. S12). This effect may have a larger influence for slow velocities. Slow velocities occurred mainly at the surface of the OCF. Hence, it was possible that for two consecutive voxels one contained a velocity value while the other one was inside of OCF structure and had a value of zero. This effect was not present for fast velocities because they did not occur at the surface. As discussed, when increasing the threshold value above the optimum value (45% of the highest velocity), the similarity index decreased. This can be described by the contribution of displacement errors in MRV data, which increases with the velocity and thus decreases the similarity index. Displacement errors in the MRV data are discussed in John et al. (2020). During the MRV measurements, a displacement of spins from one to another voxel can occur and thus velocities would be assigned to a wrong position. The maximum displacement (John et al. 2020) is calculated by: Here t = 4.033 ms is the time between the center of the velocity encoding gradients and the center of acquisition. For foams with flow rates of 500 mL/min (corresponding to v max = 83mm∕s ), 625 mL/min ( v max = 111mm∕s ), and 705 mL/min ( v max = 129mm∕s ) a maximum displacement of 0.64, 0.85 and 0.99 voxels was calculated. Even though these values are below 1, displacement errors may occur, e.g., for spins located at the border of a voxel, which need to move a shorter distance to reach next voxel. Therefore, displacement errors may contribute to the reduced similarity index.
Potentially, all sources of error discussed above could influence the similarity index by 12-16% (see Table S1). Additionally, there were other sources of deviation like fluctuation of the utilized peristaltic pump and mass flow controllers, which may interrupt steady flow within monoliths. Further, in pores small compared to the voxel size, segmentation of velocity fields (described in Sect. 2.2) did not accurately distinguish the fluid phase from the monolith structure. Furthermore, the utilized glass reactor in MRV measurements may not be perfectly parallel to the flow direction and leading to slightly tilted MRV velocity fields. Despite the improved similarity there is still a systematic offset, as can be seen in Fig. 8, which shows that the MRV porosity data is lower compared to the calculated ones, with an averaged and maximum offset of 4% and 10%, respectively. As discussed, the porosity from MRV data was calculated after applying a threshold filter equal to 8% of inlet velocities. This threshold value was determined to achieve the highest agreement between CFD and MRV velocity profiles shown in Fig. 5 which is 95%. Judging from this, the obtained CFD and MRV porosities have a higher deviation compared to their velocity profiles and velocity fields.

Conclusion
In this work, the flow of PEG within catalytic monoliths was analyzed by full-field CFD simulations and MRV measurements. Using PEG as operating fluid allowed to obtain higher spatial resolution of MRV velocity fields in comparison to gaseous methane, which was used in our previous work. In the present work, the investigated monoliths with perfectly round outer geometry were generated by 3D printers to avoid measurement errors introduced by flow bypass and back flow caused by geometrical artifacts. Additionally, here we use a 10 PPI OCF, which has larger pore and window diameters compared to the 20 PPI OCF used in our previous work. This allowed to better resolve the foam morphology by MRV measurements. The CFD and MRV averaged velocity profiles within 10 PPI OCF showed excellent agreement with a precision of 95%, which was 5% points higher than the best achieved agreement reported for velocimetry of methane flow within 20 PPI OCF in our previous work. In case of flow within the 3D-printed honeycomb, the averaged velocity profiles showed excellent agreement as well with a precision of 97%. The good agreement between the velocity fields from CFD and MRV data was quantified by determining their similarity index. The attained similarity index, which is 1.0 for identical images, showed values around 0.83. Several error sources and their contribution to the similarity index were determined. The total contribution of these sources into the reduction of similarity index were quantified to be between 12 and 16%. Therefore, removing such sources of error would allow us to achieve higher similarity indices from 0.95 to 0.99. There were other sources of error such as fluctuation of the utilized pump and displacement of spins through voxels during measurements, but we could not determine their contribution on the similarity index.
The results of this work underpin the reliability of the MRV technique as a tool for velocimetry of fluid flows