20th OpenFOAM Workshop Book of Abstracts

Compiled and edited by Mirza Popovac 20th OpenFOAM Workshop Book of Abstracts 30. June – 4. July 2025 Vienna, Austria doi: 10.5281/zenodo.17560862 1 Preface This Book of Abstracts compiles the submissions presented at the 20th OpenFOAM Workshop. Abstracts appear as provided by the authors, with only minimal formatting and organizational editing. As chairman of the 20th OpenFOAM Workshop, I had the pleasure and honor of organizing the largest and longest-running OpenFOAM® conference, held in Vienna from June 30th to July 4th, 2025. This special anniversary edition marked two decades of a tradition that continues to bring together the global OpenFOAM® community. For the first time, the workshop took place in Austria—and fittingly, in its capital, Vienna—a city renowned for its cultural richness, creative energy, and its setting along the Danube, where history and innovation converge. Hosting the event here felt both timely and symbolic, highlighting Austria’s growing role in open-source computational fluid dynamics (CFD) and the collaborative ethos that has shaped the OpenFOAM® community over the years. The Austrian Institute of Technology – AIT, as host of the 20th OpenFOAM Workshop, embodied this spirit remarkably well. As Austria’s largest non-university research organization, AIT’s mission blends scientific excellence with societal impact. By fostering open-source approach in research and technological development, AIT provided an ideal setting for a workshop dedicated to the world’s most widely adopted open-source CFD platform. Through this gathering, AIT proudly showcased not only its own advancements in simulation and modeling but also the breadth and vitality of Austrian research in fluid mechanics and its diverse applications. This year’s program reflected that diversity: prominent professors from Austria’s top universities shared national expertise in multiphase, particulate, and laser-induced flow phenomena. Meanwhile, esteemed keynote speakers from institutions across Europe, America, and Asia enriched the global dialogue, oering perspectives on topics ranging from machine learning in fluid dynamics to the evolution of CFD in fields as varied as maritime engineering and biomedical flows. Collectively, these contributions underscored the field’s ongoing expansion—from fundamental physics to data-driven innovation. The workshop welcomed over two hundred participants from academia and industry, promoting a truly international and interdisciplinary exchange of ideas. With technical sessions, poster presentations, and a dedicated high-performance computing session, the event highlighted how OpenFOAM® thrives at the intersection of community, computation, and creativity. Beyond the formal program, the workshop oered a valuable chance to reconnect in person, share experiences, and reflect on the remarkable progress of OpenFOAM® and the open-source movement through years of collaboration. My sincere thanks go to all participants, contributors, sponsors, and volunteers who made this event possible. Organizing the 20th OpenFOAM Workshop has been a privilege and a rewarding opportunity to celebrate not only technical progress but also the strength of an international community united by curiosity, transparency, and shared purpose. The 20th OpenFOAM Workshop was the perfect opportunity to connect, learn, and celebrate the spirit of innovation together! On behalf of the OpenFOAM Workshop Committee, I thank you all for your contribution to this exciting event. Mirza Popovac Conference Chairman 2 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NUMERICAL ANALYSIS OF HIGH-PRESSURE CAPILLARY VISCOMETER MOHAMMAD POURHOSSEINIAN1,2*, BAHRAM HADDADI1 1 Competence Center CHASE GmbH, Linz, Austria 2 TU Wien, Institute of Chemical, Environmental & Biological Engineering, Vienna, Austria * [email protected] Keywords: High-pressure capillary viscometer, non-Newtonian fluids, Heat dissipation, OpenFOAM® High-pressure capillary viscometers are essential tools for characterizing the rheological properties of non-Newtonian fluids under controlled pressure and temperature conditions. These devices are particularly useful for understanding the viscosity behavior of polymers, elastomers, and other complex materials under high shear rates. These non-Newtonian fluids exhibit viscosity changes based on shear rate and temperature variations. To describe this behavior accurately, established models such as Power Law, Cross, and Carreau with their temperature extensions are employed. Using numerical simulations in OpenFOAM®, it is possible to capture the effects of shear-rates, temperature dependence, and viscous heating, providing insights into the material characterization devices such as High-pressure capillary viscometers. In this work, a new incompressible solver for OpenFOAM® v2006 is developed. Power-law Arrhenius, CarreauArrhenius, and Cross-Arrhenius viscosity models are integrated to effectively capture the shear-dependent and temperature-dependent behavior of shear-thinning non-Newtonian fluids [1]. The viscosity parameters are modularly defined within the transport properties dictionary, allowing adjustments to key parameters and viscosity limits for each model. These implementations enhance OpenFOAM®’s usability for engineering applications requiring precise rheological modeling. Additionally, viscous dissipation is incorporated to account for heat generation due to shear heating. The solver has been rigorously validated against analytical solutions and experimental data, ensuring its accuracy and robustness. Figure 1a illustrates the computational geometry of high-pressure capillary viscometer setup consisting of piston section (upper part of the geometry) and outlet section (die); in a second geometry, the last 3/4 of the die section is considered for analytical validation. The die has a length of 20 mm and a diameter of 2 mm, giving a length-to-diameter (L/D) ratio of 10. The piston-cylinder section has a diameter of 15 mm. A fixed velocity corresponding to the piston speed is applied at the inlet, while the outlet is set to atmospheric pressure. The walls are maintained at a fixed temperature for consistency with measurements, and no-slip conditions are applied. Computational geometry is prepared using OpenFOAM®’s blockMesh utility, which can be seen in Figure 1b. After a detailed mesh study, an optimally structured grid consisting of 44,400 cells provides a good balance between accuracy and computational efficiency. Figure 1: (a) Geometry and (b) generated mesh with side, top and bottom view of the die section The power law Arrhenius model with n = 0.139, k = 16526 Pa·s0.139 and E = 6023 J/mol is used to model the nonNewtonian behavior of the elastomer for the simulations. These parameters are extracted from the same measurement series that are used for comparison. Additionally, dissipation effects are included to account for viscous heating, significantly improving the accuracy of simulations where temperature variations impact viscosity and flow characteristics. a) b) 3 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria To validate the solver, two separate validation cases at different shear rates were performed: (1) Analytical validation: The pressure drop in the last 3/4 of the die was compared against analytical predictions following equation 1 which describes a fluid that obeys power law behavior in a circular pipe [2]: ∆𝑃=2𝑘𝐿 𝑟[�D.3+1 𝑛 ⁄ 𝜋𝑟3]𝑛 (1) (2) Experimental validation: The full high-pressure capillary viscometer (HPCV) setup was simulated, and pressure drop results were compared with experimental data. Figure 2 presents a comparative analysis between experimental measurements, analytical predictions, and CFD results, demonstrating good agreement. Figure 2: Simulation pressure drop results comparison with (a) Analytical and (b) Experimental data in 80 °C Additionally, Figure 3 provides contour plots of pressure, viscosity, and temperature within the system, illustrating the impact of shear-thinning effects and viscous heating. The incorporation of viscous dissipation allows for the evaluation of temperature variations, which is particularly important in measuring output temperature or in applications where temperature changes significantly impact viscosity and flow behaviour. The implemented solver accurately captures pressure drop, viscosity variations, and heat dissipation, confirming the solver’s accuracy in predicting non-Newtonian flow behavior. This result confirms the robustness and reliability of the proposed implementations, making them valuable additions for the studies with high-pressure capillary viscometers and complex non-Newtonian fluid simulations. Figure 3: Pressure, dynamic viscosity and temperature contours at 80 °C and 100 1/s Acknowledgements The authors acknowledge financial support through the COMET Centre CHASE, funded within the COMET − Competence Centers for Excellent Technologies programme by the BMIMI, the BMWET and the Federal Provinces of Upper Austria and Vienna. The COMET programme is managed by the Austrian Research Promotion Agency (FFG). The authors would like to express their gratitude to Semperit Technische Produkte GesmbH, especially Silke Koch and Florian Arthofer for providing the necessary data and financial support for this research. References [1] V. Kimmel, E. Ercolin, R. Zimmer, M. Yörük, J. Winck, and M. Thommes, “Measuring and Modeling of Melt Viscosity for Drug Polymer Mixtures,” Pharmaceutics, vol. 16, no. 3, p. 301, Feb. 2024, doi: 10.3390/pharmaceutics16030301. [2] H. A. Barnes, A handbook of elementary rheology. Aberystwyth: University of Wales Institute of Non-Newtonian Fluid Mechanics, 2000. 10 100 1000 10 100 1000 10000 Pressure Drop [bar] Shear Rate [1/s] Analytical Simulation 10 100 1000 10 100 1000 10000 Pressure Drop [bar] Shear Rate [1/s] Experiment Simulation 4 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria BENCHMARKS IN THE TENDERING PROCESS FOR HPC CLUSTERS CARSTEN THORENZ1, FABIAN BELZNER2 1Federal Waterways Engineering and Research Institute, Karlsruhe, Germany, carsten.thoren[email protected] 2Federal Waterways Engineering and Research Institute, Karlsruhe, Germany, [email protected] Keywords: HPC, Cluster, Benchmark, CFD, IO The Federal Waterways Engineering and Research Institute (BAW) is performing consultancy work for the German Federal Ministry of Transport and Digital Infrastructure and the executing bodies of the ministry. For this purpose, highperformance computing (HPC) systems have been used for several decades, starting in the 1990s with Cray multiprocessor systems and have continually advanced since that time. As a public entity, the BAW must comply with strict rules in the tendering process. The following provides an overview on the benchmarks used in the past tendering processes and the lessons learned from it. Particular emphasis is placed on the role of OpenFOAM® in this process. The BAW uses HPC resources mainly for computational fluid dynamics simulations of various kinds. These range from large-scale geophysical models that cover the North Sea with the estuaries of northern Germany to small-scales such as gap flows measuring just a few centimetres. As a consequence, different numerical codes are used on the same computer systems. In the following, only the small-scale CFD models will be regarded. In the BAW, several CFD codes with a connection to the work of Professor Gosman’s group at Imperial College London have played and continue to play an important role. This started with “Comet” in 1999, developed by Professor Peric, which was later superseded by StarCCM+ in 2005. After a short intermezzo with the university code NaSt3DGPF, which is based on cartesian grids, the first trials with OpenFOAM® began in 2009. OpenFOAM® was slowly introduced into the production work in the following years, after substantial effort was put into making it suitable for waterway consultancy. A major effort was the implementation of a suitable set of boundary conditions for hydraulic structures engineering. These were recently published [1] and are now available to the OpenFOAM® community. The available computing resources at BAW (Karlsruhe) have continuously increased over the past three decades, doubling the number of CPU cores roughly every 2.5 years (Figure 1). Figure 1: Installed HPC CPU cores in BAW (Karlsruhe branch) This evolution has been accompanied by architectural changes: starting with vector processors in 1995, evolving into large shared memory machines, and, since 2007, into distributed memory cluster systems. As hardware advanced, application software progressed in parallel. Consequently, the software used for benchmarking in the tender process had to be adapted. Over the years, a mix between synthetic benchmarks and application benchmarks has always been used. As a long-standing standard LINPACK [2] was used. While not matching the application's performance, the LINPACK performance has the significant advantage to deliver consistent and comparable results over decades. Later file system tests and more recently HPCG were introduced into the benchmark suite. On the application side, OpenFOAM® was introduced into the mix in 2010 - at first only as a compatibility test, but since 2012 also as a benchmarking test for the 1 10 100 1000 10000 100000 1995 1999 1999 2005 2005 2007 2010 2012 2015 2017 2020 2023 2025 Cores Year 5 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria tendering process. The usage of OpenFOAM® in the benchmarks has revealed some difficulties in the last decade. The benchmarking is based on the interFoam solver, as this is the primary solver used at BAW. Due to required changes of the versions used - triggered by the close link between the OpenFOAM® version and the version of the compiler toolchain - results are barely comparable. We started with interFoam in version 1.6.1, which used a “p” formulation, changed than to 2.3.1, which used a “p_rgh” formulation, and later versions saw further changes in the algorithm. Thus, the benchmark results lost their direct comparability. In the benchmarks for OpenFOAM®, the tendering vendors were always requested to perform tests on both the throughput as well as scalability of the jobs. Additionally, BAW performs extended in-house tests after the commissioning of the systems. These tests are not targeted on finding the peak performance, but at identifying the highest throughput. To this end, parallel efficiency is used as a measure. Over the last 15 years, we have observed a significant drop in the number of cells per MPI partition that results in the optimal efficiency. Starting with approximately 200,000 cells/partition in 2010, we now see the optimum at around 20,000 cells/partition (depending on the topology of the problem). It is notable that the scaling performance depends significantly on the preconditioner used (AMG, DIC) and on the linear solver (PBiCGStab, AMG). Our tests have shown that the apparent scaling (in terms of speedup) is better for the DICpreconditioned PBiCGStab, resulting in the highest peak performance. However, the performance of the AMGpreconditioned PBiCGStab is only slightly worse at high core counts, while it performs better at lower core counts. Thus, AMG is the preconditioner of choice at BAW. Apart from the CPU performance, the input/output (I/O) performance is also relevant for simulations. Here, we observe several trends. Firstly, ever-growing simulation sizes. Secondly, a shift from RANS-based simulations, which require little output, to DES or LES-based simulations, which require the output of much more data. Together with the increasing computing power available and the uncollated writing format, this results in a spectacular increase in the number of files written in a typical simulation. Unfortunately, the performance of the used file systems has not increased in the same way as CPU performance. We observed an increase in LINPACK performance by a factor of 100 over the last 15 years. For the same machines, I/O throughput has increased by a factor of only around 3-4 - significantly less. Even worse is the situation for OpenFOAM®: Here we additionally observe an increase in the number of written files, as the optimal number of cells per MPI partition is decreasing, as mentioned above. Thus, even for the same jobs, the number of written files has increased by a factor of 10 in the last 15 years. This increase must then be multiplied with the increase in the computing performance (which is not precisely known for OpenFOAM®, but is approximately 100 for LINPACK) and set in relation to the increase in the metadata performance of the file system. For the latter, we observed an increase by a factor of approximately 3 - several orders of magnitude lower than the increase in output generated by the CFD simulations. This effectively makes the file system metadata throughput a significant bottleneck for the use of OpenFOAM®. References [1] C. Thorenz, “Boundary Conditions for Hydraulic Structures Modelling with OpenFOAM” in Proceedings of the 10th International Symposium on Hydraulic Structures (ISHS 2024), pp 76-86, ETH Zürich, 2024, Available: https://doi.org/10.3929/ethz-b-000675949 [2] J. J. Dongarra, “Performance of Various Computers Using Standard Linear Equations Software”, University of Tennessee, Knoxville TN 37996-1301, Oak Ridge National Laboratory, Oak Ridge TN 37831, CS-89-85, 2011, Available: http://www.netlib.org/benchmark/performance.ps 6 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria DEVELOPMENT OF A VOLUME-AVERAGED DIRECT SURFACE DESCRIPTION METHOD FOR MODELING WAVE INTERACTION WITH POROUS STRUCTURES MOHSEN FARHADZADEH1, ERIK DAMGAARD CHRISTENSEN2, 1Department of Civil and Mechanical Engineering, Technical University of Denmark, [email protected] 2Department of Civil and Mechanical Engineering, Technical University of Denmark, edc[email protected] Keywords: Porous medium, Turbulence, Free Surface Modeling, Floating Wind Turbines, Nature-Inclusive Design Modeling the interactions between surface waves and marine structures with porous components is essential for predicting hydrodynamic behavior, structural loads, and energy dissipation across a wide range of coastal and offshore engineering applications. Porous materials are increasingly used to enhance the performance and environmental integration of marine systems, from breakwaters and scour protection layers to innovative nature-inclusive designs for floating offshore wind turbines (FOWTs)[1]. This study develops and implements a robust numerical framework that extends the Direct Surface Description (DSD) method [2] by incorporating the volume-averaging technique[3] to accurately simulate wave–porous media interactions. The new volume-averaged DSD solver integrates sharp-interface free-surface tracking with macroscale averaged momentum equations that account for porosity, permeability, and resistance effects, while efficiently solving only the water-phase flow. The volume-averaged DSD framework is validated through benchmark comparisons with existing OpenFOAM methods (interFoam, interIsoFoam) and demonstrated in two-dimensional and three-dimensional showcase cases featuring more complex and realistic scenarios. The focus of this study is on the development, validation, and performance assessment of the volume-averaged DSD model, which provides an advanced simulation tool for a broad range of marine engineering problems involving structures with porous components and wave interactions. Acknowledgments This study is supported by the INF4INiTY project, funded by the European Union through Horizon Europe under Grant Agreement No. 101136087. References [1] C. Windt, N. Goseberg, S. Schimmels, H. Rusche, F. Adam, W. Swidzinski, K. K.-F. Kazimierowicz-Frankowska, V. Kirca, B. Sumer, T. Petersen, N. Rojas, M. Penalba, G. Bracco, G. Girogi, P. Troch, M. Streicher, E. H.-W. HalvorsenWeare, H. Braaten, E. Christensen, and J. Sørensen, Integrated designs for future floating offshore wind farm technology – Towards nature inclusive innovations. CRC Press, 2024, p. 775–783. [2] J. R. K. Qwist and E. D. Christensen, “Development and implementation of a direct surface description method for free surface flows in openfoam,” Coastal Engineering, vol. 179, 2023. [3] Y. Zhai, D. R. Fuhrman, and E. Damgaard Christensen, “Numerical simulations of flow inside a stone protection layer with a modified k-ωturbulence model,” Coastal Engineering, vol. 189, 2024. 7 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria COMPARISON OF HODMD MODES OBTAINED FROM NUMERICS AND EXPERIMENT ON THE EXAMPLE OF CYLINDER IN CROSSFLOW AT Re ≈5000 MAREK BELDA1, TOMAS HLAVATY2, MARTIN ISOZ3 1Department of Fluid Dynamics and Thermodynamics, Faculty of Mechanical Engineering, Czech Technical University in Prague, Technick´ a 4, 160 00 Prague 6, Czech Republic, mar[email protected] Institute of Thermomechanics, Czech Academy of Sciences, Dolejˇ skova 5, 182 00 Prague 8, Czech Republic 2Institute of Thermomechanics, Czech Academy of Sciences, Dolejˇ skova 5, 182 00 Prague 8, Czech Republic, [email protected] 3Institute of Thermomechanics, Czech Academy of Sciences, Dolejˇ skova 5, 182 00 Prague 8, Czech Republic, [email protected] Keywords: CFD, Coherent structures, Dynamic mode decomposition, Higher-order dynamic mode decomposition. In many industrial applications, the dynamical behavior of the flow field is of great importance, as it determines the excitations acting on any solid bodies submerged in the liquid, the noise generated by the flow field, etc. As the computing power of modern computers increases, fully dynamical simulations are making their way to engineering practice, generating the need for new advanced validation methods tailored specifically for such simulations. Since fully dynamical simulations generate vast amounts of data, data decomposition methods must be employed. There exists a broad spectrum of data decomposition methods applicable to fluid flow problems; see, e.g., [1] and [2]. The use of Proper Orthogonal Decomposition (POD) for the validation of global dynamics has been proposed in [3]. This contribution aims to expand this idea by using HigherOrder Dynamic Mode Decomposition (hoDMD). HoDMD is one of the cutting-edge frequency-based methods, first proposed around 2016 by [4]. One of the advantages of all DMD-based methods is that they extract single-frequency modes, allowing a frequency-specific search to be employed. This in turn enables focusing the validation on the frequency band of interest, e.g. around eigenfrequencies of solid structures. HoDMD has been proved to be robust and reliable when accurate extraction of important data from a highly nonlinear underlying system is needed [5, 6]. In the case presented in this contribution, the simulation data is obtained from a digital twin of the 25 ×25 ×300 cm experimental test section. A circular cylinder, which spans the entire test section, is placed at the beginning of the test section perpendicular to the flow; see Figure 1. The inlet velocity is uin =5m/s, which results in a Reynolds number of Re = 4815. Inlet u= (uin; 0; 0)T n· ∇p= 0 k=kin ω=ωin Virtual walls Slip condition for u n· ∇p= 0 n· ∇k= 0 n· ∇ω= 0 Outlet n· ∇u= (0; 0; 0)T p= 0 n· ∇k= 0 n· ∇ω= 0 Walls u= (0; 0; 0)T n· ∇p= 0 Wall function for k Wall function for ω Cylinder u= (0; 0; 0)T n· ∇p= 0 k= 0 n· ∇ω= 0 Plane of interest Figure 1: Sketch of the computational domain with boundary conditions and part of the wind tunnel. Discrepancy between the domain and wind tunnel (virtual walls) in dark blue, cylinder in light blue, plane of interest in yellow. Image inspired by [3]. 8 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria The numerical solution is based on the finite-volume method and is computed using OpenFOAM. The mesh was prepared in such a manner that most of the region of interest is simulated with a resolution close to direct numerical simulation (DNS). The rest of the region of interest is well inside the large eddy simulation (LES) region. The use of DES approach allows us to set well-defined boundary conditions on the test section boundaries. The boundary conditions used are mostly standard; see Figure 1. For an exact case description, as well as any additional reference and case setup, the reader is referred to [3]. The experimental data used were measured in the plane of interest using time-resolved particle image velocimetry (TR-PIV) and are the same as in [3]. The numerical and experimental datasets contain 4000 snapshots of the flow field sampled at fs= 2 kHz and were pre-processed using the inherent symmetry of the case as in [3]. The main results are shown in Figure 2. The results compare the main antisymmetric modes computed from the experimental data and the OpenFOAM results. Figure 2 shows great qualitative agreement between experiment and computation, as expected, hinting at the possibility of using the hoDMD decomposition as a validation tool, provided it is combined with suitable validation metrics. Good agreement can also be seen between POD, hoDMD, and the theoretical vortex shedding frequency of fSt = 69.4 Hz. The POD frequencies were taken as the location of the highest peak in the power spectrum of the respective mode. (a) y x (b) fDMD = 69.0 Hz fPOD = 69.5 Hz fDMD = 71.1 Hz fPOD = 71.4 Hz negative positive Figure 2: Streamwise components of first hoDMD mode from (a) Experimental data and (b) OpenFOAM calculation Acknowledgments The authors acknowledge the financial support provided by the Ministry of Education, Youth, and Sports of the Czech Republic via the project No. CZ.02.01.01/00/23 020/0008501 (METEX), co-funded by the European Union. The work was financially supported by the institutional support RVO:61388998, by the Czech Science Foundation (GA 25-17815S), and by the grant project with No. TN02000069/001N of the Technology Agency of the Czech Republic. This work was supported by the Grant Agency of the Czech Technical University in Prague, grant No. SGS25/127/OHK2/3T/12. References [1] P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, and W. Schilders, Model Order Reduction. Berlin, Germany: de Gruyter, 2021, vol. 2. [2] K. Taira, S. L. Brunton, S. T. M. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: An overview,” AIAA Journal, vol. 55, no. 12, 2017. [3] T. Hlavat´ y, M. Isoz, M. Belda, V. Uruba, and P. Proch´ azka, “Is the proper orthogonal decomposition suitable to validate simulation of turbulent wake?” Journal of Wind Engineering and Industrial Aerodynamics, vol. 255, p. 105953, 2024. [4] S. Le Clainche and J. M. Vega, “Higher order dynamic mode decomposition,” SIAM Journal on Applied Dynamical Systems, vol. 16, no. 2, pp. 882–925, 2017. [5] E. Rodr´ ıguez-L´ opez, D. W. Carter, and B. Ganapathisubramani, “Dynamic mode decomposition-based reconstructions for fluid–structure interactions: An application to membrane wings,” Journal of Fluids and Structures, vol. 104, 2021. [6] A. Corrochano, G. D’Alessio, A. Parente, and S. L. Clainche, “Higher order dynamic mode decomposition to model reacting flows,” International Journal of Mechanical Sciences, vol. 249, p. 108219, 2023. 9 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NON-SPHERICAL DEM CONTACT MODEL CONSISTENT WITH HERTZIAN FORMULATION FOR LARGE-SCALE CFD-DEM OND ˇ REJ STUDEN´ IK1, MARTIN KOTOU ˇ Cˇ SOUREK2MARTIN ISOZ3, 1,3Czech Academy of Sciences, Institute of Thermomechanics, CZ 1,2Department of Chemical Engineering, University of Chemistry and Technology, Prague, CZ 1[email protected], [email protected], 3[email protected] Keywords: non-spherical particles, discrete element method (DEM), CFD-DEM Suspensions are prevalent in both nature and industry. They are present in diverse contexts, from riverbed transport and water seepage in gravel soils to industrial applications such as fluidized bed reactors and pharmaceutical microdosing [1]. Various approaches can be used to describe such systems, the most common being empirical, experiment-based correlations [2], which are standard in engineering practice. A more detailed alternative is the utilization of two-fluid models, where the particles are considered as the second fluid, and the fluid-”fluid” interaction is phenomenologically modeled. Lastly, a high-fidelity alternative is the combination of Computational Fluid Dynamics (CFD) for the description of flow, and Discrete Element Method (DEM) for the description of particles. The coupling of CFD-DEM is commonly done using one of three strategies: (i) one-way coupling, where only the fluid influence on the particle motion is considered; (ii) two-way coupling, which accounts for the mutual interactions between fluid and particles but neglects the particle-particle interactions; and (iii) four-way coupling, which additionally considers particle–particle interactions. Recently, we have been developing a four-way coupled CFD-DEM solver called openHFDIB-DEM [3], which is implemented as a custom extension for the OpenFOAM C++ framework. The software is detailed in [4]. It is designed as a CFD-DEM library for both pure DEM and CFD-DEM applications; the CFD-DEM library primarily focuses on four-way coupled particleresolved DNS simulations for coarse-grained slurry flows. In the present contribution, we will focus on pure DEM because the code was extended by a new formulation of a contact model for simulations of granular matter consisting of arbitrarily shaped polyhedral particles represented by triangulated surfaces (STL). In its original formulation, DEM was proposed only to account for spheres, with easily described geometry using only two parameters: center position and radius, while contact forces were calculated analytically using Hertzian contact theory [5], scaling them proportionally to the length of the particle overlap. However, since spheres are not prevalent, more realistic geometrical representations with varying levels of fidelity are required. The most widely used shape models can be divided into three types: (i) the multisphere model, (ii) superquadrics, and (iii) polyhedral-based particles [6]. Models (i) and (ii) are particularly suited for particles with smooth, curved surfaces, benefiting from the standard contact model with minor modifications. On the other hand, the polyhedral model is the most general and ideal for representing particles with sharp edges and flat faces [7]. On the other hand, this approach demands estimation of particle properties such as mass, center of mass, and moment of inertia, along with significant adjustments to contact detection and force calculation with respect to particle overlap characterization. The main novelty of this contribution is focused on the recent reformulation of the contact model and its alignment with the LIGGGHTS DEM solver [8]. As a demonstration, we provide simulation results for repose angle tests. We first analyze spherical particles, comparing the results of openHFDIB-DEM with those of LIGGGHTS, as shown in Figure 1. The comparison is then extended to non-convex polyhedral particles (Figure 2). These cases were chosen because of the significant role of tangential forces in the overall behavior of the system. The differences in behavior between spherical and non-spherical particles are apparent. Thus, a correct DEM implementation is of utmost importance for large CFD-DEM simulations of coarse-grain slurry flows. More details and results will be presented at the conference presentation. 16 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria (a) µ= 0 ⇝˜α= 0 o(b) µ= 0.3⇝˜α= 14.69 o(c) µ= 1.0⇝˜α= 24.78 o openHFDIB-DEM LIGGGHTS initial state x y Figure 1: Repose angle (˜α) study results, displayed as overlay in end time t= 2 s for the two tested solvers and three different friction coefficients (µ). The study contained 1700 particles. (a) µ= 0.1⇝˜α= 4.85 o(b) µ= 0.5⇝˜α= 9.70 o(c) µ= 0.7⇝˜α= 12.05 o initial state particle detail x y Figure 2: Results of repose angle (˜α) study, displayed in end time t= 2 s for the non-convex STL-based particles and three different friction coefficients (µ). The study contained 770 particles. Acknowledgments This research was co-funded by the European Union under the project Metamaterials for thermally stressed machine components (reg. no. CZ.02.01.01/00/23 020/0008501). The work was financially supported by the institutional support RVO:61388998, by the grant project with No. 25-17815S of the Czech Science Foundation of the Czech Republic and by the grant project with No. TN02000069/001N of the Technology Agency of the Czech Republic. References [1] H. Ma, L. Zhou, Z. Liu, M. Chen, X. Xia, and Y. Zhao, “A review of recent development for the cfd-dem investigations of non-spherical particles,” Powder Technology, vol. 412, p. 117972, 2022. [2] Z. Bi, Q. Sun, F. Jin, M. Zhang, and C. Zhang, “Temporal correlation of force and position in granular materials,” Particuology, vol. 10, no. 3, pp. 306–309, 2012. [3] Institute of Thermomechanics of the CAS, techMathGroup (M. Isoz, O. Studen´ ık, M. Kotouˇ cˇ Sourek, P. Koˇ c´ ı), “OpenHFDIB-DEM.” [Online]. Available: https://github.com/techMathGroup/openHFDIB-DEM [4] O. Studen´ ık, M. Isoz, M. Kotouˇ cˇ Sourek, and P. Koˇ c´ ı, “OpenHFDIB-DEM: An extension to OpenFOAM for CFD-DEM simulations with arbitrary particle shapes,” SoftwareX, vol. 27, 2024. [5] D. Antypov and J. A. Elliott, “On an analytical solution for the damped hertzian spring,” Europhysics Letters, vol. 94, 2010. [6] W. Zhong, A. Yu, X. Liu, Z. Tong, and H. Zhang, “DEM/CFD-DEM modelling of non-spherical particulate systems: Theoretical developments and applications,” Powder Technology, vol. 302, pp. 108–152, 2016. [7] J. Chen, “Understanding the discrete element method: Simulation of non-spherical particles for granular and multi-body systems,” Ph.D. dissertation, 2012. [8] C. Kloss, C. Goniva, A. Hager, S. Amberger, and S.Pirker, “Models, algorithms and validation for opensource DEM and CFD–DEM,” Progress in Computational Fluid Dynamics, an International Journal, vol. 12, pp. 140–152, 2012. 17 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Modeling High-Viscosity Spray Injection Using Volume-of-Fluid Method and Large Eddy Simulation Po-Han Chen, Alberto Ceschin, Junjun Guo, Hong G. Im Clean Energy Research Platform, King Abdullah University of Science and Technology, Thuwal, 23955-6900, Kingdom of Saudi Arabia, pohan.c[email protected] Keywords: High-viscosity spray injection, Spray Atomization, Geometric Volume-of-Fluid Method, Large Eddy Simulation, Adaptive Mesh Refinement, Dynamic Load Balancing Introduction Despite the rise of renewables, fossil fuels—particularly heavy residues like heavy fuel oil (HFO) and vacuum residue oil (VRO)—remain essential as light crude supplies decrease. In entrained-flow gasifiers, these high-density, cost-effective feedstocks are converted into syngas via drying, pyrolysis, and char conversion, with efficiency governed by reaction kinetics and the interfacial surface area of liquid sprays. This study models the spray interface evolution of a VRO-like liquid to quantify its dynamics and enable a more environmentally friendly energy conversion. Methods The Volume-of-Fluid (VOF) approach is a simulation technique that models the interface between two immiscible fluids [1]. It is particularly effective at capturing the breakup of a liquid into ligaments and droplets as VOF tracks the liquid–vapor interface by computing the liquid volume fraction, α, in each computational cell. The plic-RDF algorithm reconstructs the interface accurately [2], and the isoAdvector algorithm geometrically advects and preserves interface sharpness [3]. The governing equations for Large Eddy Simulation (LES) for compressible flows, namely conservation of mass, momentum and energy, along with the advection equation for the liquid volume fraction, are given by: ∂ρ ∂t +∂ ∂xj (ρ˜uj)=0 (1) ∂ ∂t (ρ˜uj) + ∂ ∂xj (ρ˜ui˜uj) = −∂p ∂xi +ρgi+fsurf,i +µ∂2˜ui ∂xj∂xj −∂τR ij ∂xj (2) ∂(ρT) ∂t +∂ ∂xjρ˜ujT=∂ ∂xjαT ∂T ∂xj−α cv,l +1−α cv,g ∂ ∂xj (˜ujp) + ∂ ∂t ρ˜ui˜ui 2+∂ ∂xj˜uj ρ˜ui˜ui 2 (3) ∂ ∂t(αρ) + ∂ ∂xj (αρ˜uj)=0 (4) where the ˜ ·denotes Favre-filtered quantities, and the ·indicates filtered properties. tis time, uis velocity, xis the spatial variable, pis pressure, Tis temperature, and gis gravitational acceleration. The specific heat capacity at constant volume is cv, and αTis thermal diffusivity. The subgrid-scale stress tensor, τR ij , is closed using the wall-adapting local eddy-viscosity (WALE) model. µand ρare the weighted-average dynamic viscosity and density, respectively, calculated using the volume fraction. The surface tension force, fsurf , is represented via the continuum surface force (CSF) model [4]. In deriving the simplified energy equation (Equation 3), we neglect terms of minor magnitude, specifically pressure jump at the interface, viscous dissipation, and surface tension work. The framework is derived from the solver of compressibleInterIsoFoam in OpenFOAM. To enhance robustness and efficiency over a wide range of Mach numbers [5], we modify the pressure equation to be Equation 5, by treating the gas density as temperature-dependent via the ideal-gas law and assuming the liquid phase to be incompressible. (1 −α)ψg ρg Dp Dt−1 T DT Dt+∂˜uj ∂xj = 0 (5) where ρgis the gas density, and ψgis the gas compressibility in OpenFOAM. 18 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Computational Setup and Results Atomization is achieved using an air-blast atomizer. The liquid is injected axially at an average velocity of 0.096 m/s, while the gas entered the chamber through eight inclined side pipes, reaching a peak speed of 260 m/s. Along with the large velocity difference, the high viscosity and density differences also pose challenges for this multiphase flow simulation. Using the setAlphaField utility in OpenFOAM, we initialize a pre-existing liquid column in the atomization zone with a velocity of 0.096 m/s to mitigate delays in spray formation caused by the low injection speed. The computational grids for this study are generated using blockMesh and snappyHexMesh tools in OpenFOAM. Adaptive mesh refinement (AMR) is applied to capture interface dynamics with a minimum cell size of approximately 62.5 microns. To balance the uneven workload induced by AMR, dynamic load balancing redistributes cells across processors for optimal performance. Figure 1 shows the spray evolution. High liquid viscosity produces elongated ligaments that subsequently fragment into droplets, and their lateral protrusion promotes a potential breakup into clusters of droplets. Figure 1: Time evolution of the air-blast spray. Liquid is injected axially through the central tube, while gas enters via eight surrounding pipes. The liquid–gas interface is colored by velocity magnitude, and the gray surface shows the Mach number 0.2 isosurface in the gas phase. Acknowledgments This work was sponsored by King Abdullah University of Science and Technology (KAUST). Computational resources were provided by the KAUST Supercomputing Laboratory (KSL). References [1] C. W. Hirt and B. D. Nichols, “Volume of fluid (vof) method for the dynamics of free boundaries,” Journal of computational physics, vol. 39, no. 1, pp. 201–225, 1981. [2] H. Scheufler and J. Roenby, “Accurate and efficient surface reconstruction from volume fraction data on general meshes,” Journal of computational physics, vol. 383, pp. 1–23, 2019. [3] J. Roenby, B. E. Larsen, H. Bredmose, and H. Jasak, “A new volume-of-fluid method in openfoam,” in Marine vi: Proceedings of the vi international conference on computational methods in marine engineering. CIMNE, 2017, pp. 266–277. [4] J. U. Brackbill, D. B. Kothe, and C. Zemach, “A continuum method for modeling surface tension,” Journal of computational physics, vol. 100, no. 2, pp. 335–354, 1992. [5] D. Fuster and S. Popinet, “An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension,” Journal of Computational Physics, vol. 374, pp. 752–768, 2018. 19 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria VERIFICATION OF A FINITE VOLUME SOLVER FOR ACTIVE AND PASSIVE CARDIAC MATERIAL BEHAVIOUR AARON MULLEN-HALES, PHILIP CARDIFF School of Mechanical and Materials Engineering, University College Dublin, Ireland, aar[email protected], philip.car[email protected] Keywords: Finite volume method, solids4foam, Cardiac mechanics, Hyperelasticity solids4foam [1] is an open-source computational toolbox for solving solid mechanics and fluid-solid interaction simulations in OpenFOAM. The study focuses on verifying solids4foam for active and passive cardiac material behaviour based on the benchmarks presented by Land et al. [2] and Aróstica et al. [4]. The current work uses a newly developed Jacobian-free Newton-Krylov, cell-centred finite volume solid solver implemented in solids4foam [5]. To the best of the authors’ knowledge, this is the first application of a finite volume solver to these benchmarks; in contrast, all results presented by Land et al. [2] and Aróstica et al. [4] came from finite element methods. The benchmark problems (problems 2 and 3 from Land et al. [2]) consist of an idealised ellipsoidal ventricle geometry (Figure 1) loaded with a combination of internal pressure and active material contraction. The material behaviour is defined by the Guccione anisotropic hyperelastic material law [3]. A utility was developed to initialise the fibre distribution for the active contraction case (Figure 2, left). The displacement magnitude for the passive inflation case (problem 2 of Land et al. [2]) is shown in Figure 1 (right), while predictions for four different mesh refinements for this case are shown in Figure 2 (right) and coincide with the predictions from other groups presented in the benchmark paper. In addition, benchmark problems from Aróstica et al. [4] are explored, specifically, cases A and B, which use different geometry, material, pressure and activation models compared to Land et al. [2]. The setFibreField utility, activation function as well as the problems described in Land et al. [2] and Aróstica et al. [4] can be found at http://github.com/solids4foam/cardiacFoam, where the feature-petsc-snes branch of solids4foam has been used. Figure 1: Left: Idealised ventricle: undeformed geometry (103 680 cells); Right: Idealised ventricle: deformed geometry 20 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 2: Left: Idealised ventricle: fibre field distribution; Right: Idealised ventricle: predictions for the four meshes from the present work overlaid on the round-robin predictions from Land et al. [2] Acknowledgements This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 101088740). Financial support is also acknowledged from I-Form, funded by Science Foundation Ireland (SFI), Grant Number 21/RC/10295_P2, co-funded under European Regional Development Fund and by I-Form industry partners, and from NexSys, funded by SFI Grant Number 21/SPP/3756. Additionally, the authors wish to acknowledge the DJEI/DES/SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support (www.ichec.ie), and part of this work has been carried out using the UCD ResearchIT Sonic cluster which was funded by UCD IT Services and the UCD Research Office. References [1] P. Cardiff, solids4foam: A toolbox for performing solid mechanics and fluid-solid interaction simulations in OpenFOAM. Journal of Open Source Software, 2024. [2] S. Land et al., Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proceedings of the Royal Society A, 2015. [3] J. M. Guccione, K. D. Costa, A. D. McCulloch, Finite element stress analysis of left ventricular mechanics in the beating dog heart. Journal of Biomechanics, 1995. [4] R. Aróstica et al., A software benchmark for cardiac elastodynamics. Computer Methods in Applied Mechanics and Engineering, 2025. [5] P. Cardiff et al., A Jacobian-free Newton-Krylov method for cell-centred finite volume solid mechanics. arXiv:2502.17217, 2025. 21 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria OPENFOAMGPT: AN AI ENGINEER JINGSEN FENG1, RAN XU2, XU CHU1,2 1Faculty of Environment, Science and Economy, University of Exeter, Exeter EX4 4QF, United Kingdom 2Cluster of Excellence SimTech, University of Stuttgart, Stuttgart, Germany, [email protected] Keywords: OpenFOAMGPT, LLM, agent We developed a Large Language Model (LLM)-based agent to advance Computational Fluid Dynamics (CFD) workflows1,2. CFD simulations often involve intricate configurations, including solver selection, initialand boundary condition setup, and iterative adjustments, which require substantial expertise and manual effort. The proposed LLMCFD framework integrates natural language processing capabilities with CFD tools, such as OpenFOAM, to assist users in navigating these complexities more efficiently. The LLM agent employs retrieval-augmented generation (RAG) on domain-specific datasets to provide accurate, context-aware guidance. Key functionalities include interactive solver configuration, and error diagnostics, allowing users to focus on analysis and decision-making. Reliability is a central focus of this work, addressing challenges such as potential hallucinations and uncertainty in LLM-generated outputs. The framework incorporates validations and benchmarks to ensure dependable performance. The integration of LLM technology into CFD workflows has the potential to enhance productivity by automating routine tasks, reducing human error, and improving accessibility for non-specialists. It also provides opportunities for interdisciplinary collaboration by bridging the gap between domain experts and computational tools. The presentation will highlight use cases demonstrating the agent’s capabilities, discuss methods to mitigate reliability risks, and outline future directions for applying LLMs in CFD research and industry. This work aims to contribute to the ongoing discourse on LLM applications in CFD, fostering innovation in simulationbased engineering and expanding the scope of fluid dynamics research. Figure 1: Structure of the OpenFOAMGPT 22 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria References [1] OpenFOAMGPT: A retrieval-augmented large language model (LLM) agent for OpenFOAM-based computational fluid dynamics, S Pandey, R Xu, W Wang, X Chu, Physics of Fluids 37, 035120 [2] A Status Quo Investigation of Large Language Models towards Cost-Effective CFD Automation with OpenFOAMGPT: ChatGPT vs. Qwen vs. Deepseek, W. Wang, R. Xu, J. Feng, Q. Zhang, X. Chu, https://doi.org/10.48550/arXiv.2504.02888 23 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria A TEMPORALLY SECOND-ORDER ACCURATE FINITE VOLUME DISCRETIZATION OF ONE-FIELD NAVIER-STOKES EQUATIONS WITH HIGH DENSITY RATIOS JUN LIU1, TOBIAS TOLLE1, DAVIDE ZUZIO2, JEAN-LUC ESTIVALEZES2, SANTIAGO M. DAMIAN3, DIETER BOTHE1, TOMISLAV MARIC1 1 Technical University of Darmstadt, Mathematics Department, Mathematical Modeling and Analysis, Darmstadt, Germany, [email protected], bo[email protected]dt.de, [email protected] 2 Universite de Toulouse, France, [email protected], jean-luc.estiva[email protected] 3 Centro de Investigaciones en Mecanica Computacional (CIMEC), UNL/CONICET, Argentina, santiagomarquez[email protected] Keywords: two-phase flows, high density ratios, one-field Navier-Stokes equations, unstructured finite-volume method Abstract We derive conditions for consistent volume, mass, and momentum conservation in the one-field formulation of NavierStokes equations for incompressible two-phase flows without phase change. The derived consistency conditions are independent of the interface tracking method or the equation discretization scheme. Aiming at two-phase flows in geometrically complex engineering systems, we discretize the model equations and the consistency conditions using the unstructured Finite Volume method and verify and validate our approach for two very different interface capturing methods, namely the Level Set / Front Tracking method and the flux-based Volume-of-Fluid method. Verification and validation demonstrate stability when simulating challenging two-phase flows with density ratios up to 106 and viscosity ratios up to 105. We finally discuss the requirements for consistent mass-flux approximation with second-order temporal accuracy in the context of the Crank-Nicolson temporal discretization scheme. Introduction Incompressible two-phase flows that involve fluids with very different densities (high density ratios) are ubiquitous in natural and technical processes, and the one-field formulation of Navier-Stokes Equations (one-field NSE) is widely used to model such flows. The one-fluid, integral form of the two-phase NSE enforces consistency between volume, mass, and momentum conservation, which, when not ensured in the discretization, leads to instabilities. The phase indicator function further uniquely defines the mass flux, which must be consistently used in the mass conservation and momentum conservation equation. These exact consistency requirements, derived from the integral form of the onefield NSE, are independent of the method used to model the fluid interface. Research into two-phase flow simulation methods for handling high density ratios is highly active; however, focusing primarily on the discrete level. In this work, we derive exact consistency requirements at the level of the mathematical model. These conditions must be tailored to the PDE discretization and fluid interface tracking methods, but they must be upheld. Since we develop numerical methods for simulating geometrically complex engineering multiphase flow systems, we discretize the one-field NSE with the consistency requirements using the unstructured Finite Volume method and validate and verify them for the unstructured Level Set / Front Tracking [1] method and plicRDF-isoAdvector, a flux-based geometrical Volume-ofFluid method [2]. Method The solution domain Ω is split into two sub-domains filled by incompressible fluids of constant densities 𝜌±, Ω≔ Ω(𝑡)∪Ω(𝑡)\Σ(𝑡), separated by the sharp fluid interface Σ(𝑡):=𝜕Ω(𝑡) (cf. Fig. 1). The phase indicator identifies Ω±(𝑡) as 𝜒(𝒙,𝑡)=1, if 𝒙∈Ω(𝑡), and 𝜒(𝒙,𝑡)=0, if 𝒙∈Ω(𝑡), resulting in the one-field density 𝜌(𝒙,𝑡)= (𝜌−𝜌)𝜒(𝒙,𝑡)+𝜌, with constant phase-specific densities 𝜌±. Integrating 𝜌(𝒙,𝑡) over a fixed control volume Ω (cf. Fig. 1), and defining the volume fraction 𝛼(𝑡)≔ ||∫𝜒(𝒙,𝑡)𝑑𝑉 , results exactly in 𝜌(𝑡)=(𝜌−𝜌)𝛼(𝑡)+ 𝜌, if the integral ∫𝜒(𝒙,𝑡)𝑑𝑉  is exactly evaluated. This density update is used by all methods that discretize one-field NS, even those that do not actually solve a transport equation for the volume fraction or the phase indicator and, instead, replace 𝛼(𝑡) by an approximate “marker field” modeled by the fluid interface Σ(𝑡), e.g., from the Level Set function or the Front in Front Tracking methods. The relation 𝜌(𝑡)=(𝜌−𝜌)𝛼(𝑡)+𝜌, although exact, hides the following additional exact consistency requirement between the mass and volume transport, namely, with shorthand notation 𝜙≔𝜙(𝑡), 24 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria   𝜕𝜌   𝑑𝑉=𝜌− 𝜌=−(𝜌−𝜌) |Ω|  𝜒𝒗⋅ 𝒏 𝑑𝑆 𝑑𝑡   −𝜌 |Ω|  𝒗 ⋅ 𝒏 𝑑𝑆    𝑑𝑡, (1) with ∫𝒗 ⋅𝒏 𝑑𝑆 =∫(∇ ⋅ 𝒗) 𝑑𝑆 =0 for incompressible two-phase flows. In other words, the control-volume density average is updated by a phase-specific volumetric flux 𝜒𝒗⋅𝒏, scaled with the density difference (𝜌−𝜌). If the method used to track the fluid interface Σ(𝑡) does not use a phase-specific volumetric flux to track Σ(𝑡), the equivalence between 𝜌(𝑡)=(𝜌−𝜌)𝛼(𝑡)+𝜌, when 𝛼(t) is a “marker field,” and equation (1) is not ensured, and this inconsistency leads to artificial acceleration and deceleration, i.e., parasitic currents, which are catastrophic in some cases. The solution for methods that do not rely on volume fractions (phase indicator) transport is to solve an auxiliary density equation to ensure equivalence defined by (1) when solving onefield NSE. Solving the auxiliary density equation has been done for a variety of methods in the literature on the discrete level; we introduce the mass conservation equation on the level of the mathematical model. Methods that do not transport the phase indicator function require an exact reconstruction of the mass flux 𝜌𝒗(𝒙,𝑡)=(𝜌−𝜌)𝜒𝒗(𝒙,𝑡)+𝜌𝒗(𝒙,𝑡), by reconstructing exactly 𝜒(𝒙,𝑡) from Σ(𝑡). Discretely, we approximate 𝜒(𝒙,𝑡)≈𝜒(𝒙,𝑡) from the approximation of Σ(𝑡). For the unstructured Level Set / Front Tracking method [1], we approximate 𝜒(𝒙,𝑡) from signed distances. For the plicRDF-isoAdvector flux-based geometrical Volume-of-Fluid method [2] we demonstrate the equivalence of solving an auxiliary mass conservation equation and discretizing the one-field NSE with the Euler and upwind schemes when (𝜌−𝜌)𝜒𝒗(𝒙,𝑡) in the mass flux 𝜌𝒗(𝒙,𝑡) is approximated by scaling the phase-specific fluxed volume. We analyze and apply the consistency requirements in the context of the Crank-Nicolson second-order temporal discretization scheme. Finally, on the discrete level, once a consistent mass flux has been approximated, the same mass flux must be re-used when discretizing the momentum equation. Results and Discussion Challenging verification and validation test cases involving droplets and jets with density ratios up to 106 and viscosity ratios up to 105 confirm our analysis. Verification involving droplets translated with constant velocity results in solutions that are accurate up to the linearsolver tolerance and maintain sphericity (Fig. 2), while droplets transported in ambient fluid remain stable. An exceptionally challenging mixing layer case in an inviscid setting remains stable. Finally, a coarsely resolved jet in crossflow remains stable, and its trajectory is successfully validated with experimental data. In all cases, when our suggested combination of consistent schemes, or, equivalently, an auxiliary mass conservation equation, is not used, errors increase significantly in the best case, and simulations catastrophically fail in the worst-case scenario. Conclusion We demonstrate how to ensure the consistency of volume, mass, and momentum conservation prescribed by the one-field formulation of Navier-Stokes Equations for incompressible two-phase flows, for methods that do not transport a phase indicator function, and how to avoid hidden sources of inconsistency in Volume-of-Fluid methods. References [1] Liu, J., Tolle, T., Bothe, D., & Marić, T. An unstructured finite-volume level set/front tracking method for two-phase flows with large density ratios. Journal of Computational Physics, 493, 112426, (2023). [2] Liu, J., Tolle, T., Zuzio, D., Estivalèzes, J. L., Damian, S. M., & Marić, T. Inconsistencies in unstructured geometric volume-of-fluid methods for two-phase flows with high density ratios. Computers & Fluids, 106375, (2024) Figure 1 : Two - phase flow domain. Figure 2: Top left and top right images show catastrophic failures of inconsistent methods for a translating droplet. Lower left image shows artificial deformation of the droplet caused by inconsistent discretization. Lower right image shows a more accurate droplet motion given by the consistent method. 25 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria An accurate near-field aeroacoustic model for low Mach number flows Benedetto Di Paolo1, Eugene de Villiers2, Jakub Knir3, Paolo Geremia4 1ENGYS S.r.l., b[email protected] 2ENGYS Ltd., e.de[email protected] 3ENGYS Ltd., [email protected] 4ENGYS S.r.l., p.ger[email protected] Keywords: Aeroacoustics, CAA, PCWE, CFD Aeroacoustic modelling plays a vital role in modern engineering, increasingly driven by the demand for quieter and more sustainable technologies. This need extends far beyond the traditional aerospace and automotive sectors, affecting electric vehicles, renewable energy systems like wind turbines, railway infrastructure, fans, HVAC systems, turbomachinery, consumer electronics, and even marine propulsion. Regulatory frameworks such as ISO 3744, ISO 9613, IEC 61400-11 for wind turbines, ISO 5136 for ducted fans, EN ISO 3095 for railways, and ICAO Annex 16 for aviation, as well as RTCA DO-160 and EU Directive 2000/14/EC, increasingly require predictive noise control at the design stage. The growing regulatory and environmental pressure to reduce noise emissions across industries has placed Computational Aeroacoustics (CAA) at the heart of modern product development. As noise limits tighten and sustainability goals become more ambitious, accurate and efficient predictive tools are essential. CAA enables high-fidelity simulations of noise generation and propagation in complex engineering systems, providing a deep understanding of flow-induced acoustic phenomena without the need for extensive and expensive experimental campaigns. This accelerates the design process, reduces development costs, and supports early-stage decision-making. Crucially, CAA allows engineers to isolate dominant noise sources, explore mitigation strategies, and ensure compliance with evolving standards, all within a virtual, controllable environment. In this work, the Perturbed Convective Wave Equation (PCWE) is implemented in HELYX [1] to study low Mach number flows. CAA model based on PCWE has many advantages over more traditional CAA methods such as acoustic analogies (e.g. Ffowcs Williams and Hawkings) or Direct Noise Computations. The Perturbed Convective Wave Equation (PCWE) offers an efficient and versatile approach to aeroacoustic modelling. It directly computes the acoustic pressure field using a simple, easy-to-evaluate scalar source term, which reduces computational cost and memory requirements. Acoustic wave reflection and refraction can be analysed. The scalar source term is easily visualisable and can be derived directly from any incompressible unsteady flow simulation, relying only on the incompressible pressure and mean velocity fields, making it broadly applicable and easy to integrate into existing CFD frameworks [2]. The work also explores the pivotal role of advanced numerical methods in aeroacoustic simulations, with a particular focus on higher-order time integration schemes and sub-cycling strategies that enhance both accuracy and efficiency. It demonstrates how the modular USF (Unified Solver Framework) technology of HELYX enables seamless coupling between flow and acoustic fields, including the integration of synthetic noise generation. Additionally, the study presents validation and verification (V&V) efforts that support the reliability and robustness of the approach. References [1] ENGYS, HELYX - Open-source CFD Software for Enterprise. [2] S. Schoder, Aeroacoustics Theory and methods for analyzing flow-induced sound generation of technical and biological applications, 2024. 32 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria ASYNCHRONOUS EULER-LAGRANGE COUPLING METHOD SERGEY LESNIK1, SILVIO SCHMALFUß2, HENRIK RUSCHE3 1Wikki GmbH, [email protected] 2Leibniz-Institute for Tropospheric Research, Atmospheric Microphysics Departement, [email protected] 3Wikki GmbH, [email protected] Keywords: Euler-Lagrange simulations, HPC, Parallel computing, GPUs Motivation Clouds play pivotal role in weather and climate, e.g. by influencing precipitation and the Earth’s radiative budget. Understanding cloud microphysical and dynamical processes, their interactions, and their coupling with aerosol particles is a challenging task due to the turbulent and vast-scale nature of clouds. Cloud chambers offer a solution by creating controlled supersaturated environments for repeatable experiments [1]. While cloud chambers are valuable for studying cloud microphysics, they are often complemented by simulations, especially for properties without suitable measurement techniques or those challenging to measure across the entire domain. Such simulations typically employ the Euler-Lagrange (E-L) method, where the continuous phase, such as air, is represented within the Eulerian frame of reference and the dispersed phase, such as discrete droplets, within the Lagrangian frame of reference. In the case of clouds, the two phases significantly affect each other. This is depicted in the model by the so-called two-way coupling that considerably complicates simulations. A representative result of such a simulation for the Pi-Chamber [2] is shown in Figure 1. It shows particle distribution and size as well as the relative humidity in a diagonal plane. Figure 1: A snapshot of a simulation of the Pi-Chamber showing the relative humidity in a diagonal cross section and a random selection of droplets and their diameters. While Lagrangian particle tracking is well suited for parallel execution because each particle can be tracked independently, achieving efficient parallel scaling for the two-way coupled E-L method is a well-known challenge [3]. Consequently, tracking all particles in real-world situations, such as an 81 m3 cloud chamber with a particle/droplet concentration of 1000 cm−3 (equalling a total of 81·109 particles) as currently being designed, is not viable, even with the most powerful computers. The presented approach is pushing these boundaries. Algorithm Design The calculations within the Eulerian and Lagrangian frames of reference are executed asynchronously on CPUs and GPUs, respectively. In the novel algorithm, transfer of the Eulerian data that is necessary to perform the Lagrangian particle tracking is initiated as soon as it becomes available and the transfer of the Eulerian source terms calculated within the Lagrangian particle tracking is initiated as soon as the particle tracking is finished. The conservation of exchanged properties such as total liquid mass or total linear momentum are satisfied in a time-average sense. This can be achieved by predicting the Eulerian source terms based on values from previous time steps and correcting them later such that the integral matches what has been calculated by the Lagrangian tracking. Furthermore, the particles are 33 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria organised in chunks in a cache friendly way, and a bounding box around each chunk is maintained to facilitate dynamic load balancing across GPUs and to optimise data transfers between GPUs and CPUs. Results The algorithm’s scalability and efficiency are demonstrated by simulating a 81 m³ cloud chamber containing up to 256 billion particles. Results from this simulation are presented, including information on the parallel efficiency and comparison to reference cases simulated with conventional CPU-based particle tracking in OpenFOAM. Figure 2 shows strong scaling for the cloud chamber discretised into 24 million cells and containing 8 billion droplets. For every 20 cores there is on dedicated Nvidia H100 GPU. The calculations belonging to the tracking algorithm show an almost ideal strong scaling. The overall performance scales approximately up to 512 cores, which is due to scaling limitations of OpenFOAM. Figure 2: Strong scaling results showing execution time per simulated time step on the left y-axis and time-to-solution for 30'000 iterations on the right. Acknowledgements This project has received funding from the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 101118139. The JU receives support from the European Union’s Horizon Europe Programme. References [1] Niedermeier, D. et al. Characterization and first results from LACIS-T: a moist-air wind tunnel to study aerosol–cloud–turbulence interactions, Atmos. Meas. Tech. 13(4), 2015–2033 (2020). [2] Chang, K. et al. A laboratory facility to study gas–aerosol–cloud interactions in a turbulent environment: The π chamber., Bull. Am. Meteorol. Soc. 97(12), 2343–2358 (2016). [3] Bonnier, F. Algorithmes parallèles pour le suivi de particules, Doctoral dissertation, Université ParisSaclay (ComUE) (2018). 34 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria MODELING CAKE FILTRATION VIA POROUS MEDIA APPROACH JOSEF TAUSENDSCHÖN1 1Research Center Pharmaceutical Engineering GmbH, [email protected] Keywords: Porous media, Cake filtration, Modeling Cake filtration is a widespread form of solid-liquid separation in industrial applications [1]. The flow behavior in cake filtration is commonly described using Computational Fluid Dynamics (CFD), where the cake is often modelled as porous media [2]. The principle of cake filtration, in the present context, is to remove solid particles from a liquid through a porous medium that permits fluid flow while retaining solids. As a slurry approaches the medium, particles accumulate on its surface, progressively forming a filter cake. The thickness of this cake increases over time, making cake formation and growth key dynamics in filtration. Further, cake height and structure might spatially differ. Despite their significance, these dynamic behaviors of the filter cake are often overlooked in existing numerical models. Out-ofthe-box available models for porous zones have fixed static coefficients describing the pressure loss, e.g., the oftenapplied Darcy-Forchheimer model, and are not able to adequately model the dynamics of the process. To acknowledge the spatial and temporal change of cake resistance, a new porous zone model is presented. The cake height (assuming all solids are retained on the filter, no liquid remains the cake, and a constant porosity) is defined by: ℎ𝑐𝑎𝑘𝑒 =𝑐�`𝑢�`�] (1−𝑐�`𝑢�`�])(1−𝜀) 𝑉 𝐴 (1) where �I is the filtrate volume, 𝐴 the filter area, 𝑐�`𝑢�`�] the solid volume fraction in the suspension and 𝜀 the cake porosity. The pressure gradient in the porous zone is then: ∇𝑝=(𝑅𝑚+𝑅𝑐𝑎𝑘𝑒) ℎ�]𝑜𝑟𝑜𝑢�` 𝜇 𝑈 (2) Where 𝑅𝑚 is the filter media resistance, 𝜇 the kinematic viscosity of the liquid, ℎ�]𝑜𝑟𝑜𝑢�` the height of the porous zone and 𝑅𝑐𝑎𝑘𝑒 the cake resistance. The cake resistance is based on constant specific cake resistance 𝛼𝐻, is: 𝑅𝑐𝑎𝑘𝑒 =𝛼𝐻 ℎ𝑐𝑎𝑘𝑒 (3) The implementation follows the available Darcy-Forchheimer model but is limited to laminar flows, and subsequently, the quadratic contribution is removed. The resulting pressure gradient of each cell is added then as source term in the momentum equation. The presented model is validated with experimental data from two classic dead-end cake filtration experiments [3], [4]. In constant pressure difference experiments, the resulting filtrate volume and resulting cake height were compared to the simulation. Given the cake porosity and the specific cake resistance, the model can accurately simulate the collected filtrate volume over time as well as the cake height. Acknowledgements The authors thank all those involved in the organisation of this OpenFOAM Workshop and to all the contributors that will enrich this event. The Research Center Pharmaceutical Engineering (RCPE) is funded within the framework of COMET – Competence Centers for Excellent Technologies by BMK, BMAW, Land Steiermark and SFG. The COMET program is managed by the FFG. Special thanks to Mathias Glatz from Takeda Manufacturing Austria AG and Dominik Duckgeischel from STRASSBURGER Filter GmbH for the collaboration and expertise, which greatly enriched the development. References [1] H. Anlauf, Wet Cake Filtration: Fundamentals, Equipment, and Strategies. WILEY-VCH, 2019. [2] H. Neunzert and D. Prätzel-Wolters, Currents in industrial mathematics: From concepts to research to education. Springer Berlin Heidelberg, 2015. doi: 10.1007/978-3-662-48258-2. [3] F. M. Mahdi and R. G. Holdich, “Laboratory cake filtration testing using constant rate,” Chemical Engineering Research and Design, vol. 91, no. 6, pp. 1145–1154, Jun. 2013, doi: 10.1016/j.cherd.2012.11.012. 35 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria [4] T. Buchwald, “Nonlinear Parameter Estimation of Experimental Cake Filtration Data,” Technische Universität Bergakademie Freiberg, 2021. 36 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NUMERICAL STUDY OF A BLUFF BODY STABILIZED LEAN PREMIXED HYDROGEN FLAME USING ARTIFICIALLY THICKENED FLAME MODEL IN OPENFOAM RICCARDO VADI1, ANTONIO ANDREINI 1HTC Group, DIEF Department of Industrial Engineering, University of Florence, Italy [email protected] Keywords: Reacting flow, Combustion, Hydrogen, Thickened Flame Model, LES, OpenFOAM In recent years, hydrogen has garnered significant attention due to its carbon-free nature, making it an appealing candidate for sustainable energy solutions. Its unique properties as an energy vector, such as high energy density and versatility across various applications, have further highlighted its potential in the transition to cleaner energy systems. It is important to characterise numerical models to expand knowledge in hydrogen combustion. In this study, the Thickened Flame Model [1] is implemented in the opensource software OpenFOAM v8. New solver is therefore developed from reactingFoam. Thermal diffusion effects will also be investigated by integrating the Soret effect into the governing equations. A case study of a single sector hydrogen burner is selected. The test case will show the characteristics of hydrogen in terms of diffusion during the reactive process. Analysis of the atmospheric test bench operating with a lean, perfectly pre-mixed, hydrogen-air flame stabilised on a conical bluff body [2] will be carried out. LES simulations are conducted in which the impact of modelling the interaction between the flame front and turbulence through the definition of efficiency function is analysed. In addition to the one proposed by Colin [1], the efficiency function presented by Charlette [2] is implemented. The simulation is conducted in parallel using the DLBFoam library [3] to balance the load in the reaction rate computation. An extension to the original library is also proposed for the reference zonal mapping model to speed up the chemistry. In reference mapping, the mixture fraction is replaced by the progress variable and the flame sensor to map a reference solution to reference cells that satisfy a conditions. The investigation is done by the qualitative comparison of the three-dimensional flame visualisation with the experimental test case. The quantitative one is conducted by comparing detailed Particle Image Velocimetry (PIV)/ OH-Planar LaserInduced Fluorimetry (OH-PLIF) measurements, which allow the characterisation of the flame behaviour. The comparisons will make it possible to identify the optimal flame representation models among those selected. References [1] Colin, O., et al. "A thickened flame model for large eddy simulations of turbulent premixed combustion." Physics of fluids 12.7 (2000): 1843-1863. [2] Charlette, Fabrice, Charles Meneveau, and Denis Veynante. "A power-law flame wrinkling model for LES of premixed turbulent combustion Part I: non-dynamic formulation and initial tests." Combustion and Flame 131.1-2 (2002): 159-180. [3] Yahou, Tarik, James R. Dawson, and Thierry Schuller. "Impact of chamber back pressure on the ignition dynamics of hydrogen enriched premixed flames." Proceedings of the Combustion Institute 39.4 (2023): 46414650. [4] Tekgül, Bulut, et al. "DLBFoam: An open-source dynamic load balancing model for fast reacting flow simulations in OpenFOAM." Computer Physics Communications 267 (2021): 108073. 37 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Simulation based Optimization of Wetting Phenomena Clara Bernklau1, Stefan Ulbrich2 1Clara Bernklau TU Darmstadt , [email protected] 2Stefan Ulbrich TU Darmstadt, ulbric[email protected] Keywords: Wetting Phenomena, Optimization, Multiphase flow, Doctor Blading We consider a wetting process that is motivated by gravure printing, a process where a gravure cylinder is rotated, while an ink reservoir is used to wet the surface with ink and a doctor blade is then pulled over the surface by the rotation of the cylinder to remove the excess fluid. This printing process is described by the two-phase Navier-Stokes equation with surface tension, see for example [1], and a dynamic contact line model. The optimization framework and doctor blade test case were developed by Diehl in [2]. The optimization algorithm is composed of an outer initialization and optimization framework that is implemented in Matlab. The inner simulation part is performed with OpenFOAM, where a sensitivity approach or an adjoint approach can be chosen for the calculation of the derivatives. We use the multiphase solver interFoam of OpenFOAM to compute the state and have extended it to compute the derivative of the objective function. This framework is schematically depicted in Figure 1. In a first step, discretization methods for the VOF-based sensitivity equations were derived and implemented in interSensFoam by Diehl in [2] following an optimize-then-discretize approach. Here, the terms containing sensitivities of the face fluxes must be discretized by suitably selected schemes, since for example the upwind direction has to be determined by the original face flux. Furthermore, particular attention is required for the linearisation of the surface tension term which is used by interFoam. In this approach the sensitivity δα is not corrected by the MULES algorithm, but rather computed with an upwind scheme and the sensitivities of the velocity and the pressure are calculated by the PISO algorithm in the same way as the state variables, where the additional terms are treated as source terms. The development of an adjoint solver based on interSensFoam is currently in progress. Figure 1: Coupling of Matlab and OpenFOAM simulations. The geometrical setup consists of a domain with a doctor blade and homogeneous Dirichlet boundary conditions on the upper and left boundaries for the velocity uwhile on the lower wall a velocity ucyl is applied. The fluid flows out on the right border. In the simulation, two initial conditions of the fluid are considered illustrated in Figure 2. An example of an optimization problem is the effect of the viscosity of the printing fluid on the film formation occurring behind the doctor blade, where the sensitivities are calculated using the interSensFoam solver. Another optimization goal is the reduction of air bubbles inside the printing fluid to avoid printing failures. Therefore, the minimization of the vorticity inside the ink reservoir is of interest. This can be influenced by changing the inclination angle of the doctor blade. 38 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 2: On the left side are the velocity boundary conditions of the doctor blading test case. In the middle is scenario 1 and on the right side scenario 2. References [1] J. Pr¨ uss and G. Simonett, “Simulation based optimization of multiphase flow in the context of wetting phenomena,” Interfaces Free Boundaries, pp. 311–345, 2010. [2] E. Diehl, “Simulation based optimization of multiphase flow in the context of wetting phenomena,” Ph.D. dissertation, Technische Universit¨ at Darmstadt, 2023. 39 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Verification and Validation of OpenFOAM for High-Speed Aerodynamics Viola Rossano1, Tony Di Fabbio2, Eric Segalerba2, Joel Enrique Rivas Guerrero2 1Leonardo Labs, Leonardo S.p.A., C.so Castelfidardo 22, 10138 Torino, Italy 2Leonardo Labs, Leonardo S.p.A.,Via Pieragostini 80, 16151 Genova, Italy Abstract Accurate prediction of sonic booms is essential for advancing the design of future supersonic aircraft with minimized boom impact. This study focuses on validating the performance of OpenFOAM in predicting nearfield pressure signatures for the test cases from the Second AIAA Sonic Boom Prediction Workshop. The test cases included an axisymmetric equivalent area body, a JAXA wing-body configuration, and a NASA low-boom supersonic configuration with flow-through nacelles. All simulations were conducted at a free-stream Mach number of 1.6, zero degrees angle of attack, and a Reynolds number of 5.7 million per meter. The simulations were performed using viscous-mixed grids provided by the workshop committee. A comparison of the results demonstrates excellent agreement between the solutions obtained with OpenFOAM and Ansys, as well as with the results from other participants. The locations of shocks and expansions were found to be consistent across all test cases. Keywords: Applied Aerodynamics, Compressible Flows, Supersonic Flows, Sonic Boom, HiSA 1 Introduction The era of commercial supersonic civil aviation effectively ended with the final flight of the Concorde in 2003. In many countries, including under regulations such as FAR 91.817, enacted in the 1960s, the operation of civil supersonic aircraft over land is prohibited unless it can be assured that the sonic boom will not reach the ground [1]. Recently, there has been a global surge in research focused on the next generation of supersonic aircraft, with particular attention given to sonic boom regulations. Accurate prediction and mitigation of sonic boom signatures are crucial for the development of supersonic transport aircraft and the possibility of unrestricted overland supersonic flight [2]. Research has demonstrated a notable market potential for supersonic transport aircraft [3, 4]; however, the time savings provided by high-boom airliners on a significant number of routes are not substantial enough to compensate for the decreased efficiency caused by the increased drag during supersonic cruise. A low-boom airliner would enable additional flight routes and create new market opportunities, but the reduction of sonic booms must be achieved to a level deemed acceptable by the general public [5]. Recent advancements in numerical methods have significantly enhanced the design of low-boom supersonic transport aircraft and the prediction of their sonic boom signatures, making it a focal point of ongoing research [6]. In the present study, the High Speed Aerodynamic (HiSA) solver [7], was used to compute flow fields for test cases presented in the Second AIAA Sonic Boom Prediction Workshop (SBPW2). The objective was to evaluate OpenFOAM-based software’s effectiveness as a tool for predicting nearfield pressure signatures on industry-relevant geometries and to identify areas for further development to improve the accuracy of sonic boom predictions. Simulations were performed with both HiSA and Ansys Fluent v2024R2, and the results were compared to those from other participants of SBPW2. The structure of this paper is as follows. Section 2 provides a description of the SBPW2 test cases and the simulation setup. Section 3 presents the results. Finally, Section 4 offers some concluding remarks and outlook. 2 Test Cases The Figure 1 shows the three test cases analyzed in the present study, which consist of an axisymmetric equivalent area body (AXIE), a JAXA Wing Body (JWB), and a NASA low-boom supersonic configuration with a flow-through nacelle (C25F). The AXIE body was designed using Cart3D to match the near-field pressure signature of the inviscid NASA Concept 25D with a flow-through nacelle at 3 body lengths on the centerline, as referenced in [8]. JAXA provided the second test configuration [9], JWB, based on an inverse design to recover the C25F equivalent area. The JWB model has a length of 38.7 m, a reference area of 65.6 m2, and a design angle of attack of 2.3◦. The third configuration, the NASA C25F with a 40 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria flow-through nacelle, is a low-boom conceptual vehicle designed at NASA Langley [10]. It includes a 3.37◦nose rotation to achieve the design angle of attack, with a length of 38.7 m and a reference area of 37.16 m2. All test cases were conducted at Mach 1.6 and an altitude of 15,760 m. Figure 1: Geometries of the second AIAA Sonic Boom Prediction Workshop CFD Solver Settings and Numerical Approach The details of the numerical approach and the computational grids used to compute the flow field for the SBPW2 test cases are presented. The simulations were performed using SBPW2 viscous-mixed grids, scaled to full-scale meters, designed for a free-stream Mach number of 1.6 and an angle of attack of zero degrees. The CFD simulations were carried out using Ansys Fluent, version 2024R2. A cell-centered finite volume solver was employed, utilizing the fully implicit density-based solution method with the Roe-FDS flux vector splitting scheme for convective fluxes. The Green-Gauss node-based method was used to compute secondary diffusion terms, velocity derivatives, and scalars at the cell faces. Additionally, second-order upwind spatial discretization schemes were applied to the flow equations (pressure, momentum, and energy). Finally, settings were configured to enable High-Speed Numerics (HSN) to stabilize and accelerate solution convergence. The HiSA flow solver, implemented in OpenFOAM [7], was used for the simulations. This solver employs the AUSM+up upwind scheme for calculating face fluxes [11], allowing for the simulation of unsteady flows with significantly higher Courant numbers compared to explicit solvers. Additionally, it has been reported that implicit solvers are more efficient and accurate than explicit solvers at high Mach numbers [12]. The governing equations were solved using the Generalized Minimal Residual (GMRES) algorithm [11], with Lower-Upper Symmetric Gauss–Seidel (LUSGS) preconditioning [13]. The implementation of characteristic boundary conditions for non-reflective farfield and reflective walls and dual time stepping method in the HiSA solver plays a critical role in accurately simulating and predicting sonic booms. In sonic boom simulations, non-reflective farfield boundary conditions ensure that shock waves generated by supersonic flight exit the computational domain without reflecting back. This is crucial for accurate predictions, as reflected waves could distort the simulation, especially in the farfield where the sonic boom is perceived. Reflective wall boundary conditions are used when simulating scenarios like wind tunnel tests or low-altitude flights, where shock waves interact with the ground or other surfaces. These conditions allow the solver to model the reflection of shock waves from the ground, a key factor in determining the sonic boom’s characteristics. Properly handling these reflective interactions ensures accurate predictions, particularly for ground-level boom impacts. Dual time stepping is a method used to enhance the accuracy and efficiency of simulations, especially for transient phenomena such as shock waves in sonic boom predictions. It involves two key components: physical time stepping, which advances the solution in real-world time intervals, and pseudo time stepping, which introduces an artificial time parameter to iteratively update the solution within each physical timestep, speeding up convergence. Additionally, local time stepping is used to adapt the timestep size across different regions of the computational domain. Smaller timesteps are applied to regions with rapid changes (e.g., shock waves), while larger timesteps are used for less dynamic areas, improving computational efficiency. In general, dual time stepping accelerates solution convergence by iterating in pseudo-time while advancing in physical time. Pseudo time stepping stabilizes the solution and refines shock interactions. Local time stepping optimizes timestep usage, ensuring critical regions are resolved with high accuracy. Together, these methods ensure efficient and accurate simulations, particularly for capturing transient effects like sonic booms. 41 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 13: Contour of ∆p p∞: ANSYS (left), OpenFOAM (right) Figure 14: Contour of ∂ ∂x ∆p p∞: ANSYS (left), OpenFOAM (right) 48 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 15: Numerical Schlieren visualization: ANSYS (left), OpenFOAM (right) Figure 16: Nearfield Normalized Pressure Signatures at H/L = 1. Reference data by [15]. 4 Conclusions This study presents the near-field sonic boom signatures for three benchmark test cases: an axisymmetric equivalent area body, a JAXA Wing Body, and a NASA low-boom supersonic configuration with flow-through nacelles. Simulations for all three cases were conducted using the grids provided for the SBPW2. For all cases, OpenFOAM performs well in predicting near-field sonic boom pressure signatures for the benchmark test cases considered. Minor discrepancies are observed in the aft part of the signatures that can be attributed to the complex flow region at the aft part of the C25F configuration. This study lays the foundation for applying the HiSA solver to more complex configurations and challenging test cases, such as transonic conditions and delta wings, including examples like the simple delta wing VFE-2 [16, 17]. Future work could explore multiple wing configurations and other intricate geometries to further understand aerodynamic phenomena, including shock formation, leading-edge vortices, and related flow dynamics [18, 19]. 49 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria References [1] Federal Aviation Administration, Civil Aircraft Sonic Boom, 2017, vol. Title 14. [2] T. R¨ otger, C. Eyers, and R. Fusaro, “A review of the current regulatory framework for supersonic civil aircraft: Noise and emissions regulations,” Aerospace, vol. 11, no. 1, p. 19, 2023. [3] B. Liebhardt and K. L¨ utjens, “An analysis of the market environment for supersonic business jets,” in 60th Deutscher Luftund Raumfahrtkongress, 2011, pp. 617–622, regulations Paper 241457. [4] B. Liebhardt, K. Luetjens, and V. Gollnick, “Estimation of the market potential for supersonic airliners via analysis of the global premium ticket market,” in 11th AIAA Aviation Technology, Integration, and Operations (ATIO) Conference, September 2011, aIAA Paper 2011-6806. [5] R. Seebass and B. Argrow, “Sonic boom minimization revisited,” in 2nd AIAA Theoretical Fluid Mechanics Meeting, June 1998, aIAA Paper 1998-2956. [6] A. Norman, V. Viti, K. MacLean, and V. Chitta, “Improved cfd methodology for compressible and hypersonic flows using a hessian-based adaption criteria,” in AIAA SCITECH 2022 Forum, 2022, p. 0582. [7] J. Heyns and O. Oxtoby, “Modelling high-speed viscous flow in openfoam®,” 9th South African Conference on Computational and Applied Mechanics, SACAM 2014, 01 2014. [8] M. A. Park and M. Nemec, “Nearfield summary and statistical analysis of the second aiaa sonic boom prediction workshop,” Journal of Aircraft, vol. 56, no. 3, pp. 851–875, 2019. [9] A. Ueno, M. Kanamori, and Y. Makino, “Multi-fidelity low-boom design based on near-field pressure signature,” in 54th AIAA Aerospace Sciences Meeting, 2016, p. 2033. [10] I. Ordaz, M. Wintzer, and S. K. Rallabhandi, “Full-carpet design of a low-boom demonstrator concept,” in 33rd AIAA Applied Aerodynamics Conference, 2015, p. 2261. [11] J. A. Heyns, O. F. Oxtoby, and A. Steenkamp, “Modelling high-speed flow using a matrix-free coupled solver,” in Proceedings of the 9th OpenFOAM Workshop, Zagreb, Croatia, 2014, pp. 23–26. [12] X. Gao, C. Xu, Y. Dong, M. Xiong, D. Li, Z. Wang, and X. Deng, “Developing a parallel density-based implicit solver with mesh deformation in openfoam,” Journal of computational science, vol. 28, pp. 59–69, 2018. [13] S. Yoon and A. Jameson, “Lower-upper symmetric-gauss-seidel method for the euler and navier-stokes equations,” AIAA journal, vol. 26, no. 9, pp. 1025–1026, 1988. [14] F. Mazzelli, A. B. Little, S. Garimella, and Y. Bartosiewicz, “Computational and experimental analysis of supersonic air ejector: Turbulence modeling and assessment of 3d effects,” International Journal of Heat and Fluid Flow, vol. 56, pp. 305–316, 2015. [15] B. Ma, G. Wang, J. Ren, Z. Ye, and G. Zha, “Near field sonic boom analysis with huns3d solver,” in 55th AIAA Aerospace Sciences Meeting, 2017, p. 0038. [16] J. Chu, Experimental surface pressure data obtained on 65° delta wing across Reynolds number and Mach number ranges. National Aeronautics and Space Administration, Langley Rearch Center, 1996, vol. 1 - Sharp Leading Edge. [17] T. Di Fabbio, E. Tangermann, and M. Klein, “Analysis of the vortex-dominated flow field over a delta wing at transonic speed,” The Aeronautical Journal, p. 1–18, 2023. [18] A. H¨ ovelmann, A. Winkler, S. M. Hitzel, K. Richter, and M. Werner, “Analysis of vortex flow phenomena on generic delta wing planforms at transonic speeds,” in New Results in Numerical and Experimental Fluid Mechanics XII. Springer International Publishing, 2020, pp. 307–316. [19] K. Rajkumar, T. Di Fabbio, E. Tangermann, and M. Klein, “Physical aspects of vortex-shock dynamics in delta wing configurations,” Physics of Fluids, vol. 36, no. 6, 2024. 50 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria CO2 PHASE CHANGE SIMULATIONS IN A CROGENIC INSULATED TANK BY OPENFOAM SEYON AHN1, SUNHO PARK2 1Korea Maritime and Ocean University, [email protected] 2Korea Maritime and Ocean University, [email protected] Keywords: Crogenic, CO2, OpenFOAM, Phase change As Carbon Capture and Storage (CCS) technologies have recently emerged as a key part of climate change mitigation strategies, reliable and efficient storage, and transportation of liquefied carbon dioxide (CO2) has become an important area of research. Cryogenic insulated tanks are used as an essential device for CO2 storage, but the phase change phenomena that occur inside the tank directly affect boil-off gas (BOG) generation and storage efficiency. In particular, the phase change process due to changes in temperature and pressure in the tank can affect the safety and long-term operational efficiency of the storage system, so quantitative analysis is essential. However, this phenomenon is a transient problem involving multiphase flow, heat transfer, and phase change, which is limited to be analyzed experimentally. In this study, computational fluid dynamics (CFD) was used to accurately analyze the phase change phenomenon of liquid CO2 in a cryogenic insulation tank, and analyze the characteristics of BOG generation and its effect on storage performance. Based on the simulation results, the optimal operating conditions to maximize CO2 storage efficiency was derived. OpenFOAM, an open-source computational fluid dynamics (CFD) platform, was used as the simulations [1]. Figure 1 shows the geometry used in the simulations. The insulation is wrapped around the tank and the total number of unstructured meshes is 540,000. The external temperature of the insulation is set to 20°C to account for ambient temperature. In Figure 2, the initial state of liquefied carbon dioxide is set to 94.5% of the total tank volume. Figure 1: Tank shape and typical mesh Figure 2: Initial liquid CO2 filling level (94.5%) 51 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria To evaluate the accuracy of the simulation results, convergence was tested using three different grids: coarse, medium, and fine, with coarse consisting of 250,000 unstructured grids, medium consisting of 540,000 unstructured grids, and fine consisting of 1.2 million unstructured grids. The convergence test was performed by measuring the pressure at four locations as shown in Figure 3. The results in Figure 4 show that the simulation converges sufficiently, so for the sake of computational time and economy, the simulation was performed using the medium case. It was found that the heat transfer between the tank and the insulation was well simulated due to the lowered liquid and vapor interface and the temperature of the tank caused by the BOG over time. Figure 3: Prove locations Figure 4: Pressures at each prove locations Acknowledgements This research was supported by the National Research Foundation of Korea (RS-2021-NR065941) and the Korea Planning & Evaluation Institute of Industrial Technology (20025023, RS-2024-00432033) References [1] W. Jeon, S. Park, S.K. Cho, Moored motion prediction of a semi-submersible offshore platform in waves using a OpenFOAM and MoorDyn coupled solver: International Journal of Naval Architecture and Ocean Engineering, 15, 100544 52 Tesla Cybertruck Aerodynamics - A fully automated reverse engineering methodology using 3D Scanning and OpenFOAM Wouter Remmerie1and Nikola Majksner2 1CEO, AirShaper, Antwerp, Belgium 2Research and Development, AirShaper, Antwerp, Belgium E-mail: [email protected] Abstract. Analyzing and understanding the aerodynamics of existing cars is key to improving the shape and efficiency of new cars. 3D scanning is commonly applied to obtain 3D models of cars, but these models pose significant challenges for CFD (computational fluid dynamics) methods - these typically require simplified and watertight models. This paper presents key aspects of a fully automated and validated workflow to analyze the aerodynamics of 3D scans without any manual efforts. The workflow has been validated against wind tunnel test data and has been endorsed by OEMs. 1 Problem Statement Aerodynamics are key when it comes to reducing fuel consumption of ICE (Internal Combustion Engine) cars or extending the range of BEV (Battery Electric Vehicles). With many automotive players working on similarly shaped vehicles, it common for manufacturers to reverse engineer competitor vehicles. In terms of aerodynamics, CFD (computational fluid dynamics) simulations performed on accurate 3D scans of vehicles can provide a fast and reliable alternative or complement to time-consuming and costly wind tunnel testing. Such 3D scan files, however, pose an enormous challenge for CFD simulations, as they are typically composed of a large number of non-watertight and/or non-manifold volumes featuring irregular surfaces. Manually repairing these CAD models can require weeks of work, which requires a lot of time and resources. With time to market being crucial to the success of OEMs and with Research and Development budgets under constant pressure in the highly competitive automotive market, this can become a strong challenge. In this paper, a number of key elements and learnings are highlighted concerning a fully automated and proven workflow to run CFD simulations on such models. This approach uses AirShaper [1], an online CFD platform based on OpenFOAM, and the 3D scan used as an example was provided by A2MAC1 [2]. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 53 Figure 1: Visualization of the STL file of the Cybertruck 2 3D model 2.1 Mesh visualization For this paper, a 3D scan of the Tesla Cybertruck was provided by A2MAC1. A2MAC1 offers various 3D scan files: •The exterior scan of the car, which is made before disassembly. This scan features a relatively low number of components. Figure 1 illustrates the surface mesh of this STL file. •The full 3D scan of the car, for which the car is completely disassembled. Each individual component is then scanned separately and digitally put back together (as shown in Figure 2). This can result in models with around 2.000 individual components. 2.2 Challenges for CFD simulations Both methods offer highly precise data, but these 3D models come with a number of challenges for CFD simulations: •File size: the resulting 3D files (usually in STL or OBJ format) are often very large - 2GB and more. This slows down or even prohibits loading such 3D files on machines with limited memory capacity. •Number of components: some of these scans can have thousands of individual components. This requires the boundary conditions for the CFD simulations to be defined for each individual patch. Doing this manually would be very time consuming. Also, meshing each of these components separately can pose memory problems if this process is not properly dealt with. •Surface quality: 3D scan data is based on spatial measurements of points on a surface. This means the resulting surfaces are never perfectly flat or smooth. When working with scans of a quality level below that of the Cybertruck example, this can pose strong challenges for the meshing algorithm to obtain a smooth surface mesh. •Non-watertight and non-manifold geometry: when a 3D surface is not perfectly closed (imagine a sphere with a hole in the surface), then it is considered non-watertight. The edges of such a hole do not have a face on each side, and are considered non-manifold. Conventional meshing methods will not work on such models, as it is not possible define the boundary beyond which certain volumes are not supposed to be meshed. Figure 3 highlights the large number non-manifold edges of the example 3D scan file. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 54 Figure 2: 3D scan file of the Tesla Cybertruck in STL file format •Overlapping geometry: various parts of the mesh, especially when reconstructed based on individual scans, can overlap. Such overlapping geometry can again cause meshing problems, as different surface meshes cross into each other. •Geometry variations: each car is different and thus has a unique aerodynamic flow pattern. For optimal results, a unique mesh needs to be generated for each geometry, taking not only the geometry variations but also the specific flow pattern into account. •Rotating elements: 3D scanned data, by default, does not contain any intelligence as to which components belong to a rotating wheel assembly, where the axis of rotation is located, what the tire radius is, and so on. Setting this up manually can be a very time consuming and error-prone process. 2.3 Traditional solutions Traditionally, these problems can be solved as follows: •Manual CAD repair: manually repairing the CAD geometry by stitching 3D surfaces together, closing holes, etc. This can be very time consuming, which increases the time to market. •automated CAD repair: techniques such as shrink wrapping (the digital equivalent of a vacuum bag) can be used to close holes. But these techniques are arbitrary, closing holes which should not be closed and vice versa. This also leads to larger, more monolithic 3D models. With less insights into the forces per component, this also reduces the value of a CFD simulation. In the following section, a number of key aspects of the fully automated AirShaper workflow are highlighted. 3 Solution The cloud-based AirShaper platform is built on top of OpenFOAM and Paraview. OpenFOAM is the most widely used open source CFD environment and paraview is the most widely used open source scientific 3D data visualization software. Both allow for a large degree of automation and customization. The AirShaper workflow is proprietary and confidential. Nevertheless, through this paper, some of the key elements and learnings are presented on how a fully automated CFD workflow compatible with any given 3D model was built using a.o. the tools mentioned above. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 55 Figure 3: Visualization of the non-manifold edges of the STL file 3.1 3D visualization Loading 3D models in a CFD software environment can be very demanding in terms of memory. 3D files of 2GB and larger can sometimes take 15 minutes or more to load, or not load at all of the computer of the user does not have sufficient memory. In the presented workflow, 3D models are decimated in the back end, on the AirShaper servers, to reduce their file size. This decimation involves the reduction of vertices and faces of the 3D models until a certain maximum number of faces is reached. The challenge in this process is to maintain an equal surface quality across the different components of the 3D scan files without merging them into one large 3D file, losing the component structure. Therefore, the workflow features methods to apply the decimation across the entire model while still grouping the resulting decimated faces per component. The compressed 3D models are then further compressed for online viewing using the Draco compression, again preserving the component structure. The result is a compressed model which, to the user, doesn’t look significantly different from the original model. This allows them to set up the simulation on a simple laptop with little memory, or even on a mobile phone. Once a simulation is launched, the uncompressed model will be used on the AirShaper servers for the actual CFD simulation. 3.2 Large number of components Having thousands of individual components requires an automated approach to naming these components, so that afterwards, the forces for each component can be calculated. This requires separate surface meshes for each component. This can lead to very large memory requirements, especially when the meshing process is executed in parallel across different processors. To avoid loading the full geometry on each processor, the automated processes uses the distributedTriSurface functionality to reduce the memory overall memory load. 3.3 Component splitting Very often, 3D files are stored in STP or OBJ file format - these remember the component hierarchy, which can then be copied by AirShaper. But many 3D scanning softwares export files in STL format. Very often, individual parts or components are not stored separately, but merged into one large mesh in this file format. However, to properly define elements such as rotating wheels or radiators, it’s essential to maintain some sort of component definition. To this end, AirShaper has implemented a component splitting algorithm. This method will automatically loop across the faces (triangles) of an STL file to identify isolated components. Via this The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 56 Figure 4: Result of the component splitting algorithm - top view method, it is often possible to regain some of the component structure that was lost during the export to STL (or other) formats. This, in turn, allows for a more precise setup of the simulation. Figures 4 and 5 illustrate the result of this process applied to the 3D STL file of the Cybertruck. As can be seen, this results in a large number of isolated components, with the wheels being separate ones. This is crucial, as this now allows for the definition of rotating wheels. 3.4 Non-watertight geometries To avoid losing assembly intelligence and errors in closing gaps and holes, AirShaper opted to not fix 3D models at all. Instead, the entire workflow has been modified in such a way that it is compatible with non-watertight 3D models. This has large impacts on pre-processing, solving and post-processing. 3.4.1 Pre-processing During the meshing process, mesh will leak to the inside of non-watertight geometries. Take the example of a car, where there is a small panel gap between the door and the A-pillar, without any rubber seal to virtually seal this gap. If the mesh is fine enough, it will connect the outside air to the inside air inside the car via this small gap. This has a number of implications: •Mesh inside the car: the inside of the car will also get meshed. This costs more in terms of meshing, but it is especially expensive when it comes to solving the simulation. The result is a static pocket of air, inside the car, being solved but not contributing to the actual aerodynamics of the car (assuming the gap is very small). •Thin surfaces: the scan usually represents a surface, and not a solid body. That means a surface will have flow on both sides - so patches need to be defined on both sides (inside and outside), each with their own wall treatment. •Surface edges: it’s important to accurately capture the mesh around these open (non-manifold) edges. That means extracting them as sharp edges to be used in surfaceFeatureExtract for local refinement. Also, extra loops in the snappyHexMesh to snap to such sharp features can be added to improve this. 3.4.2 Solving •Stability: if the gap inside mesh is connected to the external airflow via a gap in an area which is subject to large pressure variations, then this leads to large variations of pressure inside the car as The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 57 •Solving: There are hundreds of not thousands of possible combinations to solve a certain flow problem. Steady-state versus transient, RANS versus LES or DES or Lattice Boltzman, wall functions or prism layers, .. Each of these have a direct impact on the results. So when comparing a drag coefficient, it should be noted that for both a wind tunnel test and a CFD simulation, such a value is the result of a stack of errors ad uncertainties. Nevertheless, the drag coefficient obtained from the simulation is 0.344 differs by just 2.7% from the officially claimed Cd of 0.335. This certainly goes to show the values are very close. Even more so, it can reasonably be assumed that the official wind tunnel test was performed with aero covers. As these have been shown to reduce the drag on other cars, this gap is expected to be reduced even further. 4.5 Flow structures In addition to evaluating performance parameters like the drag coefficient, it’s important to also look at the flow structures of the simulation. Figure 13 illustrates the iso-surface for the total pressure coefficient for a value of 0. These clouds are very useful to understand the flow structure, as they indicate where energy is being lost. This visually corresponds fairly well with wake and flow separation areas: •Front air deflectors: the front air deflectors help to push air away around the front tires. This reduces the pressure on and flow separation around the front tires which feature a very coarse thread pattern. •Front hood: the sharp edge between the vertical nose of the car and the bonnet is a large source of flow separation. After some distance, however, the flow reattaches to the surface, as the surrounding air presses the flow against the bonnet & front wind shield. •Wheels: the front wheels feature a large wake, due to the fact that they are fairly exposed (because of the ride height, which is still fairly high) and their rough thread. The rear wheels feature a smaller wake, as they are partially in the wake of the front wheels and do not face the incoming air directly. •A-pillar: the very long A-pillars feature a sharp edge towards the side of the car, causing a very long and strong A-pillar vortex. As this vortex travels upstream along the side windows, it crosses over onto the roof & rear window of the car. This again results in a vortex, but one which rotates in the opposite direction. •Roof edge: the transition between the front wind screen and the roof of the car is again a sharp edge. The flow struggles a bit but on average remains fairly well attached to the roof. This is key, as the base shape of the Cybertruck, with its long downward slope at the rear, allows for beneficial pressure recovery - but only if the flow stays attached across the roof. •Mirrors: the mirros are quite large and as a consequence feature a local wake. •The underfoor has been optimized quite a lot to be as smooth and flat as possible. There are even covers on the suspension, so that they would form part of the larger flat underfloor and keep the flow attached all the way to the rear of the car. This accelerated airflow helps to reduce the size of the wake behind the car. For a more in-depth analysis of the aerodynamics of the Tesla Cybertruck, please refer to the detailed AirShaper video [5]. 4.6 Component forces Figure 14 illustrates how each separate component of the 3D model was meshed 5 Validation The method described in this paper has been applied to thousands of cases and has proven to be highly robust across different disciplines [6]. Specifically in the field of automotive simulations: •Tesla Model Y: A2MAC1 took the Tesla Model Y to the wind tunnel. The delta between the simulated and measured Cd value was less than 6% for a mesh of just 10 million cells [7]. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 64 Figure 13: Iso-surface of the total pressure coefficient at value 0 - 3D view Figure 14: Analysis of the forces & moments per individual component The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 65 •Rivian R1T: a reverse engineering analysis of the Rivian R1T was carried out at a resolution similar to the Tesla Cybertruck and published in a video [8]. This video was endorsed by Rivian themselves: ”It’s also great to see what we do being validated by an independent 3rd party. Wouter Remmerie & his AirShaper team correctly identify many of the subtle aerodynamic mechanisms throughout R1T that make it so efficient! Nice work.” [9] - Andrew McAllister, senior manager fluid dynamics at Rivian. •Other well known car manufacturers have been using AirShaper’s automated approach for years [10] - for example: –NIO [11]: ”We have been using AirShaper at NIO (ONVO) Design to get quicker aerodynamic simulations that are faster than the traditional ones with very similar results. Also, AirShaper’s Design Advice [12] provides us a red and blue map of areas for potential improvements on Cd values giving Design Flexibility to play with to achieve Cd targets” - Roy Morlet, Studio Engineering, NIO (ONVO) –Morgan Motor Company [13]: ”We utilized AirShaper on XP1, our electric prototype. We reduced the Cd from 0.65 to 0.42, impacting range significantly. You’ll see the effect on our future cars!” - Matt Hole, managing director. –Aptera [14]: ”AirShaper helped us keep expenses low while still getting amazing CFD work done fast” - Chris Anthony, Co-CEO. •Vecto Trucks: the exact same automated workflow has also been applied to the European Union regulations on Heavy Duty Vehicle emissions. These regulations include reference 3D models to validate CFD methods. The AirShaper CFD results all are within the allowed error bands [15]. •the automated AirShaper workflow also includes the option to automatically morph geometries using the adjoint shape optimization technique on non-watertight 3D models: –Porsche Taycan: an A2MAC1 3D scan model of the Porsche Taycan was used to demonstrate the shape morphing tools of AirShaper on non-watertight 3D models [16]. –Volkswagen ID3: here too, an A2MAC1 3D scan of the Volkswagen ID3 was used to optimize the rear roof spoiler geometry. The aerodynamic drag of the car was reduced by more than 5% [17]. 6 Further work The method is fully functional, but further work still needs to be done to further cut the time to market and increase the efficiency and accuracy for customers. Some of this work includes: •Prism layers: currently, the default is to run simulations without prism layers and use wall functions. There already is a free-of-charge option available to work with wall functions, but further beta testing is required before releasing this as the default setting. •Transient simulations: these are already available but not via the automated process. The computational cost of transient simulations is roughly 10 times higher (this can vary a lot though) compared to a steady state simulation and therefore economically less interesting. Especially because steady-state simulations provide most of the accuracy and insights required for drag reduction, the main goal of most OEMs. •Artificial Intelligence integration: there is a lot of potential to integrate artificial intelligene into the automated workflow. Current research and software tools still require hundreds of simulations to train a model within a narrow band of geometries. But fast progress is being made in this field. •Wheel rotation: currently, a tangential wall velocity is used for the wheel rotation. In some cases, however, it’s more realistic to use a MRF (moving reference frame) technique for the rims. •More features: the automated workflow also includes the option to add radiators. Other upcoming features include the option to add inlets and outlets so that air intakes and exhaust outlets can also be included in the simulation. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 66 7 Conclusion Having a fully automated workflow for CFD (computational fluid dynamics) is key to reduce the time to market and to increase the consistency and comparability between simulations. Traditional CFD softwares put strict requirements on CAD models, however, increasing the time required for CAD preparation. Especially when analyzing 3D scan files of cars, CAD repair can be highly time consuming, making the process extremely slow and costly. In this paper, a fully functional and validated automated workflow has been presented. The result is a workflow which requires just 5 minutes of set up time and yet provides accurate simulation results, even and especially on highly detailed and accurate 3D scans of cars. This allows for a thorough competitor analysis, as well as the application to in-house developed 3D models, as both can be tackled using the same automated CFD approach. The paper has highlighted the complex challenges and solutions to cope with the effects of a high number of components, non-watertight and non-manifold 3D models, geometry variations in terms of meshing and variations in force convergence behavior. Techniques such as automated component splitting, adaptive mesh refinement and automated convergence detection have been highlighted as just some of many to make this automated workflow possible. 8 Acknowledgements AirShaper would like to thank A2MAC1 for providing highly detailed 3D scan models of dozens of cars. Not only for this research, but also as part of their joint offering to provide standardized aerodynamic CFD analyses as a benchmarking tool. References [1] AirShaper. AirShaper - Aerodynamics Made Easy. https://www.airshaper.com. Accessed on Mon, Dec 30, 2024. [2] A2MAC1 - Decode the Future. https://www.a2mac1.com/. Accessed on Mon, Dec 30, 2024. [3] AirShaper adaptive mesh refinement code repository. https://github.com/airshaper/adaptive-mesh-refinement. Accessed on Fri, May 10, 2024. [4] AirShaper convergence detection code repository. https://github.com/airshaper/openfoam-convergence-detection. Accessed on Fri, May 10, 2024. [5] AirShaper. Tesla Cybertruck Aerodynamics. https://airshaper.com/videos/tesla-cybertruckaerodynamics-exclusive-3d-scan-reveals-flow-secrets-through-cfd-simulation/KPO5HEjLCbM. Accessed on Mon, Dec 30, 2024. [6] AirShaper validation cases. https://airshaper.com/validation. Accessed on Fri, Dec 31, 2024. [7] Tesla Model Y Wind Tunnel Test. https://airshaper.com/validation/tesla-model-y-wind-tunnel-test. Accessed on Mon, Dec 31, 2024. [8] Rivian R1T Aerodynamics. https://airshaper.com/videos/rivian-r1t-aerodynamics-is-the-claimeddrag-coefficient-of-030-correct/LxXCOT5ID20. Accessed on Mon, Dec 31, 2024. [9] Rivian Endorsement. https://www.linkedin.com/feed/update/urn:li:activity:6978109945569316866/. Accessed on Mon, Dec 31, 2024. [10] AirShaper Testimonials. https://airshaper.com/testimonials. Accessed on Mon, Dec 31, 2024. [11] Nio. https://www.nio.com. Accessed on Fri, Dec 31, 2024. [12] AirShaper Design Advice. https://airshaper.com/videos/reducing-aerodynamic-drag-on-a-sportscar-without-expertise/Ac469ecXQHQ. Accessed on Fri, Dec 31, 2024. [13] Morgan Motor Company. https://morgan-motor.com/. Accessed on Fri, Dec 31, 2024. [14] Aptera Motors. https://www.aptera.us. Accessed on Fri, Dec 31, 2024. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 67 [15] AirShaper Vecto Truck Validation. https://airshaper.com/validation/cfd-vecto-trailer-certification-eu-co2-reduction. Accessed on Fri, Dec 31, 2024. [16] Porsche Taycan Shape Optimization by AirShaper. https://airshaper.com/videos/aerodynamicshape-optimization-the-adjoint-cfd-method/cZAhPQFINZ8. Accessed on Fri, Dec 31, 2024. [17] VW ID3 Shape Optimization by AirShaper. https://airshaper.com/research/from-scanned-cad-to-an-optimized-car. Accessed on Fri, Dec 31, 2024. The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 68 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria A Second-Order Accurate Topology Preserving Remap for ALE Methods Ali Shayegh1,3, Philip Cardiff2 1University College Dublin, ali.shaye[email protected] 2University College Dublin, philip.cardif[email protected] 3TensorFields, [email protected] Keywords: Data Transfer, Field Mapping, Topology-Preserving Remap, Arbitrary-Lagrangian-Eulerian We present a second-order accurate data transfer technique between two arbitrary polyhedral meshes with same topology. Data transfer, is a key element in Arbitrary-Lagrangian-Eulerian (ALE) simulation methods. The term ALE method does not refer to a unique simulation procedure. In this work, we offer an advective data transfer technique suitable for topology-preserving discrete rezoning ALE procedures, in which the solution is obtained in a time-marching manner following the steps shown in Figure 1. Start Solve Lagrangian Eq End Time? Rezone Criterion Met? no Rezone yes no Remap Stop yes Figure 1: Discrete rezoning ALE Discrete Rezoning ALE: Discrete rezoning ALE, the term which we will explain shortly, is particularly useful in the simulation of severe deformation. By adopting the material displacement as the unknown in such problems, the Lagrangian momentum equation in Einstein notation reads RΩρd2ui dt2dV=RΩ0NjsjidS, where ρis the mass density, uiis the displacement, tis time, Ωand Ω0are the current and the reference material regions respectively, Njis the surface normal to the boundary of Ω0, and sji is the nominal stress tensor [1]. As detailed in [2], this equation can be solved using an implicit updated Lagrangian Finite Volume approach, which means, at each time step, the equation is converted into a system of linear algebraic equations and solved. So far, the method is Lagrangian, however, in discrete rezoning ALE, one stops time marching before the end time, when the mesh quality is not acceptable due to introducing errors or terminating the simulation. At this stage, the mesh is updated to a higher quality. Producing the new mesh is called rezoning, and since it is not every time step, but only when needed, it is discrete. When the mesh is changed, one needs to transfer all fields from the old mesh (also known as the source mesh) to the new mesh (also known as the target mesh); this phase is called remapping. After rezone+remap, the Lagrangian solution on the target mesh is obtained again. This loop is iterated until the end time is reached, see Figure 1. Method and Results: In this work, we introduce a second-order accurate advective remap procedure, the details of which follows. First, we move the source mesh points towards obtaining the target mesh, then derive and solve the equation that governs the transport of fields in such a motion. By mesh motion, we mean corresponding every mesh configuration with a real number, which we call fictitious time, denoted by τ. The fictitious-time-rate of an example scalar field Treads d dτRVρTdV= RVρdT dτdV+RVρTni(vbi−vi)dS, where Vis an arbitrary region, as the mesh does not stick to the material in this motion. Considering the fact that Tis not a function of τ, and not so is ui,dT dτ= 0 and vi=dui dτ= 0, and the transport equation reduces to d dτRVρTdV=RVρTnivbidS, where vbi=dxbi dτ, and xbistands for the position of the region boundary. Assuming a spatially uniform density, the Finite Volume (FV) discretised approximation of the equation is (T V )N−(T V )O ∆τ=PfTSiubi , where Oand Nstand for ‘old’ and ‘new’ fictitious times, and the summation is done over all faces of the cell under consideration. When the Tat the RHS is TN, the equation for all cells forms a system of linear algebraic equations with TN for all cells being the unknown vector. By solving the system, Tfield over the new mesh is obtained. The algorithm is summarized in Figure 3. According to the flowchart, we input to the remap algorithm the source mesh, e.g., 69 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 2a, and the desired target mesh, e.g., Figure 2c. Having the two meshes, we obtain the displacement field that produces the target mesh from the source mesh, which is simply obtained by subtracting the coordinates of the corresponding points of the two meshes. Then at each time step, we move the source mesh only by a fraction of the displacement field obtained in the previous step, so that the mesh Courant number remains less than or equal to a user-specified value. Finally we solve the system of the fictitious transport equations to obtain the field over the new mesh, i.e., the remapped field. In OpenFOAM syntax, this reads fvm::ddt(T) == fvm::div(mesh.phi(), T). The loop iterates until the end time. We examine our remap method by mapping an arbitrary scalar field ρ, used in a similar work [3], ρ(x, y) = 1 + sin(2πx)sin(2πy). Field contours on the source and the target meshes shown in Figures 2b and 2d show good agreement. Figure 4 is the plot of the mapping error versus the mesh size for two methods, i.e., the current method and an inverse distance interpolative map, available in OpenFOAM meshToMesh class. Comparison shows that the current method produces lower error magnitudes, and also it is second-order accurate, while the interpolative map is only first-order accurate. (a) Distorted source mesh (b) A field over the source mesh (c) Smooth target mesh (d) Target field Figure 2: Source and target meshes and fields Start Source and Target Meshes Move Points Advect Fields End Time? no Stop yes Figure 3: Advective remap algorithm 0.000010 0.000100 0.001000 0.010000 100 50 25 Error Ncells 2nd Order 1st Order Advective remap Interpolative remap Figure 4: The order of the error Acknowledgments This publication has emanated from research conducted with the financial support of Science Foundation Ireland and I-Form Advanced Manufacturing Centre under Grant number 21/RC/10295 P2. For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission. Additionally, the authors want to acknowledge project affiliates, Bekaert, through the Bekaert University Technology Centre (UTC) at University College Dublin (www.ucd.ie/bekaert), and I-Form (www.I-form.ie). References [1] P. Chadwick, Continuum mechanics: concise theory and problems. London: Allen and Unwin, 1976, no. Book, Whole. [2] P. Cardiff, “A Lagrangian cell-centred finite volume method for metal forming simulation,” International Journal for Numerical Methods in Engineering. [3] K. Lipnikov and M. Shashkov, “Conservative high-order data transfer method on generalized polygonal meshes,” Journal of Computational Physics, vol. 474, p. 111822, Feb. 2023. 70 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria DEVELOPMENT OF A MULTIPHASE EULERIAN - PBM MODEL FOR BUBBLE TRANSPORT WITH APPLICATIONS TO AN ELECTROLYSER CELL ASWATH ASHOK1, SILVANIA LOPES2, DELPHINE LABOUREUR3, GERALDINE J.HEYNDERICKX4 1von Karman Institute for Fluid Dynamics, [email protected] 2von Karman Institute for Fluid Dynamics, [email protected] 3von Karman Institute for Fluid Dynamics, delphine.labour[email protected] 2University Gent, Laboratory for Chemical Technology, Geraldine.Heynderic[email protected] Keywords: Multiphase, Eulerian, Hydrogen, Interfacial models, Population Balance Model, Method of Moments Hydrogen produced from Electrolysers utilising renewable electricity is termed as ’Green’ Hydrogen. Alkaline Water Electrolysers (ALK) and Proton Exchange Membrane Electrolysers (PEME) are the two most mature Electrolyser technologies currently in use [1]. There have been improvements in design, materials and process configuration, but their energy efficiency is limited to 50% for more extreme operating conditions. In both ALK and PEME, one of the limiting factors for operating efficiency is the transport of Hydrogen/Oxygen bubbles away from the electrodes and into the bulk flow channel (part of the Flow Field Plate). Some of the issues of bubbles in Electrolysers are highlighted in the schematic in Figure 1. Figure 1: Schematic of Bubble coverage issues In general, the formation of bubbles can be classified into 2 ways; a homogenous and a heterogenous one. In the former route bubbles nucleate in open spaces, while in the latter bubbles nucleate from the surface of the electrode/CL. In electrolysers, the nucleation of bubbles is strongly heterogenous in nature and takes place in a supersaturated solution. Bubbles nucleate at the Catalyst Layer (CL) (see Figure 2) and can grow upto an average size of 70 - 120 microns [2]. The electrogenerated bubbles then grow, detach and are transported from the CL to the bulk channel. Numerous electrogenerated bubbles are accounted by a size distribution [3]. The discrete bubbles can be modelled using a Lagrangian or a Statistical approach. Modelling this distribution through a statistical approach is achieved by the Population Balance Equation (PBE) and is known as Population Balance Modelling (PBM) [4]. The PBE is a continuity statement that is characterised by a distribution function (size distribution for example). Figure 2: Schematic of PEM Electrolyser cell (left) [5] and transport phenomena at Catalyst Layer (right) [6] 71 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria In previous studies modelling was done at the macroscopic scale. These models had limitations, for example in properly defining closure models for transport, mass transfer or capillary effects [7]. A full 3D macroscopic model was developed by Olesen et al. [8]. They did not account for interfacial forces and electrochemical reactions. Most studies analyze the bubble behaviour based on electrochemical performance rather than from a fluid dynamics perspective. The available bubble transport models were primarily developed for the Anode side, often neglecting the Hydrogen evolution [9]. The models do not account for transport within the Porous Transport Layer (PTL) (Figure 2) and from the PTL into the channel (Figure 2). Furthermore, the existing models for electrolyser performance characterise the bubbles with a uniform distribution. The objective of this current study is to develop a Eulerian - Population Balance Model (PBM) in OpenFOAMv12 based on the built-in PBM model in the dispersed multiphase flow solver. A Eulerian (Two Fluid) model for Gas-Liquid flow is considered. To more accurately account for the distribution of bubbles (size), a Population Balance Model (PBM) is considered with the Eulerian model. The electrolyser geometry is simplified to a horizontal channel with two-phase dispersed bubble flow. Only the growth and transport of the bubble is considered in the presented model. The aim of our study is to understand the system, characterise the bubble behaviour and identify where the improvements are needed in the model. These can include using the Method of Moments to replace the PBE, tuning of the Eulerian model for interfacial forces and/or inclusion of boundary conditions for bubble nucleation. The model we develop is validated against the experimental work of Tomasoni [3]. Building on this work, The overall research objective is to address the nucleation, growth and transport of the Hydrogen/Oxygen bubbles from the Catalyst Layer through the PTL into the bulk flow channel, by developing a 3D macroscopic multiphase CFD model. References [1] F. Pavan, “Electrolysers - Energy System,” Oct. 2023. [Online]. Available: https://www.iea.org/energy-system/ low-emission-fuels/electrolysers [2] T. Nierhaus, “Modeling and simulation of dispersed two-phase flow transport phenomena in electrochemical processes,” 2009. [Online]. Available: https://www.researchgate.net/publication/38229203 Modeling and simulation of dispersed two-phase flow transport phenomena in electrochemical processes [3] F. Tomasoni, Non-intrusive assessment of transport phenomena at gas-evolving electrodes. Rhode-St-Gen` ese, Belgium: Ph.D. Thesis, VKI / VUB, 2010. [4] D. L. Marchisio and R. O. Fox, Computational Models for Polydisperse Particulate and Multiphase Systems, ser. Cambridge Series in Chemical Engineering. Cambridge: Cambridge University Press, 2013. [Online]. Available: https://www.cambridge.org/core/books/computational-models-for-polydisperse-particulate-and-multiphase-systems/ A218C23D43C03DF345AF894A4859F55A [5] M. Holst, S. Aschbrenner, T. Smolinka, C. Voglstatter, and G. Grimm, “Cost Forecast for Low Temperature Electrolysis,” 2021. [Online]. Available: https://www.ise.fraunhofer.de/en/publications/studies/catf.html [6] J. Lopata, Z. Kang, J. Young, G. Bender, J. W. Weidner, and S. Shimpalee, “Effects of the Transport/Catalyst Layer Interface and Catalyst Loading on Mass and Charge Transport Phenomena in Polymer Electrolyte Membrane Water Electrolysis Devices,” J. Electrochem. Soc., vol. 167, no. 6, p. 064507, Mar. 2020, publisher: IOP Publishing. [Online]. Available: https://dx.doi.org/10.1149/1945-7111/ab7f87 [7] X. Guan, J. Bai, J. Zhang, and N. Yang, “Multiphase flow in PEM water electrolyzers: a minireview,” Current Opinion in Chemical Engineering, vol. 43, p. 100988, Mar. 2024. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2211339823000928 [8] A. C. Olesen, S. H. Frensch, and S. K. Kær, “Towards uniformly distributed heat, mass and charge: A flow field design study for high pressure and high current density operation of PEM electrolysis cells,” Electrochimica Acta, vol. 293, pp. 476–495, Jan. 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0013468618322242 [9] J. Nie and Y. Chen, “Numerical modeling of three-dimensional two-phase gas–liquid flow in the flow field plate of a PEM electrolysis cell,” International Journal of Hydrogen Energy, vol. 35, no. 8, pp. 3183–3197, Apr. 2010. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0360319910001217 72 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria SPECTRAL UNCERTAINTY QUANTIFICATION OF PRESSURE LOSS THROUGH A 90° PIPE BEND USING SST K-OMEGA WEEKS, B.1, TABOR, G.2 1University of Exeter, Babcock International, [email protected] 2University of Exeter, [email protected] Keywords: Uncertainty Quantification, RANS, SST K-Omega, Internal Flow, Polynomial Chaos, 90° Pipe Bend This work looks to apply Uncertainty Quantification (UQ), using the OpenFOAM [1] implementation of the SST kOmega [2] Reynolds-Averaged Navier Stokes (RANS) solver, to determining the pressure loss / loss coefficient of flow through a 90° pipe bend. Using the extensive numerical and experimental data that exists for this flow, this work aims to produce validated results as probability distributions, from uncertain inputs to a Computational Fluid Dynamics (CFD) solver. Spectral UQ using Polynomial Chaos (PC) Expansion, with the DAKOTA [3] toolkit, is used to perform this analysis with less computational effort than a naïve Monte-Carlo surrogate. Flow Domain The flow domain consists of a 1-inch (25.4 mm) diameter pipe with a 90° bend of radius 46.5 mm. A length of 50 pipe diameters (50D) upstream allows the flow to fully develop while providing a straight pipe length unaffected by the bend, for comparison to the downstream length. One diameter downstream, the flow is mapped back to the inlet to provide an ‘infinite length’ for flow development. The pipe is treated as hydraulically smooth. The downstream pipe length is 55 diameters (55D), as Ito [4] showed 50D to be sufficient to resolve all post-bend effects, and the flow domain is sampled 5D from the outlet to exclude any variations induced by the fixed value boundary condition. The flow medium considered is water. The meshing used Salome [5] to create inputs for SnappyHexMesh [1]. Mesh Sensitivity & Validation A two-part mesh sensitivity analysis was used. An arbitrary ‘inner mesh’ resolution was used for the initial ‘boundary layer’ sensitivity analysis, and the resulting boundary layer resolution used in the subsequent inner mesh sensitivity. The boundary layer sensitivity analysis plotted pressure losses against average y+, and the inner mesh analysis against total cell count. The losses measured are:  TotalLoss – The pressure loss from 10 pipe diameters (10D) upstream of the bend to 50D downstream.  Loss30 – The pressure loss from a 30D length of pipe, taken upstream of the bend, starting 15D from the inlet to ensure fully developed flow and to avoid effects from the bend.  BendLoss – The total loss attributable to the bend, including post-bend effects. Calculated by subtracting 2 x Loss30 from TotalLoss. The primary Quantity of Interest (QoI). The pressure losses were plotted against empirical / calculated values over a range of Re for the chosen mesh, to allow comparison to other research data. The Darcy-Weisbach equation was used for Loss30 and Ito’s [4] equations for BendLoss. Empirical TotalLoss is calculated from these. BendLoss varied from the empirical values by less than 2.3 % for Re values between 34,000 and 60,000, validating the method’s accuracy for calculating the bend loss. Figure 2 - BendLoss vs. Reynolds Number Figure 1: Flow Domain and Mesh 73 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria TOWARDS EFFICIENT CHEBYSHEV POLYNOMIAL SMOOTHING FOR MULTIGRID SOLVERS MATEJ ˇ CORAK1, TESSA UROI ´ C2, HRVOJE JASAK3, 1Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] 2Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] 3The Cavendish Laboratory, Department of Physics, University of Cambridge, UK, [email protected] Keywords: Finite volume method, OpenFOAM, Polynomial smoothing, Chebyshev polynomials, Highly parallel algorithms The efficient solution of large sparse linear systems is essential in computational fluid dynamics (CFD) applications, where it serves as the foundation for computing numerical solutions of partial differential equations (PDEs). These numerical solutions are required for simulating fluid flows, heat transfer, and combustion in real-world scenarios. In CFD applications, the accurate modeling of physical phenomena relies on the discretization of continuous PDEs into algebraic systems, typically resulting in large sparse linear systems. The size of these systems is determined by the spatial discretization, where the number of finite volumes corresponds to the dimension of the coefficient matrix. As computational domains become more refined, the resulting sparsity patterns of the linear systems can significantly impact the efficiency of iterative linear solvers [1]. The algebraic multigrid (AMG) method [2] has proven to be an efficient and scalable parallel solver [3], but modern high-performance computing platforms introduce new challenges for its parallel performance. This study presents an algorithm that utilizes polynomial smoothing as a simple and effective way to enhance the efficiency and scalability of the AMG solver. As suggested in the literature, there is no universally “best” polynomial; the optimal choice depends on the eigenvalue distribution, which is rarely known a priori [4]. In particular, we demonstrate that the Chebyshev polynomial is attractive for parallel architectures due to its avoidance of memory-intensive inner product calculations [5]. The newly developed algorithm will be tested for single-phase, incompressible, laminar, and turbulent flows. The expected outcome of the proposed algorithm is a faster solution phase. We anticipate that the findings of this study will help address the challenges associated with using recent high-performance computers and that the proposed approach can serve as an effective solution procedure for solving large sparse linear systems in CFD applications. References [1] Tessa Uroi´ c and Hrvoje Jasak, Block-selective algebraic multigrid for implicitly coupled pressure-velocity system, Computers & Fluids, vol. 167, pp. 100–110, 2018, Elsevier. [2] Achi Brandt, Algebraic multigrid (AMG) for sparse matrix equations,Sparsity and its Applications, pp. 257–284, 1984, Cambridge University Press. [3] Allison H Baker, Robert D Falgout, Tzanio V Kolev, Ulrike Meier Yang, Scaling hypre’s multigrid solvers to 100,000 cores,High-performance scientific computing: algorithms and applications, pp. 261–279, 2012, Springer. [4] Steven F. Ashby, Thomas A. Manteuffel, James S. Otto, A comparison of adaptive Chebyshev and least squares polynomial preconditioning for Hermitian positive definite linear systems,SIAM Journal on Scientific and Statistical Computing, vol.13, no. 1, pp. 1–29, 1992, SIAM. [5] Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo, Charles Romine, Henk Van der Vorst, Templates for the solution of linear systems: building blocks for iterative methods, 1994, SIAM. 80 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria ADAPTIVE AND MULTI-OBJECTIVE LOAD BALANCING FOR DETAILED CHEMISTRY SIMULATIONS IN OPENFOAM ALBERTO CUOCI CRECK Modeling Lab, Department of Chemistry, Materials, and Chemical Engineering, Politecnico di Milano (Italy), [email protected] Keywords: detailed chemistry, reactive flows, load balancing, parallel efficiency The accurate simulation of reactive flows with detailed chemical kinetics, involving hundreds of species and thousands of reactions, remains a major computational challenge in computational fluid dynamics (CFD). When operator-splitting approaches are adopted, the chemical integration step often accounts for more than 95% of the total computational cost. Moreover, the highly localized and non-linear nature of chemical activity leads to severe load imbalances across parallel processes: cells near flame fronts or ignition kernels demand significantly more CPU time than inert regions. Standard domain decomposition techniques in CFD solvers, including OpenFOAM, partition the computational domain based on geometric criteria and assume uniform computational load per cell. Consequently, they fail to address the heterogeneity introduced by detailed chemistry, resulting in processors finishing at different times, idle cores, and poor parallel efficiency. As the number of MPI processes increases and the number of cells per core decreases, this imbalance further amplifies, severely limiting scalability. To mitigate this issue, several chemistry-specific dynamic load balancing (DLB) strategies have been developed. Notable efforts include DLBFoam [1], which dynamically redistributes chemistry tasks during runtime and applies reference mapping to exploit local thermochemical similarity; the chemistry load balancing library by Gärtner et al. [2], which optimizes task redistribution for both standard and TDAC chemistry models through efficient point-to-point MPI communications; and the method of Zirwes et al. [3], focusing on simple dynamic offloading of chemical computations between paired processors based on measured load. Recently, Van den Oord et al. [4] proposed a multi-level dynamic balancing approach targeting both chemical kinetics and real-fluid thermodynamics, achieving remarkable performance improvements through modular load redistribution stages. However, the existing techniques present certain limitations. Many rely on pairwise work exchange strategies that may become suboptimal at high core counts or suffer from increased communication overhead when scaling to massively parallel architectures. Furthermore, some approaches, while effective, are tightly coupled to specific solvers or require careful tuning of rebalancing parameters. In this context, we implemented and systematically evaluated several dynamic load balancing algorithms within the OpenFOAM framework, specifically targeting the chemical integration step in reactive flow simulations involving detailed kinetic mechanisms. The tested strategies include task migration schemes inspired by previously published methods, such as sender/receiver pairing [3], cost-sorted redistribution [4], and an optimization-based scheduling algorithm using a 1.5-approximation scheme [5]. All techniques were adapted to the operator-splitting framework of OpenFOAM and assessed in terms of load balancing effectiveness, communication overhead, and overall simulation performance. The objective is to achieve a more uniform chemical workload distribution across MPI processes, reduce idle time, and improve parallel efficiency in simulations where chemistry dominates the computational cost. Methodology The proposed dynamic load balancing technique is specifically conceived to address the challenges arising in reactive flow simulations involving very large detailed kinetic mechanisms, where the number of species can easily exceed several hundreds. In such cases, the chemical integration cost becomes extremely high and tends to scale more than quadratically with the number of species, making it the dominant component of the total computational time. The method builds upon the optimization-based scheduling strategy originally introduced by Wu et al. [5], and adapts it to the OpenFOAM framework, introducing key improvements to maximize performance and robustness under these specific conditions. Rather than attempting to rebalance the complete computational workload, the methodology focuses exclusively on the chemical step. Each computational cell is dynamically assigned a weight proportional to its estimated or measured chemical integration cost. The objective is to redistribute cells among MPI processes to minimize the maximum cumulative chemical load per process, thereby achieving near-uniform load balancing across the domain. The cell redistribution is formulated as a scheduling optimization problem, where the goal is to assign tasks (cells) to processors in a way that approximates the optimal solution within a known performance bound. A parallel 1.5approximation algorithm is used: heavy-weight cells are first allocated optimally to minimize peak load, while lightweight cells are assigned through a greedy procedure to fill remaining gaps. Compared to previously proposed OpenFOAM-specific approaches, the methodology introduces two major innovations. 1) Adaptive Rebalancing Frequency: rather than performing rebalancing at fixed intervals, the method monitors the chemical load imbalance dynamically during the simulation. A rebalancing operation is triggered only when the 81 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria monitored imbalance exceeds a predefined threshold (e.g., a maximum load imbalance ratio 𝜂 greater than 1.2). This adaptive strategy minimizes unnecessary data migration during quasi-steady-state phases, while responding rapidly to dynamic events such as flame propagation or ignition. 2) Multi-Objective Balancing Strategy: during the redistribution phase, the optimization does not only minimize the chemical load peak but also penalizes excessive data migration between processors. A simple multi-objective cost function is introduced, combining the chemical load with a migration penalty term, thus favouring partitioning solutions that achieve good load balance while containing communication costs. This ensures that the improvements in load balance do not come at the expense of excessive overhead, although in cases dominated by large chemistry the impact of communications remains secondary relative to the chemical integration cost. The methodology is especially suitable for simulations characterized by a moderate number of computational cells and a limited-to-moderate number of MPI processes. Unlike classical load balancing strategies that target massively parallel simulations with millions of cells, the focus here is on maximizing efficiency for highly detailed chemistry cases, where the primary challenge is the enormous computational weight of the chemical step rather than the sheer size of the mesh. Key features of the implementation include: i) lightweight data structures for storing cell-specific chemical load estimates; ii) efficient construction of send/receive maps during rebalancing phases; iii) non-blocking MPI communications for migrating thermochemical state data. Importantly, the method preserves the existing mesh decomposition for solving transport equations, ensuring compatibility with mesh-based load balancing strategies, adaptive mesh refinement (AMR), and minimizing intrusive code modifications. Results The proposed dynamic load balancing methodology has been tested in two challenging reactive flow configurations characterized by the use of highly detailed kinetic mechanisms, pushing the computational load and scalability limits typical of operator-splitting approaches. The first test case consists of a pulsating laminar coflow flame, where the fuel stream is periodically modulated to induce unsteady flame dynamics. The chemical kinetics are modelled using a detailed mechanism comprising ~160 species and ~2,500 reactions, including an extended NOx formation sub-mechanism. The simulation was carried out using ~0.1 million cells distributed across 32 MPI processes. In the absence of dynamic load balancing, the chemical step exhibited a load imbalance ratio 𝜂 exceeding 3, leading to substantial idle times among faster processes. By applying the proposed dynamic rebalancing strategy, the chemical load imbalance was reduced to below ~1.2, resulting in a total chemistry step time reduction of approximately 50%. The overall simulation speed-up, accounting for both chemistry and transport steps, was over 40%. Importantly, the dynamic rebalancing was able to maintain stability and accuracy without introducing numerical artifacts. The second test case involves the simulation of a 2D turbulent non-premixed flame in decaying isotropic turbulence [6]. The chemical mechanism includes detailed soot formation and growth modelling via a Discrete Sectional Method (DSM), involving ~320 species and more than ~14,000 reactions. The computational domain contained ~0.25 million cells, and simulations were performed on 128 MPI processes. Due to the highly localized nature of soot chemistry and the evolving turbulence structures, the chemical step load imbalance initially exceeded a ratio of 4. With the dynamic load balancing activated, the imbalance was reduced to ~1.15. Consequently, the chemical integration wall time decreased by more than 50%, and the overall simulation efficiency improved by approximately 50%. The dynamic method effectively adapted to the evolving turbulent structures, consistently redistributing the chemical load as the flame and soot production regions evolved. Acknowledgements The author acknowledges EuroHPC Joint Undertaking for awarding access to Discoverer at SofiaTech, Bulgaria (Project. EHPC-DEV-2024D12-015). References [1] B. Tekgül, P. Peltonen, H. Kahila, O. Kaario, V. Vuorinen, “DLBFoam: an open-source dynamic load balancing model for fast reacting flow simulations in OpenFOAM”, Comput. Phys. Commun., vol. 267, 2021, Art. no. 108073. [2] J.W. Gartner, A. Shamooni, T. Zirwes, A. Kronenburg, “A chemistry load balancing model for OpenFOAM”, Comput. Phys. Commun., vol. 205, 2024, Art. no. 109322. [3] T. Zirwes, F. Zhang, P. Habisreuther, J.A. Denev, H. Bockorn, D. Trimis, “Optimizing Load Balancing of Reacting Flow Solvers in OpenFOAM for High Performance Computing”, 6th ESI OpenFOAM User Conference 2018, Hamburg, October 23-25, 2018. [4] G. van der Oord, V. Azizi, M. Fathi, S. Hickel, “Dynamic multi-level load balancing for scalable simulations of reacting multiphase flows”, Int. J. High Perform. Comput. Appl., pp. 1-13, 2025. [5] H. Wu, P.C. Ma, M. Ihme, “Efficient time-stepping techniques for simulating turbulent reactive flows with stiff chemistry”, Comput. Phys. Commun., vol. 243, pp. 81-96, 2019. [6] F. Bisetti, G. Blanquart, M.E. Mueller, H. Pitsch, “On the formation and early evolution of soot in turbulent nonpremixed flames”, Comb. Flame, vol. 159, pp. 317-335, 2012. 82 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria DESIGN OF LOW REYNOLDS NUMBER AIRFOILS CONSIDERING MANUFACTURABILITY HYOJUN KIM1, CHANKYU SON2 1Cheongju University, [email protected] 2 Cheongju University, cso[email protected] Keywords: Airfoil design, Low reynolds number, Manufacturability Introduction In recent years, the application of multi-rotor unmanned aerial vehicles(UAVs) has expanded significantly across both civilian and military sectors, including logistics, surveillance, and reconnaissance missions. These UAVs typically operate in low Reynolds number regimes, where achieving high aerodynamic efficiency often requires complex airfoil geometries or numerous design parameters [1, 2]. However, due to structural constraints such as short propeller diameters and thin blade sections, implementing such complex airfoil designs demands extremely high manufacturing precision. Previous studies have primarily focused on evaluating the aerodynamic performance of airfoils and enhancing performance through the exploration of various design parameters. However, this approach leads to a rapid increase in the number of configurations to be analyzed, resulting in substantial computational resource demands. Furthermore, relatively little attention has been given to practical considerations such as manufacturability and repeatability in mass production. To address these limitations, the present study proposes a set of simplified design variables that simultaneously consider both aerodynamic performance and manufacturability. Design Parameters In this study, an airfoil design approach was proposed based on geometric parameters with constant thickness and curvature, while incorporating angular definitions at the leading and trailing edges, in consideration of manufacturability. Based on prior research[3], which identified shape parameters significantly affecting aerodynamic performance, four primary design variables were selected: maximum camber, camber location, thickness, and leading edge angle, as illustrated in Fig. 1. In addition, the trailing edge angle was included as a design variable, as it is considered to influence aerodynamic characteristics such as flow separation, wake structure, and vortex generation. Accordingly, airfoil analysis was conducted based on these five design variables under the conditions presented in Table 1. The shape definition was carried out by constructing the camber line using a B-spline method based on maximum camber and camber location, while the airfoil outer geometry was generated using thickness, leading edge angle, and trailing edge angle. Figure 1: Airfoil design parameters Table 1: Geometric parameter ranges and constraints Re Mach θleading θ trailing Camber location Max camber Thickness 1 ×104 0.1 20°–80° 20°–180° 30 %–70 % 2 %–10 % 4 %–7 % Validation of Numerical Approach and Results The computational grid was generated using Pointwise with an O-grid topology. The initial grid height was set to 1.69×10⁻⁴ to ensure that y⁺ < 1 was satisfied on all surfaces of the airfoil. To account for laminar-to-turbulent transition effects, the Langtry–Menter k-omega Shear Stress Transport (SST) turbulence model[4] was employed. For the validation of the numerical approach, simulations were conducted on the NACA 0012 airfoil, for which experimental data under a Reynolds number of 4×10⁴ are available[3]. The results are presented in Fig. 2. This study analyzed the effects of five design parameters used in airfoil development on aerodynamic performance. The results showed that a thinner airfoil profile led to a reduction in form drag, thereby improving overall aerodynamic efficiency. In cases where the leading edge was sharp, a tendency for drag reduction was observed at low AoA. However, at high AoA, performance degradation due to early flow separation was identified. Additionally, when the camber location was located closer to the leading edge, flow attachment was better maintained, which resulted in improved L1.5/D performance. Notably, the highest efficiency was achieved under the condition of 6% maximum 83 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria camber. Regarding the trailing edge angle, the 180° configuration helped suppress severe flow separation and strong vortex formation as the flow passed over the trailing edge. This phenomenon is interpreted as the result of a trailing edge stagnation point being formed to satisfy the Kutta condition. Based on these findings, the airfoil configuration with 4% thickness, 50° leading edge angle, 180° trailing edge angle, 30% camber location, and 6% maximum camber demonstrated superior L¹.⁵/D performance across most angles of attack. As shown in Fig. 3, the proposed airfoil improved performance over the conventional 6% cambered airfoil (with 1% thickness) across a wide range of angles. In particular, it exhibited peak performance at α = 4°, showing a more efficient lift-to-drag ratio under the same test conditions. These results suggest that the proposed airfoil provides enhanced aerodynamic characteristics in low Reynolds number regimes, owing to its optimized curvature and geometric configuration—even with increased thickness. Figure 2: Numerical validation with NACA 0012 (a) Cl comparison, (b) Cd comparison Figure 3: L1.5/D comparison of proposed and 6% cambered[3] Acknowledgements This research was supported by the Stratospheric Drone Technology Development Program through the National Research Foundation of Korea(NRF) and Stratospheric Drone Technology Development Center funded by the Ministry of Science and ICT, the Republic of Korea(No. 2022M3C1C7090617) References [1] H. Hu, M. Tamai. “Bioinspired Corrugated Airfoil at Low Reynolds Numbers,” Journal of Aircraft, Vol. 45, no. 6, pp. 2068-2077, 2018. [2] H. Sobieczky, "Parametric Airfoils and Wings," Notes on Numerical Fluid Mechanics, vol. 68, pp. 71–88, 1998. [3] J. Winslow, et al. “Basic Understanding of Airfoil Characteristics at Low Reynolds Numbers (104–105),” Journal of Aircraft, Vol. 55, no. 3, pp. 1050-1061, 2018. [4] R. Langtry, F. Menter, “Correlation-Based Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes”, AIAA journal, Vol. 47, No. 12, pp. 2894-2906, 2009. 84 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria HYBRID VOF TO LAGRANGIAN CFD MODEL FOR JET BREAKUP AND SPRAY TRANSPORT ULRICH HECK, MARTIN BECKER DHCAE Tools GmbH, [email protected] Keywords: Volume of Fluid (VOF), Lagrangian, multiphase flows, atomization, liquid jet in crossflow For a complete CFD modelling of jet breakup processes as in atomization followed by spray propagation, a transition model from an interface tracking consideration with the volume of fluid method (VoF) to Lagrangian particles is presented. The model was integrated by the authors into the open-source CFD toolbox OpenFOAM, validated on various benchmarks, and finally used to clarify disintegration processes in molten metal atomization. Modelling approach To represent the interactions between dynamic, viscous and surface forces during the disintegration, the initial fluid dispersion is modelled using a volume-of-fluid approach. Here, adaptive meshes on the liquid-gas interface are used to geometrically resolve the lamellae and ligament structure. If the individual cell regions of disintegrated ligaments form separated and spherical structures, they are converted into Lagrangian particles. The transformation is required, because the many droplets produced during atomization are often too small for a further spray modelling using interface tracking methods. Typical criteria to identify VoF sections valid for a transition to Lagrangian particles are the sphericity, size and position of separated fluid cell regions. After the transformation the mesh is locally coarsened again and the interactions between the droplets and continuous flow are realized via Lagrangian source terms. In the following, the droplets can further disintegrate according to their respective Weber number in a secondary breakup model. The model has been integrated into OpenFOAM as cloud-functions. Thus, it can be used with both incompressible and compressible solvers, with compressible solvers being particularly suitable for twin-fluid atomization at high gas velocities. Validation For validation the fuel jet in cross flow benchmark is set up [1] [2]. A LES model is used to represent the turbulence structure and an iso-advector interface tracking algorithm is applied [3]. The horizontally flowing air atomizes the liquid whereby the final drop sizes at a greater distance are dominated by the secondary break-up of the Lagrange droplets according to the ReitzDiwakar model. Good agreement with experimental results is achieved for the resulting particle size in the control planes at different distances from the injection nozzle. Figure 1: Fuel injection in cross flow – Experiment [1] [2] and simulation Additional applications, such as pressure atomization from round nozzles, also yield accurate results for droplet sizes and velocities, see Figure 2. The comparison is made here with experimental data from the thesis of E. Deux [4]. The simulation in OpenFOAM is also carried out with an LES model. In contrast to the fuel injection benchmark, a singlecomponent atomisation is used. For this reason, the secondary break-up (also modelled with the model according to ReitzDiwakar) is of minor importance. 85 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 2: Atomisation from a hole nozzle - experimental droplet sizes and velocities [5] at a distance of 25mm from the exit in comparison to simulation Molten metal atomisation The model developed in this way is used to explain processes in the atomisation of metal melts during powder production [5]. In this process, a pressure swirl atomiser is first applied to the molten metal, which causes the lamella to primarily disintegrate. Subsequently, the droplets are atomised into small particles by gas jets at high speed, see Figure 3. Due to the transonic gas velocity, compressible modelling is used. With this model, significant effects on the breakup processes could be better understood. In particular, the interaction of the gas jets with the ligaments and effects on the particle size could be analysed. Figure 3: Atomisation of molten metal in twin fluid configuration References [1] Gopala Y., Zhang P, Bibik O, Lubarsky E, Zinn B. “Liquid fuel jet in crossflow trajectory correlations based on the column breakup point”. In: 48th AIAA aerospace sciences meeting. Orlando, Florida; 2010. [2] Sekar J., Rao A., Pillutla S., Danis A., Hsieh S-Y. “Liquid jet in cross flow modelling” In Proceedings of ASME turbo expo 2014: turbine technical conference, Düsseldorf, Germany; 2014 [3] Roenby J., Bredmose H, Jasak H. “A computational method for sharp interface advection”. Royal Society Open Science, 2016, 3(11): 160405. [4] Deux E. „Berechnung der turbulenten Zerstäubung von Flüssigkeiten durch Kombination eines Zweifluidmodells mit dem Euler-Lagrange-Ansatz“, Dissertation Halle-Wittenberg, 2006 [5] Kamenov D, et al., “Investigating the Atomizer Performance within Aluminium Melt Atomization” The European Conference on Liquid Atomization & Spray Systems (ILASS), 2022. 86 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria CLOUD-BASED AUTOMATED CFD DESIGN TOOLS USING THE EXAMPLE OF THE HEAT TREATMENT OF TITANIUM COMPONENTS ULRICH HECK, MARTIN BECKER DHCAE Tools GmbH, [email protected] Keywords: CHT, Coupling, Multi-region, Workflow automation While CFD methods for simpler flow conditions have already found their way into product and process development, complex problems such as multiphase or multiphysics applications often still require intensive model validation and, due to the usually long computing times, model optimisation. One way to provide the developer with a design tool for these more demanding CFD tasks is to create an automated CFD tool for a specific process that is validated and optimised for the application. This makes it particularly easy to use, especially in conjunction with cloud resources, as the frequently required high computing resources can be made available as needed. This means that HPC systems can be addressed in the background under Linux and the administration of hardware systems on site is no longer necessary. Modelling approach Using the example of modelling the cooling behaviour of titanium components during heat treatment, the creation and application of such a tool to reduce creep processes during cooling is demonstrated. Natural convection, heat conduction, radiation exchange and energy release due to microstructural transformation are taken into account in a transient conjugate heat transport analysis. The CFD model was created in OpenFOAM, validated by numerous experiments and optimised in terms of computing time. Figure 1: CHT model in CFD analysis Validation Figure 2 shows the normalized temperature curve and the cooling rate of some comparison positions. In general, a good agreement between simulation and experiment is achieved. Various effects can be reproduced in the simulation, especially for the more important cooling rates in the right-hand picture: the cooling rates are characterized by an initial sharp drop due to the high energy dissipation to the environment. In the second step, there is a strong decrease in the cooling rates, since energy is released in the component due to the structural transformation. In the following, the cooling rate in the different positions depends strongly on the environmental conditions: Figure 2: Temperature and cool ingrate vs. relative time for cooling in quiescent air in different positions of the component 87 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria At the outermost component position 8 of the thermocouple, cooling occurs faster than in the centre positions (e.g. position 5), because on the one hand there is less thermal mass here and on the other hand the radiation exchange with the grate is reduced due to the smaller cross-sectional area. The effects of the position-dependent cooling rates in the experiment are thus very well reproduced by the simulation. This is the crucial prerequisite for optimizing the cooling process, since the creep of the component is finally attributed to uneven cooling of the component. In addition to cooling in quiescent air, cooling under forced convection in a rapid cooling chamber was also analysed [1]. For creep modeling, the time-dependent temperature field from the CHT analysis is transferred to the Abaqus structural solver, where the creep analysis is performed. Here, it has been shown that unidirectional coupling is sufficient. Automation The CFD workflow is fully automated. The geometry exported in STL format for grate, supports and component serves as input. With a parameter file for the individual boundary conditions, the case is automatically meshed via Python and shell scripting, provided with boundary conditions and finally calculated. For the meshing, the OpenFOAM-internal mesh generator snappyHexMesh is used, which generates a hexahedron-dominant polyhedral mesh. The meshing and simulation are carried out on a local server or in the cloud under Linux. A monitoring tool allows the simulation progress to be checked from a Windows desktop. Finally, the mapping of the temperature data to Abaqus is also done via the monitoring tool, see Figure 3. Figure 3: Automated workflow for the creation and calculation of the CFD model Application The tool is used by the end user to optimise the cooling conditions of titanium components in the preliminary design phase. In particular, it was shown that component distortion due to creep is minimised if the component is cooled as evenly as possible on the top and bottom surfaces. In this context, the grate on which the component lies is of particular importance, as a large amount of energy is stored here, which is exchanged with the component over a long period of time through radiation. As a result, the component cools down more slowly on the underside. Numerical simulations can be used to test various measures in advance in order to compensate for these effects and ensure uniform cooling of the component References [1] Ulrich Heck, Martin Becker, Ralf Paßmann, Volker Hardenacke A coupled flow, heat and structural analysis tool to reduce creep during heat treatment of titanium components, NAFEMS Multiphysics Conference, Munich, 14.1.5 Nov, 2023 88 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria MULTISCALE FIXED BED REACTOR MODELING WITH INTERNAL PORE DIFFUSION: COMPARING 3D PRCFD AND 2D POROUS MEDIA SIMULATIONS ERIC DAYMO1, MATTHIAS HETTEL2 1Tonkomo, LLC, [email protected] 2 Karlsruhe Institute of Technology (KIT), [email protected] Keywords: Fixed bed reactor; packed bed reactor, pore diffusion; porous media; PRCFD Fixed (aka packed) bed reactors with porous catalysts are commonly used in the chemical process industries, hence the accurate modeling of these reactors is of interest. The catalytic pellets that comprise the bed are often highly porous to increase the internal surface area, which also increases the number of active catalytic centers. These active catalytic centers are generally accessible only via tortuous and narrow (e.g., micron scale or smaller) passages through which flow is controlled by diffusion, not convection (Fig. 1). Such diffusion limitations can lead to a difference in concentration between the surface of the pellet and the interior of the pellet, which can affect the observed reaction rates. To improve simulation accuracy, internal mass transfer (i.e., intra-particle pore diffusion resistances) must often be considered when modeling fixed bed reactors. Two of the most common models for internal pore diffusion are the effectiveness factor method and the 1D reaction-diffusion method. Figure 1: A schematic of the internal pore structure of a typical fixed bed reactor catalyst pellet The effectiveness factor is the ratio of the reaction rate with and without internal diffusion limitations [1]. The rates for each specie i are modified by the effectiveness factor, 𝜂 , such that 𝑠 # !,#$$ = 𝜂%𝑠 # ! . Here, 𝑠 # !,#$$ is the effective surface reaction rate accounting for concentrations within the particle differing from concentrations at the surface of the catalyst pellet, and 𝑠 # ! is the surface reaction rate calculated with concentrations at the (outer geometrical) surface of the catalyst. When there are no external mass transfer resistances across the boundary layer between the bulk flow and the catalyst surface, the surface concentrations are known from the CFD solution of the species transport equation. 𝜂 can then be calculated using an analytical solution. However, when a more rigorous approach is needed to address internal mass transfer resistances, the 1D reaction-diffusion equation (Eq. 1, [1]) can be solved directly (on a subgrid) at each catalytic face (for a surface reaction) or cell (for reactions in porous media). For specie i, 𝐷!,#$$ is the effective diffusivity [1], and 𝐶! is the specie concentration. −𝛻 ∙ + 𝐷!,#$$ ∙ 𝛻𝐶! , = 𝑅! (1) The chemistry library DETCHEM [2] solves for the effective surface reaction rate using either the effectiveness factor or 1D reaction-diffusion equation approaches. DETCHEM is coupled with OpenFOAM [3] via the DUO solver (DETCHEM und OpenFOAM), allowing for access to these pore diffusion approaches within OpenFOAM. In prior work [4],[5],[6], DUO modeled packed bed reactors either with a porous media model (using packed bed heat transfer correlations) or with Particle Resolved Computational Fluid Dynamics (PRCFD, where each individual catalyst pellet is resolved in the grid). However, prior work with DUO did not consider internal pore diffusion effects. Also, it was not known whether PRCFD would still match 2D porous media models when pore diffusion modeling is enabled. The present work addresses this question by comparing PRCFD with a 730-particle bed against porous media models for methane Catalytic Partial Oxidation (CPOX) (Fig. 1), Dry Reforming of Methane (DRM), and Steam Methane Reforming (SMR). CPOX is exothermic, but DRM and SMR are endothermic reactions. Also, microkinetics models are used for CPOX and DRM, while a global reaction mechanism is used for SMR. Examples CPOX results are shown below. Fig. 2 showcases 3D PRCFD results while Fig. 3 compares 3D and 2D (porous media) results along the symmetry axis of the bed when pore diffusion modeling is disabled. Meanwhile, in Fig. 4. the Diffusion through Boundary Layer Active Centers Catalyst particle 89 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria References [1] S. Bali, Ž. Tuković, P. Cardiff, A. Ivanković, and V. Pakrashi, "A cell-centered finite volume formulation of geometrically exact Simo–Reissner beams with arbitrary initial curvatures," International Journal for Numerical Methods in Engineering, vol. 123, no. 17, pp. 3950-3973, 2022, doi: https://doi.org/10.1002/nme.6994. [2] A. Taran, S. Bali, Ž. Tukovic, V. Pakrashi and P. Cardiff, “Development and validation of a finite volume Simo-Reissner beam method for moored floating body dynamics”, doi: https://doi.org/10.48550/arXiv.2504.18248 [3] M. Campbell, G. Baruah, M. Karimirad, M. Van, and V. Pakrashi, “Experimental assessment of a semisubmersible floating solar platform subjected to wave and wind loading.” Innovations in Renewable Energies Offshore: Proceedings of the 6th International Conference on Renewable Energies Offshore (RENEW 2024) pp. 119-128, 2024, doi: https://doi.org/10.1201/9781003558859-14 96 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria APPLICATION OF ADVANCED ACTUATOR LINE METHOD FOR COMPRESSIBLE PROPELLER SIMULATIONS IN OPENFOAM LOUIS FLIESSBACH1, MOHAMED SUHAIL ERAMASRAYILLATH ABDULSEMATH2, HENDRIK HETMANN1, THILO KNACKE1, GUILLAUME MARTINAT2, CHARLES MOCKETT1, NILAKANTHA SAHOO1 1Upstream CFD GmbH, [email protected] 2FLYING WHALES, [email protected] Keywords: Propeller, Actuator Line Method, Aerodynamics, OpenFOAM, Installation Effects. For propeller design, accurate predictions of blade forces are crucial for enhancing propeller performance and optimising structural design. Especially to capture complex effects related to propeller-airframe interaction and propellers operating outside normal design conditions CFD simulations are a valuable tool in the propeller design process. While fully resolved blade simulations offer high accuracy, their computational expense restricts their practical use to later design phases, necessitating less expensive methods for early-stage evaluations. The Actuator Line Method (ALM), originally developed by Sørensen and Shen [1], is such a lower fidelity approach where rotor blades are efficiently modelled by imposing body forces into the fluid at discrete actuator points along virtual lines representing the blade axis. The local body forces are determined by sampling local velocity vectors, calculating the local angle of attack and using tabulated airfoil polars. Initially employed for wind turbines in incompressible flows, ALM's proven capability in capturing rotor wakes and tip vortices has encouraged its application to other rotor types, such as helicopter rotors and propellers. Subsequent enhancements, notably the Advanced Actuator Line Method (AALM) by Churchfield et al. [2], have improved the accuracy of velocity sampling and force projection, refining predictions of tip vortices and near-wake phenomena. In previous work [3], we have evaluated approaches of different modelling fidelity implemented in the publically available Simulator for Offshore Wind Farm Applications (SOWFA) OpenFOAM library for simulating incompressible flow around wind turbines, including the ALM approach. In this work, the AALM implementation of SOWFA was extended to compressible propeller flows. Density sampling was introduced into the body-force terms, which were originally derived from incompressibility assumptions. Secondly, these terms are consistently incorporated into both momentum and energy equations. Validation is conducted using a proprietary test rotor geometry provided by Flying Whales alongside experimental benchmark data acquired at Aero Concept Engineering’s wind-tunnel facility in Magny-Cours. The numerical approach employs unsteady RANS (URANS) simulations with the k-omega SST turbulence model. A customised version of the rhoPimpleFoam solver is used, featuring a consistent pressure-velocity coupling to avoid checkerboard pressure fields and increased numerical stability. The aerodynamic polars required as input for the AALM were generated with Upstream CFD’s RANS-based 2D airfoil suite, where values for stall conditions were extrapolated using a flat plate assumption. The computational domain represents the wind tunnel test section with the nacelle geometry mounted on a cylindrical arm. Flow speeds of 0 m/s to 30 m/s and rotation speeds of 6038 rpm to 6078 rpm were simulated, leading to maximum Mach numbers of 0.63 at the propeller tips, thus necessitating a compressible approach. A comprehensive best practice study was conducted to determine efficient settings for solver, numerical schemes, grid resolution and time step size. Table 1: Difference of time-averaged performance quantities to experimental benchmark data �|∞[𝒎/𝒔] Rotor AoA [°] Pitch [°] RPM 𝚫 𝑻 [%] 𝚫𝑸 [%] 𝚫𝜼 [%] 30 0 +15 6078 −2.2 +6.9 −8.5 30 20 +15 6078 +2.4 +6.5 −3.8 20 0 −10 6043 −6.0 −16.8 +13.0 20 20 −10 6038 −13.2 −18.4 +6.3 Validation simulations were conducted for varying inflow speeds, flow angles and blade pitch angles to simulate both forward and reverse thrust conditions. The exact operating conditions are listed in Table 1 together with the percentual differences of thrust 𝑇, torque 𝑄 and the aerodynamic efficiency 𝜂 compared to the experimental benchmark data. Flow visualisations for two representative cases are shown in Figure 1. For forward thrust conditions, thrust predictions were within ±3% and torque values within 6% to 7% of the experimental benchmark, so that the resulting aerodynamic efficiency predicted by CFD was between 3% to 9% lower. Under more challenging reverse thrust conditions, the 97 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria difference between CFD and experiment increased to 6% to 13% for the aerodynamic efficiency. In general, reverse thrust conditions can lead to highly unsteady flow fields and are thus also a challenging application case for experimental test facilities. The developed AALM implementation showed strong numerical robustness for these cases and delivered results within the same magnitude of the experimental benchmark, which is highly encouraging for future studies given the level of modelling fidelity applied in the approach. Future work will include further validation and direct comparison with fully resolved propeller simulations. The application of AALM in combination with scale-resolving methods such as 𝜎-DDES is planned as well as the combination with the EBL turbulence model, which is better suited for flows with strong coherent vortices. Finally, productive use cases of the AALM with simulations involving multiple rotors, such as full-scale airships, drones/eVTOL or wind farms are targeted. Figure 1: Flow visualisation of two representative validation cases. Iso-Q surfaces coloured with velocity magnitude + iso surface of force projection (left) and time-averaged velocity magnitude on z-slice through rotor centre (right). Acknowledgements The authors gratefully acknowledge Flying Whales for funding the project and granting permission to publish the results presented. References [1] J. N. Sørensen and W. Z. Shen, “Numerical modelling of wind turbine wakes”, Journal of Fluids Engineering, vol. 124, no. 2, pp. 393–399, 2002. [2] M. J. Churchfield, S. J. Schreck, L. A. Martínez-Tossas, C. Meneveau, and P. R. Spalart, “A large-eddy simulation of wind-turbine wakes using an advanced actuator-line method”, Journal of Physics: Conference Series, vol. 1037, art. 072032, 2017. [3] N. Sahoo, A. Busse, H. Hetmann, “Insights into OpenFOAM’s rotor simulation and modelling capabilities across different fidelity levels”, Presentation at the 7th French/Belgian OpenFOAM Users Conference, 22nd23rd May 2024, Paris / France. a) Iso-Q surface: forward thrust, 20° b) Mean velocity magnitude: forward thrust, 20° c) Iso-Q surface: reverse thrust, 0° b) Mean velocity magnitude: reverse thrust, 0° 98 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria HYBRID COMPUTATIONAL AEROACOUSTIC SIMULATIONS OF UAV PROPELLER NOISE FRAN PIPUNI ´ C1, LUKA BALATINEC2, and TESSA UROI ´ C3 1Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected].unizg.hr 1Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] 2Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] Keywords: Computational aeroacoustics (CAA), CFD, Propeller noise, Propeller interaction, OpenFOAM Noise prediction for small unmanned aerial vehicles (UAVs) is a critical factor in the design and optimization of modern aerial systems. The increasing prevalence of UAVs in both civil and military applications has highlighted the need for accurate, efficient noise prediction methods. Hybrid methods in computational aeroacoustics (CAA), which decouple the sound generation and sound propagation mechanisms, have emerged as an effective strategy for reducing the computational cost and complexity of noise analysis. Hybrid CAA methods are based on modeling acoustic source terms using results from computational fluid dynamics (CFD) simulations. In this work, a hybrid computational framework is presented for predicting noise generated by a pair of interacting UAV propellers. Theoretical foundations of hybrid aeroacoustic modelling, including the assumptions and approximations involved, are first outlined. The aerodynamic fields used to model the acoustic sources were obtained through transient CFD simulations, with particular attention given to the aerodynamic and acoustic interactions between two rotating propellers. Prior to conducting transient CFD simulations, steady-state simulations using the Multiple Reference Frame (MRF) approach were carried out to provide initial conditions for the unsteady calculations. All simulations were performed using the opensource CFD toolkit foam-extend, a fork of the OpenFOAM environment. The results demonstrate the capability of the presented framework to capture complex propeller interactions and provide detailed input for subsequent acoustic analysis. The developed methodology offers a robust and flexible tool for propeller noise prediction in UAV applications, supporting further research into noise reduction strategies and aerodynamic design optimization. References [1] A. Azeni´ c, “Implementacija i validacija modela propagiranja zvuˇ cnih valova,” 2016. [Online]. Available: https://urn.nsk.hr/urn:nbn:hr: 235:925035 [2] E. Sj¨ oberg, “Implementation of aeroacoustic methods in openfoam,” 2016, unpublished work or internal report. [3] Advanced Precision Composites, “Apc propeller geometry data,” 2024. [4] H. Bu, H. Wu, C. Bertin, Y. Fang, and S. Zhong, “Aerodynamic and acoustic measurements of dual small-scale propellers,” Journal of Sound and Vibration, vol. 511, p. 116330, 2021. 99 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria VALIDATION OF A PRESSURE-BASED COUPLED SOLVER FOR HIGHLY COMPRESSIBLE TURBOMACHINERY FLOWS PAOLO GEREMIA1, GEORGIOS KARPOUZAS2, APOSTOLOS KRASSAS3, SALVATORE RENDA4, THOMAS SCHUMACHER5 1ENGYS S.R.L., [email protected] 2ENGYS Ltd., [email protected] 3ENGYS Ltd., [email protected] 4ENGYS Ltd., [email protected] 5ENGYS GmbH, [email protected]om Keywords: Turbomachinery, compressor, coupled solver, highly compressible flows. The accurate and robust simulation of highly compressible flows in turbomachinery remains a critical challenge in the field of Computational Fluid Dynamics (CFD). In high-speed axial compressors, where Mach numbers frequently exceed unity, numerical instability, convergence difficulties, and sensitivity to discretisation schemes are common issues that compromise both the reliability and predictive capability of CFD models. To address these challenges, coupled solvers are often preferred over segregated solvers due to their improved numerical stability and convergence properties in compressible regimes. Unlike segregated approaches, where pressure and velocity are solved sequentially, coupled solvers resolve all flow variables simultaneously, yielding a tighter and more consistent solution coupling that is particularly advantageous at high Mach numbers. Commercial solvers such as ANSYS CFX, ANSYS Fluent, and NUMECA Fine™/Turbo have set industry standards in this regard, offering robust algorithms and proprietary features that enhance solver stability and performance. These tools often outperform open-source alternatives like OpenFOAM, which, despite their flexibility and extensibility, still face difficulties when applied to high-speed turbomachinery flows due to limitations in solver design and numerical schemes. This work focuses on the validation and application of the pressure-based coupled solver available within the HELYX CFD software [1], which has previously been validated for incompressible and mildly compressible flows. The objective of this study is to evaluate its performance and establish best practices for solving highly compressible turbomachinery problems. Particular emphasis is placed on solver stability, mesh and time step sensitivity, and accurate capturing of transonic flow features such as shock waves and tip-leakage vortices. By applying the solver without structural changes to its underlying numerical algorithms, the study investigates the extent to which an existing industrial-grade, open-source-based coupled solver can handle the stringent demands of high-speed axial compressors. The validation campaign includes two well-established axial compressor benchmarks: the NASA Rotor 67 [2] and the TUDa GLR open stage [3]. These cases feature open-stage configurations with publicly available experimental data, allowing for a detailed and transparent validation process. The validation focuses on replicating key performance parameters such as total pressure ratio and isentropic efficiency across the operating range, while also comparing local flow quantities (e.g., total pressure and velocity profiles) at measurement stations. Results demonstrate that the improved solver maintains stable convergence across the full operating map and achieves good agreement with experimental data, with total pressure ratio and isentropic efficiency errors within 2%, confirming its predictive capability in complex, highly compressible environments as shown in Fig. 1. Figure 1: Total pressure ratio results for the NASA Rotor 67 (left) and TUDa GLR open stage (right) 100 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria In conclusion, this work outlines a set of best-practices for achieving numerical robustness in pressure-based coupled solvers for turbomachinery applications and illustrates their effectiveness through rigorous validation. Future work will focus on expanding the validation suite to include additional compressors and multistage configurations, thereby further strengthening the generality and industrial relevance of the proposed solver framework. Acknowledgments The authors thank all those involved in the organisation of this OpenFOAM Workshop and to all the contributors that will enrich this event. References [1] ENGYS, HELYX - Open-source CFD Software for Enterprise. [2] A. J. Strazisar, J. R.Wood, M. D. Hathaway, K. L. Suder. Laser Anemometer Measurements in a Transonic Axial-Flow Fan Rotor, NASA Technical Paper, 1989. [3] F. Klausmann et al., Transonic compressor Darmstadt - Open test case Introduction of the TUDa open test case, Journal of the Global Power and Propulsion Society, vol. 6, 2022, pp. 318-329, 2022. 101 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NUMERICAL INVESTIGATION OF COAL AND BIOMASS COMBUSTION USING ARRHENIUS-BASED MODELS AND SENSITIVITY ANALYSIS IVAN HEINRICH1, MATEJ ˇ CORAK2, TESSA UROI ´ C3, 1Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected].unizg.hr 2Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] 3Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] Keywords: Combustion, Arrhenius Kinetics, Emission Reduction, Sparse Chemical Kinetics The transition toward more sustainable energy production requires a detailed understanding of how alternative fuel mixtures behave during combustion. In this study, a comprehensive sensitivity analysis examines the effects of fuel composition and key operating parameters on the combustion behaviour of coal and biomass mixtures. The numerical simulations are performed using the open-source computational fluid dynamics (CFD) software OpenFOAM, employing Arrhenius-based reaction [1] models to represent the combustion kinetics [2]. The primary objective of the study is to assess how variations in the biomass-to-coal ratio, reaction kinetics parameters, and operational conditions influence the final composition of combustion products and overall combustion efficiency. By systematically altering these input variables, the sensitivity of the combustion process is quantified, allowing the identification of the most influential factors. This approach provides valuable insight into the interplay between fuel properties and operating conditions, and how they collectively affect flame stability, pollutant formation, and energy release. The combustion model implemented in OpenFOAM relies on simplified chemical kinetics governed by Arrhenius equations, which capture the essential characteristics of the oxidation reactions involved. Specific attention is given to monitoring parameters such as the concentrations of major combustion products (e.g., CO22, CO, andNOx) [3] as well as temperature profiles and burnout efficiency. Mesh independence studies and validation against available experimental data are conducted to ensure the robustness of the simulation results. The outcomes of this sensitivity analysis aim to provide guidelines for optimizing combustion systems utilizing coal and biomass blends, offering potential pathways to enhance performance while minimizing environmental impact. Additionally, the study highlights the flexibility and capability of OpenFOAM for conducting parametric studies and advanced combustion simulations. Future work will involve extending the approach to more complex chemical mechanisms and a broader range of biomass types. References [1] Versteeg, H.K., Malalasekera, W. (2007). An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson Education. (A classic CFD book — very often cited, especially when using OpenFOAM or similar solvers.) [2] Demirbas, A. (2004). Combustion characteristics of different biomass fuels. Progress in Energy and Combustion Science, 30(2), 219–230. (Good for referencing general knowledge about biomass combustion behavior.) [3] Turns, S.R. (2012). An Introduction to Combustion: Concepts and Applications (3rd ed.). McGraw-Hill Education. (This book covers combustion theory, including Arrhenius equations, reaction kinetics, and pollutant formation.) 102 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria TOWARDS A COMPETITIVE OPEN-SOURCE GPU SOLVER FOR OPENFOAM HENDRIK HETMANN1, GREGOR OLENIK2, CHARLES MOCKETT1 1Upstream CFD GmbH, [email protected] 2Technical University of Munich, [email protected] Keywords: GPU computing, High-Performance Computing, OpenFOAM Benchmarking, EXASIM. Over the past decade, progress in High-Performance Computing (HPC) has shifted the state-of-the-art in Computational Fluid Dynamics (CFD) decisively towards GPU-accelerated solvers. Commercial GPU-enabled codes already deliver game-changing performance improvements in simulation time, cost, and energy efficiency [1]. To safeguard the immense strategic value of OpenFOAM as an industrially applicable open-source solution, it must match or exceed the capabilities of commercial software. A high-performance GPU derivative of OpenFOAM is therefore urgently required. The EXASIM project systematically evaluated the performance of OpenFOAM against modern GPU solvers using clearly defined key performance indicators (KPIs) such as the equivalent number of CPU cores giving the same performance as one GPU (“cores per GPU”, CpG), time-to-solution factor (TTSF) and cost-to-solution factor (CTSF). An automated HPC benchmark suite, based on the OpenFOAM Benchmark Runner (OBR), a derivative of the signac framework, was developed to support reproducible, scalable performance tests across diverse hardware platforms. This suite draws on recent best-practice work in benchmarking and software testing strategies, as outlined by Gärtner et al. [2]. Figure 1 shows an excerpt of the analysis of one of the microbenchmarks, where the time spend in different parts of the solution process was investigated. It becomes clear that offloading only the pressure solver to the GPU as done in hybrid approaches like OGL (described in the following) is not sufficient, as for this type of simulation most of the time is spent in the velocity matrix assembly. Figure 1: Excerpt of performance analysis for the Windsor Body, which should be representative for automotive external aerodynamics simulations. Several initiatives to develop GPU-ready OpenFOAM derivatives were reviewed:  OGL (OpenFOAM-Ginkgo-Layer): a hybrid approach offloading linear solvers to GPUs, developed within EXASIM [3][4]  OpenFOAM-UM (and its predecessor prototype zeptoFOAM): simulations run entirely on GPU via a minimalinvasive porting approach, emphasizing quick feature extension and portability [5]  NeoN: a new, generally applicable CFD backend, designed from the ground up for GPU-first architectures and full C++20 compliance [6], linked to OpenFOAM via an interface called FoamAdapter Benchmarking data for Delayed Detached-Eddy Simulations (DDES) of automotive aerodynamics are summarised in Table 1 (missing data will be assessed in time for the workshop). A clear benefit of executing the entire simulation on GPU (zeptoFOAM) vs. offloading only the linear solver to GPU (OGL) is seen, as expected. Initial data published for OpenFOAM-UM [5] states a CpG value of around 220, indicating an appreciable performance improvement over zeptoFOAM. However, analysis of performance reported for a comparable commercial solver suggests there is much a) Visualisation of Windsor Body microbenchmark in weak-scaling study mode b) Share of wall-clock time per timestep used by different stages of the linear solution process, comparing OpenFOAM running on CPU only and with the pressure equation solved on GPU using OGL. MA: matrix assembly, U: velocity, P: pressure, OH: overhead (additional solver steps e.g. turbulence equations) 103 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria more potential (approx. 3x), which is currently untapped by OpenFOAM implementations. To achieve this, it is believed that the core numerical methods must be structured from the ground up to fully exploit GPU architectures. This is the objective of the NeoN initiative, which targets the game-changing levels of performance demonstrated by the commercial reference code. As such, each OpenFOAM initiative represents a different balance between risk and reward: from the practical benefits of minimal-invasive porting to the performance targeted by revolutionary redesign. Table 1: Key performance indicators for different GPU solver approaches using NVIDIA A100 GPUs for automotive aerodynamics DDES. Higher numbers are better for all KPIs. Solver CpG CPU cores required to give the same performance as 1 GPU TTSF Time-to-solution reduction factor CTSF Cost-to-solution reduction factor OGL 15.9-19.5 0.15-1.03 0.42-0.51 zeptoFOAM 105-140 0.55-1.13 1.22-1.64 OpenFOAM-UM to be evaluated to be evaluated to be evaluated NeoN to be evaluated to be evaluated to be evaluated Commercial ref. 680 18 6.3 This contribution will summarize the EXASIM findings, highlight the OBR benchmark framework’s capabilities, and compare the different GPU initiatives currently shaping the future of OpenFOAM. Acknowledgements The work presented in this talk has been executed within the EXASIM project. The project has been funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor European commission can be held responsible for them. References [1] H. Owen, “AutoCFD HPC”, Presentation at the 4th Automotive CFD Prediction Workshop, Belfast, UK, 2024. Available online: https://autocfd4.s3.eu-west-1.amazonaws.com/presentations/HPC-Website/HPCOwen.pdf [2] W. Gärtner, G. Olenik, M. E. Fadeli, L. Petermann, A. Kronenburg, H. Marschall, and H. Anzt, "Testing Strategies for OpenFOAM Projects", OpenFOAM Journal, Vol. 5, pp. 115–130, 2025. DOI: https://doi.org/10.51560/ofj.v5.134. [3] Olenik, G., Koch, M., Boutanios, Z. et al. Towards a platform-portable linear algebra backend for OpenFOAM. Meccanica (2024). https://doi.org/10.1007/s11012-024-01806-1 [4] Hartwig Anzt, Terry Cojean, Goran Flegar, Fritz Göbel, Thomas Grützmacher, Pratik Nayak, Tobias Ribizel, Yuhsiang Mike Tsai, and Enrique S. Quintana-Ortí. 2022. Ginkgo: A Modern Linear Operator Algebra Framework for High Performance Computing. ACM Trans. Math. Softw. 48, 1, Article 2 (March 2022), 33 pages. https://doi.org/10.1145/3480935 [5] CINECA and EuroCC Italy, "OpenFOAM-UM: GPU acceleration in OpenFOAM via minimal-invasive porting," EuroCC White Paper, April 2025. Available online: https://euroccitaly.it/wpcontent/uploads/2025/04/EuroCC_White-Paper_eng.pdf [6] https://github.com/exasim-project/NeoN, accessed 28.04.2025 104 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria USING EXPRESSION TEMPLATES FOR MEMORY LOCALITY ILYA POPOV1, DMITRY ASTAKHOV2 1ISTEQ BV, ilya.pop[email protected] 2ISTEQ BV, astakho[email protected] Keywords: Performance, Memory Locality, Expression Templates, Bandwidth In this work we present further development of the expression templates technique for memory locality optimization in OpenFOAM-based CFD codes. This work consists of two parts, a general purpose expression templates library ET and its adaptation for OpenFOAM called FVE. ET library follows general design of the Boost.YAP expression templates library. The idea is to split working with expression templates into spearate phases of expression construction, transformation, and evaluation. An expression is a simple struct that does not know anything about how it is about to be used. template <typename Op, typename . . . Args> st r uct expr ; Here, Op is an operation (a function object to be applied to the arguments), and Args... are arguments, which can be expressions themselves or terminals (leaf nodes, values or references to values, such as scalars or fields). Therefore all kinds of expressions are instatiations of the same template with different template arguments. This is in contrast with alternative designs where each operator has it’s own type. The advantage is that where is much less duplication of code, and much easier to handle all kinds of operators in uniform way. For example, expression auto e1 = et : : make_terminal (3) + 5; has type expr<std : : plus <>, int , int> Note that at this point the expression is purely syntactic, similar to an AST (abstract syntax tree) used in compilers. Currently, ET library supports all C++ operators (+,−,∗,/,%, unary −,&&,||,!,&,|, , , ==,! =,<,<=,>, >=) plus identity function and select (analog of ?: operator in C++). Once an expression is constructed, we can manipulate it in various ways. For this, ET library provides utility fucntions et : : transform_matching ( expr , transform ) that call a function object transform on the nodes of the graph that match the signature of the transform. The transform would typically return a transformed copy of the expression. (Expressions themselves are expected to be lighweight and passed by value). Another utility function et : : transform_terminal ( expr , transform ) applies a transform to all terminals. Note that matching and manuipulating trees of expression is similar to how compilers (such as GCC or Clang) operate. This is also similar to how the AI frameworks (such as TensorFlow, Torch, etc.) operate. Finally, we can evaluate an expression using evaluate(expr) function, that just applies each nodes’s operation Op to its arguments. Additional utilities provided by ET include pretty-printing expressions, drawing them as graphs in graphviz format, using std::placeholders, and computing derivatives. FVE uses facilities provided by ET to manipulate expressions found in OpenFOAM-based codes. First, it defines a number of new expression operations, associated with OpenFOAM’s elementwise functions and with fvc:: functions. It also associates a location (either cells or faces) with every expression. Additionally, it allows getting 105 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria using single-step, Arrhenius-type combustion kinetics regulated by the combustion progress variable Z, following prior work [1]. Because the inlet supplies premixed, unburnt reactants, Z= 0 at the inlet. The fluid is treated as a calorically perfect gas with constant specific heat ratio γand a (non-dimensional) ideal gas constant of 1. The remaining model constants - the activation energy Ea, heat release rate, C-J post-shock (von Neumann) temperature Tign, and Damk¨ ohler number Da - are further discussed in our prior work [3]. Note that the flow time scale is t0≡y0/u0. This is the fixed convective time scale resulting from the velocity and length scales used for non-dimensionalization (near the sound speed and detonation height, respectively). The time evolution of the system was computed by numerically solving the aforementioned system of equations. This was performed using a custom application based on blastFoam, an open-source solver based on OpenFOAM-9 and developed by Synthetik Applied Technologies [4]. This solver has been extensively validated for simulating highly compressible, reacting flows, particularly those with detonations. The variables were estimated at the cell faces using the quadratic MUSCL scheme in blastFoam, which uses the gradients and Hessians (that is, the 1st and 2nd spatial derivatives) to perform high order reconstruction. The convective fluxes are estimated using the HLLC scheme for shock capturing. Temporal integration was performed using the 4th order, strong-stability-preserving Runge-Kutta (RK4SSP) scheme. An adaptive time step was used to fix the CFL at 0.75. To directly compute steady attractors and construct the aformentioned bifurcation diagram we developed the open-source software tool “LOCAFOAM.” This code combines the accessibility and widespread familiarity of OpenFOAM with the power of the numerical continuation package LOCA (“Library of Continuation Algorithms”), part of the scientific computing suite Trilinos developed by Sandia National Laboratories. Besides having a long and successful history of use for fluid mechanical problems, LOCA has the added advantage of a highly modular structure; all of its algorithms operate on a discretized dynamical system through abstract “model interfaces.” In other words, one can set up and discretize a system in OpenFOAM without ever interacting with LOCA or other Trilinos software. As OpenFOAM is exceedingly easy to use and widely popular among combustion researchers (and Trilinos far less so), this pairing is highly advantageous. Once a system is discretized in OpenFOAM, it is manipulated by LOCA through Trilinos’ own linear algebra data structures within an abstract “model evaluator” in the Trilinos package Thyra. The user has the option of either hand-coding a Jacobian or using an approximate one estimated via finite differencing. The latter enables greater abstraction and ease of use but can be computationally expensive. To mitigate this downside, we employ a graph coloring approach that takes advantage of the structural orthogonality of the Jacobian to drastically reduce the number of function evaluations needed for its estimation. All aspects of the continuation problem are subsequently handled by LOCA, including the computation of solutions and the eigenvalue problem for stability and bifurcation analysis. Although we developed LOCAFOAM specifically for this study, we emphasize that it is quite general and can be used for any PDE continuation problem. Acknowledgments This research was partially supported by the Air Force Office of Scientific Research under award number FA9550-23-1-0222 (contract monitor Dr. Chiping Li). The authors are especially thankful to Dr. Eric Phipps and Dr. Andrew Salinger for their help using Trilinos software. References [1] D. E. Paxson, “Computational assessment of the impact of wave count on rotating detonation engine performance,” AIAA SciTech Forum 2023, 2023. [2] J. Koch and J. N. Kutz, “Modeling thermodynamic trends of rotating detonation engines,” Physics of Fluids, vol. 32, no. 126102, 2020. [3] T. Kickliter, V. Acharya, E. Young, and T. Lieuwen, “Asymptotic dynamics of rotating detonation engines,” AIAA SciTech Forum 2025, 2025. [4] S. A. T. LLC., “blastFoam: A solver for compressible multi-fluid flow with application to high-explosive detonation,” 2020. [Online]. Available: https://github.com/synthetik-technologies/blastfoam 112 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria A NUMERICAL FRAMEWORK FOR WEAR ANALYSIS OF ROUGH SURFACE CONTACTS ACROSS LUBRICATION REGIMES LUKA BALATINEC1, TESSA UROI ´ C2and HRVOJE JASAK3 1Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] 2Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, [email protected] 3University of Cambridge, Department of Physics, The Cavendish Laboratory, [email protected] Keywords: Numerical wear analysis, Lubricated contacts, Rough Surfaces, Archard model, Finite Area Method, foam-extend The presence of wear in lubricated tribological systems, and its effect on mechanical component performance, is welldocumented. Traditional approaches to studying wear often rely on experimental methods, which, although effective, are costly and complex compared to numerical alternatives. Numerical simulations offer enhanced efficiency and provide deeper insights into wear mechanisms under various lubrication regimes. Archard [1, 2] defined wear as the progressive loss of material due to mechanical interaction at surface asperities during relative motion, typically sliding or rolling. Adhesive wear is one of the predominant forms of wear in mechanical systems [3] and, as such, is the focus of this research. The Archard Wear Model, widely recognized for its general applicability [4], serves as the foundation for the wear model used in this work. In this research, a numerical framework based on the Finite Area Method (FAM) [5] and implemented in foam-extend is presented for simulating wear of rough surfaces under varying lubrication regimes. The model couples a wear algorithm based on Archard’s law with a deterministic asperity contact model, enabling the use of actual measured surface profiles. The lubrication effects are captured by solving a modified Reynolds equation with cavitation modelling [6], formulated within the Finite Area framework. Surface evolution due to wear is modeled by iteratively updating the surface geometry based on contact and lubrication parameters. Validation was conducted using several tribological test cases, including Pin–On–Disc and Ring–On–Block configurations. Numerical results demonstrated excellent agreement with data from the literature in terms of contact pressure distribution, wear depth, and surface evolution. Additional test cases, such as Ball–On–Flat and Ring–On–Ring setups using real surface scans, further confirmed the ability of the model to predict wear progression under both dry and lubricated contact conditions. Finally, the framework was applied to simulate wear in a lubricated Ball–On–Disc configuration using Shell Turbo T68 oil, accurately capturing lubrication regime transitions from mixed to near-boundary lubrication. The developed wear framework shows robust performance across various tribological scenarios, supporting its applicability for advanced wear prediction and surface evolution studies in lubricated contact problems. Acknowledgments This work was supported by the Croatian Science Foundation (project number DOK-2020-01). References [1] J. F. Archard, “Contact and Rubbing of Flat Surfaces,” Journal of Applied Physics, vol. 24, no. 8, pp. 981–988, Aug. 1953. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.1721448 [2] J. F. Archard and W. Hirst, “The Wear of Metals under Unlubricated Conditions,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 236, no. 1206, pp. 397–410, 1956. [Online]. Available: http://www.jstor.org/stable/99967 [3] M. M. Khonsari and E. R. Booser, Applied tribology: bearing design and lubrication, third edition ed. Hoboken, NJ: John Wiley & Sons Inc, 2017. [4] H. C. Meng and K. C. Ludema, “Wear models and predictive equations: their form and content,” p. 15. [5] ˇ Z. Tukovi´ c, “The finite volume method on domains of changeable shape,” Ph.D. dissertation, 2005. [6] V. ˇ Skuri´ c, “Numerical simulation of lubricated wire rolling and drawing,” Ph.D. dissertation, 2019. 113 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Generative AI for Extreme Ocean Waves for Ship CFD Kevin Maki1, Benedetto Di Paolo2, Adam Sak1, Eric Heilshorn1, Paolo Geremia2 1University of Michigan, Department of Naval Architecture and Marine Engineering, USA, [email protected] 2Engys, S.R.L., Trieste, Italy Keywords: Generative AI, extreme waves, transient maneuvering The next-generation ocean-going vessels will rely increasingly on alternative energy sources for propulsion, and high-fidelity simulation is essential to design platforms that can operate efficiently and safely. Traditional methods—model tests, nonlinear potential-flow simulations [1], or idealized design-wave studies[2] can only partially capture the most complex wind-waveship interaction. Although CFD has advanced steadily over the past several decades [3], it is still too expensive to use for Monte-Carlo type simultion of long-time exposure in the ocean. To address this, we introduce GenWave [4], a generative-AI approach that learns from random JONSWAP seas and then synthesizes focused wave fields likely to produce extreme responses. Trained in a 30-minute Colab session on TPUs, GenWave generates unlimited unique seaways whose peak elevations exceed four times the standard deviation. We apply GenWave to the ONR Tumblehome in sea state 6 (Hs= 7 m, Tp= 15 s), using HELYX-Marine with adaptive mesh refinement on the air-water interface (3.9M cells, five prism layers, 0.625 mm propeller resolution) and a sliding-grid propeller model. GenWave is used to bootstrap a sequence of language models, and 10 uniqued focused wave cases, alongside a calm-water baseline, are simulated at model scale 1-to-32. The unique wave environments are used to assess the stopping maneuver in following seas. The time history of the ten wave fields is shown in Fig. 1. Figure 1: Ten unique design wave time histories at the location of encounter with the ship. An image of the ship during backing after the large wave passes is shown in Fig. 2. The ship is already backing, and the ventilated propeller is shown pushing air-water mixture towards the bow. Figure 2: Ship as large wave passes stern and causes propeller ventilation Results show that while calm-water thrust peaks more rapidly, wave-induced torque can be 15 times higher (6,440 kN·m vs. 419 kN·m) with a 609 kN·m standard deviation. Trajectory scatter reveals complex propeller–wave interactions, including intermittent emergence and ventilation (Fig. 3). This demonstrates GenWave’s potential to uncover extreme loads that traditional methods may miss. 114 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 3: Visualization of the propeller in crashback shown by a contour of the Qcriterion Acknowledgments The authors thank all those involved in the organisation of this OpenFOAM Workshop and to all the contributors that will enrich this event. References [1] K. Weems and V. Belenky, “Reduced-order model for ship motions incorporating a volume-based calculation of bodynonlinear hydrostatic and froude-krylov forces,” Ocean Engineering, vol. 289, p. 116214, 2023. [2] P. A. Anastopoulos and K. J. Spyrou, “Determining ship manoeuvrability failures in extreme waves using the critical wave groups concept,” 2024. [3] K. M. Silva and K. J. Maki, “Towards a computational fluid dynamics implementation of the critical wave groups method,” Ocean Engineering, vol. 235, p. 109451, 2021. [4] W. Xu, “A machine learning framework to model extreme events for nonlinear marine dynamics,” Ph.D. dissertation, University of Michigan, June 2020. 115 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria TOWARD FIRE DYNAMICS SIMULATIONS OF THERMAL RUNAWAY PROPAGATION IN BATTERY ENERGY STORAGE SYSTEMS DANYAL MOHADDES1, YI WANG2, 1FM, Research Division, [email protected] 2FM, Research Division, [email protected] Keywords: Lithium-ion batteries, Thermal runaway, Fire dynamics Introduction Battery Energy Storage Systems (BESS) are a prevalent and growing means of stabilizing power grids with high penetration of intermittent renewable energy sources [1]. BESS typically rely on lithium-ion batteries (LIBs) as their basic building block. LIBs are packed together in an enclosure to form a module. Multiple modules are assembled into a rack, and these racks are housed within a container to form one of many units in a BESS. LIBs pose a safety hazard if abused thermally, mechanically or electrically due to a failure mechanism known as thermal runaway (TR) [2]. TR is a localized breakdown of the LIB internal structure, causing exothermic chemical reactions that lead to a propagating reacting front within the LIB. It further results in the production of large quantities of flammable gas, which increase the LIB’s internal pressure and lead to its rupture, releasing the flammable gas.This rapid self-heating thermally abuses neighboring LIBs within a module, driving them into TR, initiating a cell-to-cell cascade known as TR propagation (TRP). Since the LIBs are enclosed within a module, the large volume of vented gas produced by TRP is driven out of vents and openings on the module body [3], whereupon it mixes with ambient air and poses a fire and explosion hazard. Fire dynamics thus come into play at the module scale, due to the potential for flame heating of the module exterior by a vented gas-fueled fire. To this end, the present study will discuss recent model developments and their implementation in an OpenFOAM solver at the module scale, toward building a predictive capability for analyzing and predicting rack and unit-scale TRP events and associated fire dynamics. Methodology The developments we consider in this work utilize FireFOAM, an LES solver in the OpenFOAM family developed originally at FM [4] for simulation of fire dynamics and fire suppression. FireFOAM is capable of simulating complex fire scenarios, handling the coupled physics arising from gaseous flows and combustion, water sprays, films, radiation and solid-phase pyrolysis. FM maintains a branch of FireFOAM with specific capabilities for handling practical fire scenarios [5]. In this work, we implemented TRP as a region model derived from the base pyrolysis model within FireFOAM. In particular, our implementation targets recent module-scale [3] and rack-scale [6] TRP fire dynamics experiments conducted at FM using large pouch-format Li-NMC532 (nickel-manganese-cobalt oxide)-graphite cells. The mathematical formulation of the TRP model is based on our earlier work [7] and is only described at a high level here. The thermo-chemical behavior of each LIB is described by one-dimensional reaction-diffusion equations for mass, energy and species to consider multispecies reactions and mass loss due to venting. A one-dimensional treatment is reasonable due to the high in-plane thermal conductivity of pouch-format LIBs, and thermal and kinetic parameters are obtained via parameter optimization against singlecell experiments [7]. To analyze TRP, we apply an inter-cell thermal resistance as a boundary condition between LIBs, the value for which we infer from the module-scale experimental data [3], and assemble a stack of LIBs as a single object representing a module within the region model. We include treatments for the coupled boundary condition such that vented gas exits the module at physically realistic locations based on the module geometry. The experiments use a torch to ensure the vented gas burns to permit studying of the fire dynamics, and thus we likewise assume that all vented gas burns in the ambient air. Results Here we show a sample of the results to be presented at the Workshop. The geometry considered is based on twomodule TRP tests conducted at FM [3] with twelve LIBs per module. Each module is a rectangular prism with dimensions 116 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 1: Left: Snapshot at a time t2< t < t3, showing the stoichiometric isosurface; Center: Space-time evolution of temperature in Module 0, with white lines showing boundaries between LIBs; Right: Coupled boundary temperature over time. 0.85 m×0.45 m×0.11 m. The lower module, denoted Module 0, is located 1.1m above the ground, and Module 1 is located 0.055 m above Module 0. The simulation begins by heating the non-coupled boundary of Module 0 to TR, which induces TRP at t=t0, where times are marked by vertical lines in the right subfigure. The left subfigure shows an instantaneous snapshot from the simulation, showing that the vented gas burns in two turbulent, buoyant plumes due to the presence of vent holes on either side of each module. The fire imposes a heat flux on the coupled boundary of both modules, increasing their temperature, as shown in the right subfigure. This initiates a second TRP front in Module 0 at t1and a third in Module 1 at t2. The bi-directional progression of TRP in Module 0 is shown in the center subfigure. TRP completes in Module 0 and Module 1 at t3and t4, respectively. The fire is thus largest for t2< t < t3, when there are three TRP fronts propagating simultaneously; the snapshot shown in the left subfigure corresponds to this time period. These results are phenomenologically in agreement with the experimental results [3], and further analysis of the simulation results, along with comparisons to experimental data, will be provided during the Workshop presentation. Acknowledgments The authors acknowledge the invaluable insights provided by the scientists at FM Research, in particular Drs. J. Cuevas, A. Krisman, X. Lu, N. Ren, T. Xiao, G. Xiong, D. Zeng and Mr. W. Brown, and the management support of Dr. S. Becz. Funding was provided by FM Research in the context of the Strategic Program for Battery Safety. References [1] Y. Yang, S. Bremner, C. Menictas, and M. Kay, “Modelling and optimal energy management for battery energy storage systems in renewable energy systems: A review,” Renew. Sust. Energ. Rev., vol. 167, p. 112671, 2022. [2] Q. Wang, B. Mao, S. I. Stoliarov, and J. Sun, “A review of lithium ion battery failure mechanisms and fire prevention strategies,” Prog. Energ. Combust., vol. 73, pp. 95–131, 2019. [3] L. Gagnon, D. Zeng, R. S. Barlow, and Y. Wang, “Detailed measurements of thermal runaway propagation and fire spread for li-ion battery multi-module setup,” J. Energ. Storage, submitted. [4] Y. Wang, P. Chatterjee, and J. L. de Ris, “Large eddy simulation of fire plumes,” Proc. Combust. Inst., vol. 33, pp. 2473–2480, 2011. [5] FM Research, “FireFOAM.” [Online]. Available: https://github.com/fmglobal/fireFoam [6] J. Cuevas, D. Zeng, and Y. Wang, “Insights on thermal runaway and fire propagation in a lithium-ion battery energy storage system,” 11th International Seminar on Fire and Explosion Hazards, 2024, submitted. [7] D. Zeng, D. Mohaddes, L. Gagnon, and Y. Wang, “Modeling initiation and propagation of thermal runaway in pouch li-ion battery cells: Effects of heating rate and state-of-charge,” Proc. Combust. Inst., vol. 40, no. 1, p. 105316, 2024. 117 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Numerical Investigation of Heat Transfer in Nanoparticle-Enhanced Taylor Flow Using thermalMPPICInterIsoFOAM Fahad S. Al-Gburi1, Faysal Khaleel1,Yulong Ding2, Jason Stofford1 1School of Engineering, University of Birmingham, Birmingham B15 2TT, United Kingdom, [email protected], [email protected], j.staffor[email protected] 2School of Chemical Engineering, University of Birmingham, Birmingham, B15 2TT, United Kingdom, y[email protected] Keywords: Taylor flows, Heat Transfer Enhancement, Microreactors, MPPIC, Volume of Fluid, Interfacial force, Iso advection This study presents a numerical investigation of heat transfer in three-phase Taylor flows containing nanoparticles within millimetric channels, using a newly developed OpenFOAM solver, thermalMPPICInterIsoFOAM. This solver extends MPPICInterFOAM by introducing a thermal energy equation and enabling dynamic coupling of fluid, bubble, and solid phases under laminar flow. Modifications include: (1) full incorporation of temperature-dependent thermophysical properties; (2) additional thermophoretic and Brownian force models acting on the Lagrangian phase; (3) replacement of the MULES interface capturing method with the geometrically conservative isoAdvection scheme [1] to suppress unphysical currents at phase boundaries; and (4) correction of the interface force submodel to eliminate spurious particle behavior near the gas-liquid interface. The improved force formulation is confined to the interfacial region using a gradient threshold filter on ∇α, and the interface force coefficient is dynamically scaled with particle properties, as detailed in [2]. The accuracy and stability of the thermalMPPICInterIsoFOAM solver were validated by comparing results with prior studies on nanoparticle-based fluid flows [3] and interface dynamics benchmarks [4]. This framework offers a scalable and extensible tool for studying nanoparticle-enhanced cooling strategies in microreactors, electronics, and biomedical systems [5, 6] Simulations were performed at Re = 200 in a capillary channel using two configurations: (1) two-phase Taylor flow (water and air), and (2) three-phase flow with suspended 100 nm Al2O3nanoparticles at a volume fraction of ϕ= 1 ×10−4. The solver computes the thermal field, particle transport, and interface evolution using a coupled VOF–MPPIC framework with iterative property updates at each time step. Heat transfer was analyzed using local and average Nusselt numbers across the unit cell. Figure 1: Temperature distribution (isotherms) along the axial direction of the channel for two configurations at Reynolds number Re = 200: (Upper) Two-phase Taylor flow without nanoparticles; (Lower) Multiphase flow with suspended nanoparticles (volume fraction ϕ= 1 ×10−4, particle diameter dp= 100 nm). Results show that nanoparticles increase the liquid film thickness between the bubble and the channel wall (Figure 1), reducing wall shear and accelerating bubble movement. This change results from the altered bubble cross-section and lower 118 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria flow resistance due to particle-induced film thickening. As the film becomes thicker, it buffers sharp thermal gradients typically observed near the wall, enhancing thermal stability. This effect is reflected in the Nusselt number distributions (Figure 2). In the two-phase case (without nanoparticles), the local Nusselt number shows a steep variation near the wall at the transition from the bubble to the liquid slug: a sudden drop followed by a sharp rise. This behavior is attributed to the very thin liquid film that forms in the two-phase case, allowing heat to transfer inefficiently from the wall to the low-conductivity gas phase. Once the flow transitions into the liquid slug, heat is transferred into the water, resulting in a rapid increase in the Nusselt number. In contrast, in the nanoparticle-enhanced case, the thicker liquid film reduces direct wall-to-gas heat transfer and moderates the thermal transition into the slug, resulting in a smoother and more gradual Nusselt number profile. Although particle accumulation near the wall may slightly reduce local peaks due to increased viscosity and wall friction, the overall heat transfer performance improves. The smoothed and extended heat transfer profile enhances the average Nusselt number by approximately 6% at Re = 200 and ϕ= 1 ×10−4, without disrupting the beneficial internal recirculations. Figure 2: Local Nusselt number distribution as a function of normalized axial position X/Lc, where Xis the channel length and Lc represents the unit cell length (slug + bubble). Acknowledgments This research was supported by the Higher Committee for Education Development in Iraq (HCED) Scholarship program, a Royal Academy of Engineering/Leverhulme Trust Research Fellowship [LTRF2122-18-108], and an EPSRC research grant [EP/Z002540/1]. References [1] J. Roenby, H. Bredmose, and H. Jasak, “A computational method for sharp interface advection,” Royal Society Open Science, vol. 3, no. 11, p. 160405, 2016. [2] F. A. Khaleel, F. S. Al-Gburi, Y. Ding, and J. Stafford, “Modelling three-phase taylor flows using a combined mppic-vof approach,” International Journal of Multiphase Flow, vol. 184, p. 105117, 2025. [3] J. Zhang, Y. Diao, Y. Zhao, and Y. Zhang, “Thermal-hydraulic performance of sic-water and al2o3-water nanofluids in the minichannel,” Journal of Heat Transfer, vol. 138, no. 2, p. 021705, 10 2015. [Online]. Available: https://doi.org/10.1115/1.4031699 [4] Y. Liao, Q. Wang, U. Caliskan, and S. Miskovic, “Investigation of particle effects on bubble coalescence in slurry with a chimera MP-PIC and VOF coupled method,” Chemical Engineering Science, vol. 265, p. 118174, 1 2023. [5] A. R. Betz and D. Attinger, “Can segmented flow enhance heat transfer in microchannel heat sinks?” International Journal of Heat and Mass Transfer, vol. 53, no. 19, pp. 3683–3691, 2010. [6] B. Mehta, D. Subhedar, H. Panchal, and Z. Said, “Synthesis, stability, thermophysical properties and heat transfer applications of nanofluid – a review,” Journal of Molecular Liquids, vol. 364, p. 120034, 2022. 119 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NUMERICAL MODELING OF ALUMINUM ANNEALING IN INDUSTRIAL FURNACES USING OPENFOAM NIKOLAOS CHAMAKOS 1*, GEORGE PASHOS 1, DIMITRIOS KORTSELIS 2, DIMITRIOS TSITSIRIDAKIS 2, ANDREAS MAVROUDIS 2, IOANNIS CONTOPOULOS 1, 3 1ELKEME SA, 61st km AthensLamia National Road, GR 32011, Oinofyta, Viotia, Greece 2ElvalHalcor SA, 61st km AthensLamia National Road, GR 32011, Oinofyta, Viotia, Greece 3Research Center for Astronomy and Applied Mathematics, Academy of Athens, GR 11527 Athens, Greece * Corresponding author, E-mail: [email protected] Keywords: Computational Fluid Dynamics, Heat Transfer Simulation, Annealing Process Optimization The annealing process is a crucial step in aluminum production, ensuring the desired mechanical and metallurgical properties of the final product [1]. Achieving uniform temperature distribution during annealing is essential to optimize product quality and minimize energy consumption. In this work, a numerical model was developed to simulate the temperature evolution of stacked aluminum plates during the annealing process in an industrial furnace. The model was implemented in OpenFOAM, using the chtMultiRegionFoam solver. This solver enables the coupling of heat transfer in both solid and fluid domains, allowing for an accurate representation of conjugate heat transfer mechanisms. The study focuses on an industrial annealing furnace consisting of multiple heating zones and an air recirculation system driven by dedicated fans. The computational model incorporates key furnace design features, including the influence of air recirculation, steel support frames, and thermal anisotropy caused by air micro-gaps in the stacked aluminum plates. The computational mesh was generated using snappyHexMesh, while simulations were performed using parallel computing to handle the complex flow and heat transfer calculations. Model validation was performed by comparing numerical predictions with experimental measurements from thermocouples placed at different positions within the aluminum stacks. The results demonstrated a strong agreement between simulated and measured temperatures, confirming the model’s accuracy. Beyond the validation, the model can also investigate different furnace geometries and loading configurations to evaluate their effect on heat distribution and process efficiency. By simulating various scenarios, it becomes possible to optimize furnace operation, reducing cycle times and energy costs while maintaining high-quality aluminum products. Figure 1 illustrates the airflow patterns within the furnace, showcasing the recirculation behavior that influences temperature distribution. This visualization highlights the complex interactions between airflow and stacked aluminum plates, reinforcing the importance of accurate modeling for optimizing furnace performance. Figure 1: Airflow patterns inside the annealing furnace, illustrating the recirculation behavior around the stacked aluminum plates. The recirculating flow significantly affects heat distribution and temperature uniformity during the annealing process. Overall, this work demonstrates the potential of OpenFOAM in enhancing industrial annealing processes through advanced CFD modeling. By utilizing open-source simulation tools, manufacturers can achieve improved process control and better energy efficiency highlighting the role of computational modeling in modern industrial applications. References [1] G. E. Totten and D. S. MacKenzie, eds., Handbook of Aluminum: Vol. 1. Physical Metallurgy and Processes, CRC Press, 2003. 120 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Progress in Design and Analysis of Fusion Power Plant Breeder Blanket Systems using OpenFOAM Aleksander Dubas1, Benedict Smith2, Pranav Naduvakkate2, Rupert Eardley-Brunt2, Gerasimos Politis2, Andrew Davis2 1UK Atomic Energy Authority, aleksander[email protected] 2UK Atomic Energy Authority Keywords: Fusion, Blankets, Multiphase, MHD, Thermal-Hydraulics Fusion power plants are a potential key technology on the road to achieving a net zero, sustainable future. However, they must be commercially viable and a significant engineering effort is essential to realise the potential of fusion energy. One of the most vital components to a successful magnetic confinement fusion concept is the breeder blanket. This component is multi-purpose and needs to be able to cope with the heat loads of close plasma proximity and those generated by fusion neutrons, as well as generating a sufficient amount of tritium to complete the fuel cycle and do all the above with maximal efficiency. The design and analysis of breeder blankets is a challenging multiphysics problem. From a fluid dynamics perspective they include a high heat flux thermal-hydraulics problem, which might include phase changes, coupled bi-directionally to the neutronics through doppler broadening. Additionally, where liquid metals are used as a breeder or multiplier material, magnetohydrodynamic effects are also significant due to the strong magnetic field and flowing conductive fluid. To successfully model such a coupled system computationally, a scalable solver is required. As part of our wide ranging digital engineering program at UKAEA, we use OpenFOAM as a scalable and well established finite volume code. The open source nature of the framework and file formats enables easy coupling into our high performance engineering software ecosystem. In particular here we present our progress on advanced heat transfer modelling and magnetohydrodynamics using publicly available OpenFOAM solvers. For a heat transfer problem representative of the high heat flux conditions experienced by breeder blankets in proximity to the plasma, we present our progress on modelling multiphase flow through a Hypervapotron[1]. As a design that leverages phase changes this presents an excellent target for both development and demonstration of capability. Here we incrementally build complexity, starting with fluid flow, adding heat transport and finally inclusion of multiphase boiling models. In our investigation of OpenFOAM’s magnetohydrodynamics capability, we look first at verification and validation against the known solutions of Shercliff and Hunt. This is then extended to comparison with experimental magnetohydrodynamics rigs at UKAEA, including a cross-comparison and benchmarking against other available codes. Our investigation and development of OpenFOAM has found it to be scalable and suitable for use on our high performance computing clusters[2]. This has significant benefits in enabling engineering analysis in greater detail, higher accuracy and physical complexity. Not only is it possible to scale up, it is also trivial to scale out design simulations onto any available hardware, allowing for simulation driven design that fully incorporates complexity and uncertainty quantification from the outset. Acknowledgments This work has been funded by the Fusion Futures Programme. As announced by the UK Government in October 2023, Fusion Futures aims to provide holistic support for the development of the fusion sector. The authors thank all those involved in the organisation of this OpenFOAM Workshop and to all the contributors that will enrich this event. References [1] J. Milnes, A. Burns, and D. Drikakis, “Computational modelling of the hypervapotron cooling technique,” Fusion Engineering and Design, vol. 87, no. 9, pp. 1647–1661, 2012. [Online]. Available: https://www.sciencedirect.com/ science/article/pii/S0920379612003171 [2] R. W. Eardley-Brunt, A. J. Dubas, and A. Davis, “On scalable liquid-metal mhd solvers for fusion breeder blanket multiphysics applications,” Plasma Physics and Controlled Fusion, vol. 66, no. 1, p. 015015, dec 2023. [Online]. Available: https://dx.doi.org/10.1088/1361-6587/ad100a 121 The 20th OpenFOAM Workshop (OFW20), 30 June - 4 July 2025, Vienna, Austria ADVANCES IN LIQUID FILM SIMULATIONS USING GEOMETRIC VOLUME-OF-FLUID TIAGO CAMPOS1,2, GEORG BR ¨ OSIGKE2, LIONEL GAMET1, OLIVIER LAGET1, VERONIQUE PENIN1, PASCAL ALIX1, JENS-UWE REPKE2 1IFP Energies nouvelles, Corresponding author: [email protected] 2Technical University of Berlin, Process Dynamics and Operations Group Keywords: Volume-of-Fluid, isoAdvector PLIC-RDF, dynamic contact angle, partial slip, Navier slip In chemical engineering, liquid films are often encountered in many process equipments like structured packings or fixed beds of catalysts, either operated in ascendant flow or as trickle beds. The interaction between liquid films and gases is thus of the utmost relevance for many industrial processes, particularly regarding mass transfers between phases. To better understand these complex phenomena, this research uses Computational Fluid Dynamics (CFD), with a focus on the Volume-of-Fluid (VOF) method. The VOF method inherently captures the interface between the phases. The geometric VOF method named isoAdvector developed by Roenby et al. [1] is employed to resolve the interface. The use of the Reconstructed Distance Function (RDF) was implemented by Scheufler and Roenby [2] to improve interface reconstruction and the computation of the interface curvature used for surface tension calculations. This method is reported to reduce spurious currents by orders of magnitude when compared to traditional algebraic interface compression methods [3]. A proper modelization of the moving contact line between the liquid, gas, and solid, is a well-known and challenging task in VOF methods. In this study, we evaluate different contact angle boundary conditions, applied to the volume fraction equation. Our current implementation of the Cox-Voinov dynamic contact angle model, based on [4, 5], correctly models the advancing and receding behavior of an inertia dominated flow. In conjunction, partial and Navier slip boundary conditions for the velocity are also assessed. To validate this implementation for surface tension dominated flows, the case of a droplet spreading over a flat plate without initial velocity is simulated and the results compared with experimental [6] and numerical [7] data. The 2D axisymmetric computational domain is shown on the left of Figure 1. The domain is discretized with 96×96 cells. The time evolution of the non-dimensional droplet radius is shown on the r.h.s. of Figure 1. Stat1 and Stat2 curves are obtained with a constant contact angle, with respectively no slip or partial slip velocity boundary conditions. λNdesigns the slip length, taken as half the grid size in the Stat2 case. The Dyn2 and Dyn3 results are obtained with the Cox dynamic contact angle model. We use the same nomenclature as Legendre [7]. One can clearly see that a static contact angle model leads to too fast droplet dynamics. Using the Cox dynamic contact angle improves drastically the results. Unsurprisingly, when a partial slip boundary condition is used, the droplet dynamics is faster and the droplet reaches its equilibrium earlier. However, when the no-slip boundary condition is used together with the Cox model (Dyn2 green dashed curve), the droplet becomes unstable and periodically contracts. This unexpected behaviour is currently under investigation. It disappears when the slip length is reduced from ∆/2 = 0.015625 to a low value like 10−8: See the magenta curve on Figure 1, which is identical to Dyn2 green dashed curve before time τ= 150. In an attempt to further validate this dynamic contact angle model, a liquid film with side contact lines is simulated flowing down a vertical flat plate under the effect of gravity, leading to the formation of rivulets. The relation of the initial film width to the number of rivulets formed and the wavelength are compared with the numerical results from [8], obtained with the JADIM CFD code. As a stepping stone for future simulations of liquid films flowing over complex geometries, a confined thin liquid film flowing over an inclined flat plate is simulated. The resulting velocity profiles and thickness of the film are measured in post-processing and compared with the known theoretical Nusselt solution. Future work includes simulating liquid films flowing down inclined plates with several different base plate geometries, incorporating of mass transfer models, and comparing the numerical results with analytical solutions and experimental data. Acknowledgments The authors thank IFP Energies nouvelles and Technical University of Berlin for providing financial support to this project. 128 The 20th OpenFOAM Workshop (OFW20), 30 June - 4 July 2025, Vienna, Austria 0.95 mm 1 mm 3 mm 3 mm Axis No-slip navierSlip 0.01 0.1 1 10 100 1000 τ = t × µV1/3/σ 0 0.2 0.4 0.6 0.8 1 r* = (r-r0)/(rf-r0) Exp. Lavi and Marmur Legendre Stat1 (λN=0) Legendre Stat2 (λN=∆/2) Legendre Dyn2 (λN=0) Legendre Dyn3 (λN=∆/2) isoAdvector PLIC-RDF Stat1 isoAdvector PLIC-RDF Stat2 isoAdvector PLIC-RDF Dyn2 isoAdvector PLIC-RDF Dyn3 isoAdvector PLIC-RDF Dyn (λN=10−8) Figure 1: Spreading droplet over a flat plate: (Left) Configuration of the 2D axisymmetric case. (Right) Non-dimensional droplet radius measured on the bottom plate. Comparison of Legendre results with isoAdvector PLIC-RDF simulations, with or without Navier-Slip boundary condition. References [1] J. Roenby, H. Bredmose, and H. Jasak, “A computational method for sharp interface advection,” Royal Society Open Science, vol. 3, no. 11, 2016. [2] H. Scheufler and J. Roenby, “Accurate and efficient surface reconstruction from volume fraction data on general meshes,” Journal of Computational Physics, vol. 383, pp. 1–23. [Online]. Available: https://www.sciencedirect.com/science/ article/pii/S0021999119300269 [3] L. Gamet, M. Scala, J. Roenby, H. Scheufler, and J.-L. Pierson, “Validation of volume-of-fluid openfoam isoadvector solvers using single bubble benchmarks,” Computers and Fluids, vol. 213, p. 104722. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045793020302929 [4] R. G. Cox, “Inertial and viscous effects on dynamic contact angles,” Journal of Fluid Mechanics, vol. 357, pp. 249–278. [5] J.-B. Dupont and D. Legendre, “Numerical simulation of static and sliding drop with contact angle hysteresis,” Journal of Computational Physics, vol. 229, no. 7, pp. 2453–2478. [Online]. Available: https://www.sciencedirect.com/science/ article/pii/S0021999109004203 [6] B. Lavi and A. Marmur, “The exponential power law: partial wetting kinetics and dynamic contact angles,” Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 250, no. 1, pp. 409–414. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0927775704004777 [7] D. Legendre and M. Maglio, “Comparison between numerical models for the simulation of moving contact lines,” Computers and Fluids, vol. 113, pp. 2–13, small scale simulation of multiphase flows. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045793014003557 [8] G. Lavalle, J. Sebilleau, and D. Legendre, “Rivulet cascade from falling liquid films with side contact lines,” Phys. Rev. Fluids, vol. 5, p. 124001. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevFluids.5.124001 129 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria FAST METHOD FOR PREDICTING RADIATIVE PROPERTIES OF PARTICIPATING MEDIA IN SOLID ROCKET MOTOR AND IMPLEMWENT IN OPENFOAM XUEFAN HAO1, HU ZHANG*, 2 1School of Aerospace Engineering, Xi'an Jiaotong University, Xi'an, 710049, P. R. China, [email protected] *,2School of Aerospace Engineering, Xi'an Jiaotong University, Xi'an, 710049, P. R. China, [email protected] Keywords: Radiative properties, Particle, Combustion gas, Solid rocket motor, Radiative heat transfer Thermal radiation has a significant impact on the internal thermal environment and plume radiation characteristics of solid rocket motor (SRM) with metallized propellants. The two-phase flow coupling with the radiative transfer process of alumina particles and combustion gas molecules determines the thermal environment of SRMs[1]. In order to achieve precise numerical predictions of the internal thermal environment and plume radiation characteristics of SRMs, it is crucial to establish an accurate and efficient prediction method for the radiative properties of alumina particles and multi-component gas. In this study, an improved method for determining the spectral and effective gray radiative properties of alumina particles and multi-component gas at the microscale was proposed. It combines Mie theory and line-by-line calculation method to establish a comprehensive database of spectral and effective gray radiative properties[2]. The code of implementing Mie theory and line-by-line calculation was verified by comparing with literature results[1, 2]. As show in the Figure 1, the database includes several important parameters such as absorption efficiency, scattering efficiency, scattering phase function of individual particle and absorption coefficient of individual gaseous species. Based on this database, a fast prediction method for the radiative properties of combustion products in SRM at macroscale was developed. As show in the Figure 2, it utilizes the data interpolation method and Planck-mean absorption coefficient global model for fast calculation of the spectral and effective gray absorption coefficient, scattering coefficient and scattering phase function of mixture of particles cloud and multi-component gas. The computational time was reduced by two orders of magnitude compared to applying the Mie theory and line-by-line methods directly. It is especially suitable for simulating the radiative heat transfer within metallized SRM, where the radiative properties of particles and gases need to be calculated frequently. Finally, using this fast prediction method, a OpenFOAM solver coupled the combustion gas flow (continuous phase), particle flow (discrete phase), and radiative heat transfer of participating media (radiation field) was developed. The calculation flowchart of the solver for one time step is summarized in Figure 3. As show in Figure 3, The continuous phase and discrete phase was coupled by source terms, and the radiation field is coupled with continuous phase and discrete phase by radiative source term. In the solution process of radiation field, the radiative properties of continuous phase and discrete phase for constructing radiative transfer equation was calculated the fast prediction method, instead of directly taking constants as in previous solver [3]. By combining the fast and accurate prediction method of radiative properties of combustion product and multiphysics field solver, a higher simulation fidelity level could be provided for the internal thermal environment in solid rocket motor in the future. Acknowledgements The authors thank all those involved in the organisation of this OpenFOAM Workshop and to all the contributors that will enrich this event. References [1] X. F. Hao, et al., “Radiative properties of alumina/aluminum particles and influence on radiative heat transfer in solid rocket motor,” Chinese Journal of Aeronautics, vol 35, pp. 98-116, Feb. 2022. [2] X.F. Hao, H. Zhang*, “A fast method for predicting radiative properties of participating media in solid rocket motor, from microscale to macroscale,” in 7th Micro/Nanoscale Heat and Mass Transfer International Conference, Nottingham, Aug. 5-7, 2024. [3] OpenCFD, OpenFOAM: The Open Source CFD Toolbox. User Guide Version 10, OpenCFD Limited. Reading UK, July. 2022. 130 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria (a) absorption efficiency (b) scattering efficiency (c) scattering phase function (d) Pressure PMAC of species Figure 1: Radiative properties of individual alumina particle and species Figure 2: Calculation flowchart for the fast method for predicting radiative properties Figure 3: Calculation flowchart of the solver for one time step 131 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria ASYNCHRONOUS MODIFIED CONJUGATE GRADIENT METHOD Ranjith Khumar Shanmugasundaram1, Henrik Rusche 1 1WIKKI GmbH, Wernigerode, Germany Keywords: HPC, One-sided communication, Conjugate gradient method, Asynchronous Linear solvers are undoubtedly one of the key building blocks of numerical software. Their performance is a limiting factor in reducing time-to-solution and energy-to-solution of many applications. Current linear solver technology is often limited at scale by the communication cost. The ratio between the cost of reading two floats from memory and adding them is of several orders of magnitude. It is therefore paramount to reduce the need for communication as much as possible, even at the expense of extra computations. This work puts forward a concept for a novel asynchronous linear solver technology which requires neither global sums nor synchronised pairwise communication. The algorithm developed operates asynchronously from start to finish, removing all synchronisation barriers. This approach enhances the efficiency of computer systems and supports the use of energy-efficient features available in nextgeneration hardware. This is a crucial step towards future exascale supercomputers’ cost-effective and environmentally friendly use. For instance, it will have a significant impact on computational fluid dynamics, a key field for advancing energy-efficient technologies and exploring new methods of energy production. The key idea behind the asynchronous Preconditioned Conjugate Gradient algorithm (aPCG) is to use 1 Krylov-space per MPI rank so that each rank continues to work individually. This method is applied for a symmetric matrix structure and adopts the Polak/Ribiere/Polyak [1] search direction. The algorithm asynchronously calculates the global residual by monitoring boundary updates and comparing the local residual to its previous value. The global residual update happens only when significant changes occur towards termination. Figure 1and Figure 2show the execution time line of a simplified classical PCG and aPCG algorithm. PCG uses a traditional MPI-based two-sided communication model. The key observation is the presence of idle periods (highlighted in red), which occur due to the synchronisation imposed by two-sided, non-blocking MPI communication (MPI Isend and MPI Irecv). These idle periods hinder performance as it takes more time to complete data exchanges. In contrast to the PCG, aPCG utilizes one-sided communication (MPI Put and MPI Get). The small black squares in the graph indicate operations where the MPI window is being filled, meaning data is placed in memory locations accessible to other processes. Using this approach, aPCG significantly reduces idle time, as communication happens in the background while computations continue. This results in a more continuous workflow with fewer delays. Additionally, since processes do not have to wait for message exchanges to complete before proceeding with computations, the algorithm achieves better load balancing and improved parallel scalability. In this presentation, the performance of aPCG algorithm tested on the pressure solver for a DLR combustor [2] (variable density, transient, LES) (DLRCJH) with 3 and 24 million cells originating from the exaFoam project will be presented. Figure 1: The schematic execution timeline of a simplified algorithm using two-sided MPI communication (PCG). 132 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 2: The schematic execution timeline of a simplified algorithm using one-sided MPI communication (aPCG). Acknowledgments a. This project has received funding from the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 101118139. The JU receives support from the European Union’s Horizon Europe Programme. Website: https://www.inno4scale.eu/amcg/ References [1] E. Polak and G. Ribiere, “Note sur la convergence de m´ ethodes de directions conjugu´ ees,” Revue franc¸aise d’informatique et de recherche op´ erationnelle. S´ erie rouge, vol. 3, no. 16, pp. 35–43, 1969. [2] S. Lesnik and H. Rusche, “Gc2: Dlr combustor case and description,” 2023. [Online]. Available: https: //develop.openfoam.com/committees/hpc/-/tree/develop/combustion/XiFoam/DLRCJH 133 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria STUDY OF OPTIMIZED SOLUTION FOR TOPOLOGY OPTIMIZATION BASED ON THE BUOYANTBOUSSINESQSIMPLEFOAM WITH THE CONTINUOUS ADJOINT METHOD JAE SUNG YANG1, SANG DON LEE2, JUNE KEE MIN3, 1Rolls-Royce and Pusan National University Technology Centre, Pusan National University, [email protected] 2NEXTfoam Co., LTD, [email protected] 3School of Mechanical Engineering, Pusan National University, [email protected] Keywords: Topology optimization, Natural convection, Continuous adjoint In this study, a topology optimization solver based on the continuous adjoint method is developed from the buoyantBoussinesqSimpleFoam solver, which mimics the buoyant force assuming the linear variation of density with temperature. For the steady-state laminar flow, the primal governing equations (Eqs. (2) – (4)) are solved. Considering the objective function (Eq. (1)), which minimizes the difference between the temperature distribution and the desired temperature in the domain as considered in the previous study [1], the definitions of adjoint equations (Eqs. (5) – (7)) and sensitivity (Eq. (9)) are derived. Also, the volume constraints (Eq. (8)) is considered. Thus, the optimization problems with constraints are described as follows: Minimize 𝐽=0.5 ∫(𝑇−𝑇)𝑑Ω  (1) Subject to ℜ=−∇∙𝐮 ≅0 (2) ℜ𝐮=(𝐮∙∇)𝐮+∇𝑝−∇∙󰇟2𝜈𝐷(𝐮)󰇠−g1−𝛽𝑇−𝑇+𝛼(𝛾)𝐮≅0 (3) ℜ=𝐮∙∇𝑇−∇∙󰇟𝐾(𝛾)∇𝑇󰇠≅0 (4) ℜ=−∇∙𝐯 ≅0 (5) ℜ𝐯=−∇𝐯∙𝐮−(𝐮∙∇)𝐯+∇𝑞−∇∙󰇟2𝜈𝐷(𝐯)󰇠+α(𝛾)𝐯+𝑇∇𝑇 ≅0 (6) ℜ=−𝐮∙∇𝑇−∇󰇟𝐾(𝛾)∇𝑇󰇠+𝛽𝐯∙g+(𝑇−𝑇)≅0 (7) ∫𝛾 𝑑Ω −|Ω| ≅0 (8) The solid isotropic material using the penalization (SIMP) method [2] is considered for permeability (𝛼) and thermal conductivity (𝑘) to distinguish the solid and fluid regions. The sensitivity (Eq. (9)), which is the variation of augmented objective (ℒ) with respect to design variable (𝛾), is derived for the definition of objective and SIMP function, as follows: 𝜕ℒ 𝜕𝛾 ⁄=󰇟𝐮∙𝐯(𝜕𝛼 𝜕𝛾 ⁄)+ ∇𝑇∙∇𝑇(𝜕𝑘 𝜕𝛾 ⁄)󰇠∫𝑑Ω  (9) To obtain a distinct optimized structure, the Helmholtz PDE filtering (Eq. (10)) and Heaviside step projection (Eqs. (11) and (12)) techniques are used, as follows: −𝑅 ∇𝜑+𝜑=𝜑 (10) 𝛾=0.5𝑒𝑥𝑝−𝛿(1−2𝛾,)−(1−2𝛾,)𝑒𝑥𝑝(−𝛿) (𝑓𝑜𝑟 𝛾,0.5) (11) 𝛾=0.5 +0.51−𝑒𝑥𝑝−2𝛿𝛾,−0.5−2(𝛾,−0.5)𝑒𝑥𝑝(−𝛿) (𝑓𝑜𝑟 𝛾,0.5) (12) where 𝑅 and 𝜑 are filtering radius and scalar value of filtering target, respectively. In this study, the sensitivity and design variables are considered as the filtering target for a stable solution. Especially, the smooth variation of steepness 134 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 1: Algorithm and results of developed topology optimization for 𝜽 𝒅 = 0.1, |𝛀| 𝐠𝐢𝐯𝐞𝐧 = 30%, and 𝒌 𝒔 = 0.4 W/mK parameter (𝛿) in projection is considered in the form of Exponential Linear Unit (ELU) function for stable optimal solution. The optimality criteria (OC) algorithm [3] for optimizer is implemented in the proposed solver. Fig. 1 shows the proposed algorithm in OpenFOAM and the topology optimization results. Considering the twodimensional natural convection problem of Ra = 105 with Pr = 0.71 by Barakos et al. [4], the optimized structures minimizing the objective for the desired temperature are obtained, satisfying the volume constraints. Because the developed solver considers parallel computing, an optimized structure for three-dimensional or large-sized problems can be proposed using the developed solver. For dimensionless desired temperature (𝜃) 0.1 with 30% solid volume constraints, having 15 times higher thermal conductivity of solid material than that of the working fluid (air), the optimized structure of two-dimensional problem shows a 22.6% reduced objective value for desired temperature. And, for a three-dimensional problem, the 15% reduced objective value compared with the objective value before optimization is achieved with the satisfaction of solid volume constraints. Results show the physical characteristics to minimize the objective, including the protruding body and guide duct. The protruding body, at the top near the hot wall, interrupts the rising hot plume to prevent the hot region from expanding into the domain. On the other hand, the guide duct spreads the descending cold flow to extend the cold temperature region for the whole domain using the fluid flow and conduction heat transfer, due to the low value of the desired temperature. Acknowledgements This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (No. RS-2024-00449088). References [1] M. Ruberto, “An adjoint based topology optimization for flow including heat transfer,” POLITECNICO DI MILANO 2017. [2] L.H. Olesen, F. Okkels, H. Bruus, “A high-level programming-language implementation of topology optimization applied to steady-state Navier–Stokes flow,” nt. J. Numer. Methods Fluids, vol. 65, pp. 975-1001, 1998. [3] O. Sigmund, “A 99 line topology optimization code written in Matlab,” Struct. Multidiscip. Optim., vol.21, pp. 120-127, 2001. [4] G. Barakos, E. Mitsoulis, D. Assimacopoulos, “Natural convection flow in squarecavity revisited: Laminar and Turbulent models with wall functions,” Int. J. Numer. Methods Fluids, vol. 18, pp. 695-719, 1994. 135 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NUMERICAL MODELLING OF SCOUR AROUND OFFSHORE STRUCTURES Ranjith Khumar Shanmugasundaram1, Sergey Lesnik1, Henrik Rusche 1, V. S. ¨ Ozg¨ ur Kirca 2,3, B. M. Sumer2 1WIKKI GmbH, Wernigerode, Germany 2BM SUMER Consultancy & Research, ITU ARI Teknokent 1, 15, Sariyer, Istanbul, Turkey 3Istanbul Technical University, Hydraulics and Marine Sci. Res. Center, Sariyer, Istanbul, Turkey Keywords: Floating offshore wind farm, Scour, Sediment transport, Finite area method Floating offshore wind farms are increasingly recognised as a promising solution to harness renewable energy in deepwater locations. However, the stability and longevity of these structures are threatened by various geotechnical challenges, particularly liquefaction and scour processes. Liquefaction, caused by the sudden loss of soil strength due to dynamic loading, poses significant risks to the foundation integrity of floating wind turbines. Similarly, scour, the erosion of sediment, soil, or bed material from the vicinity of a structure due to hydrodynamic forces, can compromise the structural stability and overall safety of these systems. Numerical modelling of these phenomena is addressed within the EU project INF4INITYa. This presentation focuses on the numerical modelling of scour. Scour, the erosion of sediment around offshore structures due to hydrodynamic forces, poses a significant threat to the stability and longevity of marine foundations, such as those used for offshore wind turbines, bridges, and oil platforms. This happens since the flow accelerates around these structures, vortices are generated near the seabed, creating pressure gradients that lift and transport sediment. Over time, this erosion can lead to the formation of scour holes, weakening the foundation and increasing the risk of failure. Understanding the mechanisms behind scour is crucial for designing resilient structures that can withstand these erosive forces, especially as offshore installations are subjected to increasingly dynamic and unpredictable environmental conditions. Larsen [1] developed a numerical model based on solutions to Reynolds-averaged Navier-Stokes equations, coupled with additional bed and suspended load descriptions forming the basis for sea bed morphology. The model uses a mass-conserving technique developed by Jacobsen [2], incorporating scalar contributions to the bed evolution equation (Exner equation), introducing a sand sliding routine, and evaluating the accuracy and conservation properties of various interpolation approaches. The bed evolution is modelled using the finite area method from the OpenFOAM framework. The model has been validated for a scour around a vertical cylinder by Baykal et al [3] and also for backfilling processes around a circular pile by Baykal et al [4]. In this contribution, the results of the model will be presented by validating against benchmark scour cases. Also, the preliminary results with the structural configurations developed under the INF4INITY project will be discussed. Additionally, the parallel computational efficiency of the model will be evaluated, and potential improvements in both numerical modelling and performance will be addressed. Acknowledgments a. The abbreviation stands for INtegrated designs for Future Floating oFFshore wINd farm TechnologY. INF4INITY is funded by the European Union under Grant Agreement No. 101136087. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Climate, Infrastructure and Environment Executive Agency (CINEA). Neither the European Union nor the granting authority can be held responsible for them. Website: https://inf4inity.com/ References [1] B. E. Larsen, D. R. Fuhrman, C. Baykal, and B. M. Sumer, “Tsunami-induced scour around monopile foundations,” Coastal Engineering, vol. 129, no. July, pp. 36–49, 2017. [Online]. Available: https://doi.org/10.1016/j.coastaleng.2017. 08.002 [2] N. G. Jacobsen, “Mass conservation in computational morphodynamics: uniform sediment and infinite availability,” International Journal for Numerical Methods in Fluids, 2010. [3] C. Baykal, B. M. Sumer, D. R. Fuhrman, N. G. Jacobsen, and J. Fredsøe, “Numerical investigation of flow and scour around a vertical circular cylinder,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 373, no. 2033, 2015. 136 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria [4] ——, “Numerical simulation of scour and backfilling processes around a circular pile in waves,” Coastal Engineering, vol. 122, no. May 2016, pp. 87–107, 2017. [Online]. Available: http://dx.doi.org/10.1016/j.coastaleng.2017.01.004 137 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria TURBULENT GAS-LIQUID SIMULATIONS OF JET KNIFE FOR HOT-DIP GALVANIZING OF WIRES EMAD TANDIS Department of Mechatronics and Biomedical Engineering, Aston University, [email protected] Keywords: Hot-dip galvanising, Large eddy simulation, Volume of fraction, Jet knife, Turbulence modeling Introduction Hot-dip galvanising is a widely employed method for protecting steel structures from corrosion by applying a durable zinc coating. During the galvanising process, steel wires are immersed in molten zinc at high temperatures and then withdrawn at controlled speeds, dragging a layer of liquid zinc along their surfaces [1]. To regulate the coating thickness, gas jets are directed onto the moving wire, stripping excess zinc and returning it to the bath. Understanding the fluid dynamic interactions between the liquid zinc and the gas jets is critical for optimising coating uniformity and minimising surface defects, as the wire’s coating properties (uniformity and concentricity of the coating thickness) are highly dependent on jet parameters: fluctuating pressure and shear stress associated with the unsteadiness of the turbulence in the jet are underlying causes of waviness on the coating surface [2, 3, 4]. In this work, a Large Eddy Simulation coupled with a Volume of Fluid (LES-VOF) approach is employed to investigate the pressure distributions around the wire at the impingement region, as well as the resulting coating thickness and its non-uniformity at the nozzle exit. Mesh independence analysis was conducted using three grids containing approximately 4.9k, 1.9M, and 6.5M cells. A cylindrical mesh refinement was applied around the jet and near the wire to accurately resolve turbulent structures. Figure 1a shows the refined 6.5M-cell mesh. The incompressible Navier–Stokes equations were solved using spatial filtering for Large Eddy Simulation (LES), with the Wall-Adapting Local Eddy-viscosity (WALE) subgrid-scale model to capture unresolved turbulence stresses. Simulations were carried out using OpenFOAM v2312, employing a second-order Gauss linear scheme for spatial discretisation and a second-order backward scheme for time integration. Figure 1b illustrates the instantaneous velocity field in the jet region. The jet exhibits a flapping motion, which induces pressure oscillations at the impingement zone and can potentially lead to variations in the final coating thickness on the wire. (a) Refined 6.5M-cell mesh around the jet and wire (b) Instantaneous velocity field showing jet flapping Figure 1: (a) Computational mesh and (b) velocity contours at the jet impingement region. Figure 2a shows the effect of wire speed on the coating thickness at the exit. As expected, increasing the wire speed raises the average coating thickness but also leads to greater variation across the surface. On the other hand, Figure 2b presents the effect of flow rate. While an increase in flow rate results in a reduction in the overall coating thickness, the variation in thickness remains relatively unchanged. 144 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 t [sec] 60 80 100 120 140 thickness [µm] ws=0.75m/s ws=1.0m/s ws=1.25m/s ws=1.5m/s (a) Effect of wire speed on coating thickness at the exit. 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 t [sec] 60 70 80 90 100 thickness [µm] Q=0.75m³/h Q=1.0m³/h Q=1.25m³/h (b) Effect of flow rate on coating thickness at the exit. Figure 2: Time evolution of coating thickness at the exit under varying wire speed and flow rate conditions. Figure 3 shows the frequency analysis of the coating thickness and pressure fluctuations. In Fig. 3a, the frequency spectrum of the coating thickness at the exit is presented, while Fig. 3b displays the frequency of pressure at the impingement point. As observed, both thickness and pressure signals exhibit two dominant peaks around 4 Hz and 18 Hz. The presence of matching peaks in both spectra suggests a strong correlation between pressure fluctuations and coating thickness variations. This supports the hypothesis that pressure fluctuations are the main driver of thickness instability at the exit. 0 5 10 15 20 25 30 freq [Hz] 0 100 200 300 400 500 600 A [µm] (a) Frequency of coating thickness variation at the exit. 0 5 10 15 20 25 30 freq [Hz] 0 50 100 150 200 250 300 350 Amp [mbar] (b) Frequency of pressure variation at the impingement point. Figure 3: Frequency analysis of coating thickness and pressure fluctuations. Acknowledgments The authors acknowledges the computing resources provided by the TAURUS High-Performance Computing (HPC) facility at Aston University. References [1] M. Mendez, A. Gosset, B. Scheid, M. Balabane, and J.-M. Buchlin, “Dynamics of the jet wiping process via integral models,” Journal of Fluid Mechanics, vol. 911, p. A47, 2021. [2] G. C. Hocking, W. Sweatman, A. Fitt, and C. Breward, “Deformations during jet-stripping in the galvanizing process,” Journal of Engineering Mathematics, vol. 70, pp. 297–306, 2011. [3] H. So, H. G. Yoon, and M. K. Chung, “Cfd analysis of sag line formation on the zinc-coated steel strip after the gas-jet wiping in the continuous hot-dip galvanizing process,” ISIJ international, vol. 51, no. 1, pp. 115–123, 2011. [4] D. Barreiro-Villaverde, A. Gosset, and M. A. Mendez, “On the dynamics of jet wiping: Numerical simulations and modal analysis,” Physics of Fluids, vol. 33, no. 6, p. 062114, 06 2021. [Online]. Available: https://doi.org/10.1063/5.0051451 145 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria NUMERICAL INVESTIGATION OF PATIENT TECHNIQUE EFFECTS ON PMDI PERFORMANCE FEDERICO ANGIUS1, RICCARDO ROSSI1, ANDREA BENASSI2 1RED Fluid Dynamics, Cagliari, Italy, [email protected] 2Chiesi Farmaceutici, Parma, Italy, [email protected] Keywords: Drug delivery, lungs, pMDI Introduction Pressurized Metered Dose Inhalers, or pMDIs, are devices widely used in the treatment of pulmonary diseases, such as asthma and chronic obstructive pulmonary disease (COPD): using a propellant, a precise dose is administered when the device is activated by the patient itself. The effectiveness of the device could be measured in terms of the dose fraction that can reach the deeper airways. In the present work, the performance of a pMDI will be investigated through the use of numerical CFD simulations, calculating the air flow inside the airways, described by a Eulerian approach, and the drug delivery and deposition into the lungs, described through a Lagrangian approach. Then, a sensitivity analysis of the performance to several parameters related to the patient technique (inhalation profile, device position, inhalation timing) has been performed to evaluate the possible criticality for these devices. A numerical model has been developed through OpenFOAM, preliminary tested and used to perform the required simulations. Numerical model Based on the work of Spasov et al. [1] [2], the OpenFOAM model has been developed to simulate the pMDI coupled within the airways. The numerical mesh has been generated through snappyHexMesh, based on the mesh sensitivity analysis performed in [1]. The model is a Euler-Lagrange one: the air flow, due to the patient’s inhalation, is represented through the Eulerian approach, while the drug is represented through a Lagrangian approach. The Eulerian model represents an incompressible and turbulent flow, where only the mean effect of the turbulence is represented, with an unsteady approach. The employed RANS turbulence model is the kOmegaSST, as suggested by the majority of the authors involved in the field. A fixed time step of 5e-5s has been imposed to have a mean maxCo lower than 1 [1]. Lastly, the model is able to represent different inhalation profiles (steady, time-dependent and with different mean flow rates). The outlflow splitting into the lung lobes has been imposed according to [3], to represent a physiological flow distribution. The Lagrangian model is a oneway coupled and is able to represent the mean pMDI plume feature, such as plume cone angle, plume velocity and plume duration. All the parameters have been set in agreement with the experiments performed on a specific pMDI. The main forces acting on the particles have been taken into account: gravity, spherical drag, Brownian motion (due to the small diameter of the simulated particles) and the stochastic dispersion. The wall-particle interaction has been modeled through the stick algorithm of OpenFOAM: the particles that impact the boundaries are trapped there and considered deposited. The time-step employed to calculate the particles trajectories has been set to be smaller than the relaxation time of the smaller simulated particles. The particles diameters have been selected based on the plume Particles Size Distribution curves, calculated through NGI experiments and specifically modified to take into account the losses on the Induction Port, not considered when the NGI data was elaborated. The PSD has been divided into 15 characteristic diameters and for each one 10k particles have been injected. Parametric analysis A first ideal case has been simulated as a reference, considering an inhalation profile with a mean flow rate of 30 l/min, device perfectly horizontal and pMDI actuation perfectly synchronized within the patient inhalation. The performance of the device has been analyzed in terms of extra-thoracic, intra-thoracic, and peripheral deposition, as a percentage of the emitted dose. A parametric analysis has been carried out that includes the effect of the inhalation profile, actuation timing, and device position. 146 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Conclusions The simulations can predict the drug deposition fraction in each of the regions of interest. The PSD of the deposited fraction can be calculated to analyze the dimensions of the particles deposited in each region. The sensitivity analysis to the main parameters involved will allow one to understand how these parameters affect the performance on the device, showing how good practice in the use of the device can provide or not a better effectiveness of the therapy. For a more robust analysis, several phenomena have yet to be implemented in the model, such as the two-way coupling between fluid and particles and the particles evaporation, to have a more realistic representation of the particles delivery process. References [1] G. H. Spasov, R. Rossi, A. Vanossi, C. Cottini, and A. Benassi, “A critical analysis of the cfd-dem simulation of pharmaceutical aerosols deposition in upper intra-thoracic airways: Considerations on air flow,” Computers in Biology and Medicine, vol. 170, p. 107948, 2024. [2] G. Spasov, R. Rossi, A. Vanossi, C. Cottini, and A. Benassi, “A critical analysis of the cfd-dem simulation of pharmaceutical aerosols deposition in upper intra-thoracic airways: Considerations on aerosol transport and deposition,” Pharmaceutics, vol. 16, p. 1119, 2024. [3] J. Elcner, M. Chovancova, and M. Jicha, “The influence of boundary conditions to the flow through model of upper part of human respiratory system,” in EPJ Web of Conferences, vol. 67. EDP Sciences, 2014, p. 02025. 147 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria A NOVEL OPENFOAM-BASED OPTIMIZATION FRAMEWORK FOR AUTOMATICALLY BALANCING FLOW IN PROFILE EXTRUSION DIES GABRIEL WAGNER1, JOÃO VIDAL1, PIERRE BARBAT2, JEAN-MARC GONNET2, JOÃO M. NÓBREGA1 1Institute for Polymers and Composites – IPC, University of Minho, 4800-058 Guimarães, Portugal; [email protected], [email protected], mnobr[email protected]inho.pt 2Michelin Ladoux Campus RDI, 63118 Cébazat, France, Auvergne; [email protected], [email protected] Keywords: Tire Manufacturing; Extrusion Die Optimization; Bayesian Optimization; OpenFOAM; Computational Fluid Dynamics; Flow Balancing In polymeric profile extrusion, where rheology is complex and various processing conditions and post-extrusion phenomena must be considered, proper die design is essential to achieve the desired extrudate profile. Among the factors that can alter the extrudate profile—such as swelling and cooling-induced shrinkage—flow balance is the most critical for successful production. Unbalanced flow at the die outlet can cause geometric distortions in the polymer melt as it exits the extrusion channel, leading to significant dimensional deviations and defects that compromise product quality. Therefore, maintaining uniform flow distribution across all regions of the die outlet is vital. However, achieving this is often challenging due to the complex geometries required, particularly in automotive and tire manufacturing, where asymmetrical and intricate profiles are common. These challenges require special attention and remain a concern for these industries, especially since current die design methods typically rely on trial-and-error iterations, resulting in substantial material and financial waste. In response to these issues, this work is a collaboration between the Institute for Polymer and Composites at the University of Minho and the R&D team at tire manufacturer Michelin. It involves implementing a novel fully automated computational framework to improve the performance of profile extrusion dies, focusing on maximizing flow balance at the die exit and minimizing pressure drop in the channel. It integrates extrusion die flow simulation in OpenFOAM with an automatic geometry update method in FreeCAD [1] for parameterized CAD design and employs a Bayesian optimization algorithm from Scikit-Optimize [2] to propose new parameter combinations for testing. To automate these steps, a complex framework composed of Python code and Bash scripts was developed to facilitate communication between OpenFOAM simulation results and FreeCAD geometry updates based on analyzed outcomes, executing all pre-, current, and post-processing tasks throughout the optimization loop. The calculation begins with an initial parameterized geometry modeled in FreeCAD, which is subsequently exported to OpenFOAM for meshing with cfMesh [3]. A simulation is then performed for a non-isothermal, non-Newtonian flow using a custom solver designed to model temperature-dependent viscosity with the Bird-Carreau-Arrhenius model, incorporating viscous dissipation and a novel boundary condition to replicate the thermal regulation used in the experimental process. For optimization, the outlet cross-section is divided into several elemental sections to quantify outlet flow distribution and a faster convergence strategy, proposed in Vidal et al. [4], is used. Channel pressure drop (∆𝑃) and velocity uniformity (𝑈𝑢𝑛𝑖𝑓) at the outlet are quantified based on the simulation data, and the objective function (𝐹𝑜𝑏𝑗) is computed with weighting factors (𝑤1 𝑎𝑛𝑑 𝑤2) for each objective, as shown in Equation 1. 𝐹𝑜𝑏𝑗 =𝑤1∆𝑃+𝑤2(1−𝑈𝑢𝑛𝑖𝑓) (1) Finally, the gp_minimize function [5] from Scikit-Optimize is used to perform Bayesian Optimization based on Gaussian Process regression [6], estimating how the objective function behaves for different combinations of parameters and proposing the next combination to test in order to reach the minimum value. To validate the implemented framework, a case study was conducted using a representative geometry of an extrusion die for tire manufacturing, with a tire tread cross-section defined as the outlet shape of the die. To control flow distribution and achieve minimum values, an obstruction with three design parameters (A, B, and C) was studied, as shown in Figure 1a. Weighting factors for pressure drop and outlet uniformity were defined as 0.2 and 0.8, respectively, and the optimization loop was set to run for 19 trials, with convergence demonstrated in Figure 1b. From the geometries proposed by the Bayesian algorithm, the optimal value was 0.1626 (Trial 19), while the worst value was 0.4966 (Trial 5), representing an enhancement of 67.3%. The comparison of flow balance improvement between Trial 5 and Trial 19 is shown in Figure 1c. 148 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria Figure 1: Case study used to test the optimization framework: (a) geometry and obstruction parameterization; (b) objective function evolution; and (c) outlet velocity field at worst (Trial 5) and best (Trial 19) parameter combinations. Acknowledgements The authors acknowledge funding from FEDER funds through the COMPETE 2020 Programme and National Funds from FCT (Portuguese Foundation for Science and Technology) under projects UID-B/05256/2020 and UIDP/05256/2020. They also thank Alexandre Dourlat from Michelin Ladoux for his valuable recommendations. References [1] “FreeCAD: Your own 3D parametric modeler.” Accessed: Mar. 31, 2025. [Online]. Available: https://www.freecad.org/ [2] “Bayesian optimization with skopt — scikit-optimize 0.8.1 documentation.” Accessed: Mar. 31, 2025. [Online]. Available: https://scikit-optimize.github.io/stable/auto_examples/bayesian-optimization.html [3] “cfMesh introduction.” Accessed: Mar. 31, 2025. [Online]. Available: https://cfmesh.com/welcome/ [4] Vidal, J. P. O., Oliveira, M., Srivastava, A., Sacramento, A., Aali, M., Guerrero, J., & Nóbrega, J. M. (2024). Enhancing extrusion die design efficiency through high-performance computing based optimization. Meccanica, 1-12. [5] “skopt.gp_minimize — scikit-optimize 0.8.1 documentation.” Accessed: Mar. 31, 2025. [Online]. Available: https://scikit-optimize.github.io/stable/modules/generated/skopt.gp_minimize.html#skopt.gp_minimize [6] “1.7. Gaussian Processes — scikit-learn 1.6.1 documentation.” Accessed: Mar. 31, 2025. [Online]. Available: https://scikit-learn.org/stable/modules/gaussian_process.html 149 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria PORTING AND VALIDATING AN OPENFOAM SOLVER FOR URBAN MICROCLIMATE SIMULATION IN TROPICAL REGIONS GABRIEL M. MAGALH ˜ AES1, JUAN A. ACERO2, BORIS HULJAK3, J. MIGUEL N ´ OBREGA4, FRANCISCO CHINESTA2,5,6 1PIEP, University of Minho, Azur´ em Campus, 4800-058 Guimar˜ aes, Portugal,[email protected] 2CNRS@CREATE Ltd., 1 CREATE Way, CREATE Tower, 138602, Singapore, [email protected] 3CNRS@CREATE Ltd., 1 CREATE Way, CREATE Tower, 138602, Singapore, [email protected] 4IPC, University of Minho, Azur´ em Campus, 4804-058 Guimar˜ aes, Portugal, mnobre[email protected] 5Arts et Metiers Institute of Technology, 155 Boulevard de l’Hˆ opital, Paris, France, [email protected] 6Keysight Technologies, Inc., 92220 Bagneux, France Keywords: Environmental flow; urban microclimate; thermal multi-region simulation; urban simulation. Urban microclimate refers to the local atmospheric conditions within urban areas that differ from the surrounding regions, often resulting in slight to substantial temperature variations. Accurately assessing this process is crucial for defining strategies to ensure pedestrian thermal comfort. In addition to physical sensors placed strategically, numerical simulations play a vital role in this field. There are open-source codes like PALM-4U [1] and commercial software like ENVI-met dedicated to these simulations. In recent years, OpenFOAM has also been utilized in this area. One framework, urbanMicroclimateFoam (uMF) [2], is an open-source solver for modeling coupled physical processes in urban microclimates, developed by a team from the Department of Mechanical and Process Engineering at ETH Zurich. It is a multi-region solver that includes an air subdomain and subdomains for porous urban building materials. A computational fluid dynamics (CFD) model addresses turbulent, convective airflow and heat and moisture transport in the air subdomain, while a coupled heat and moisture (HAM) transport model manages the absorption, transport, and storage of heat and moisture in the porous building materials. The original version of uMF is implemented for the OpenFOAM Foundation using versions 6 to 8, incorporating its own implementation of the view factor radiation model, including solar load effects. Additionally, the cases where the code is typically applied are generally related to non-tropical areas using relatively simple meshes and domains. This work ports uMF to OpenFOAM v24.06, aiming to leverage recent implementations of the view factor radiation model and solar load calculations, enabling the simulation of complex domains at a reasonable computational cost. In uMF, special boundary conditions transfer information between regions, along with solvers for each region, preand postprocessing utilities, and libraries for material and radiation models. Since the standard view factor model in OF v24.06 includes all necessary radiation implementations, the challenge was to adapt the radiation structures in uMF, which differ significantly between OF Foundation v8 and OF v24.06. After porting all relevant components to the code, validation was performed to ensure the accuracy of the work within a complex simulation framework designed to model various interconnected physical phenomena. Figure 1: Domain used in the simulation: the red buildings represent the AOI, while the grey buildings, in the distant zone, are modeled as impermeable porous regions. The area around the sensors, where accuracy is crucial for evaluating thermal comfort, is called the area of interest (AOI). The mesh in this area is more refined. In addition to the AOI, some distant structures are relevant to the flow features. To incorporate the effects of these structures on the flow without significantly increasing computational costs, these buildings are treated as impermeable porous regions, where the velocity is null, the temperature is constant (28 °C), and the mesh is coarser. The AOI is represented by the red structures in Fig. 1, while the distant buildings modeled as porous blocks are shown in grey. For validation, an area of Singapore, a tropical country with weather conditions very different from Europe (where the original version of uMF was primarily tested), was used. Data from three sensors in the area, including temperature, humidity, and wind 150 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria speed inside the AOI from June 2020, were provided. The simulation utilized the ported uMF, coupled with the Weather Research and Forecasting (WRF) multi-layer urban canopy model (MLUCM). The coupling between the WRF-MLUCM mesoscale model and uMF is achieved through oneway coupling from the WRF simulation output into the initial and boundary conditions of uMF. (a) (b) (c) Figure 2: Ambient temperature over time in sensors (a) P1, (b) P2, and (c) P3. Table 1: Root Mean Square Error (RMSE) between uMF and sensor data Sensor RMSE [◦C] P1 1.338 P2 1.044 P3 1.134 The numerical results of air temperature using uMF are presented in Fig. 2 for the three sensors. Tab. 1 shows the root mean square error (RMSE), a quantitative measure for evaluating result quality. The RMSE is calculated considering N, the number of points over time, by RMSE =v u u t 1 N N X n=1 (Tnumerical −Tmeasured)2.(1) As shown in Fig. 2 and Tab. 1, the temperature evolution over three diurnal cycles aligns well with the sensor measurements. These results from a complex case study enhance confidence in the code porting, the standard view factor radiation model with the solar calculator in OF v24.06, and the adopted methodology (e.g., treating distant structures as impermeable porous regions). Some differences between the measurements and uMF at the beginning and end of the simulations stem from the boundary conditions of the WRF mesoscale model, which do not fully replicate the existing weather conditions (e.g., cloud cover levels). Another source of the observed minor differences is that in the final hours of the last day the sensor data indicate precipitation, which is not yet addressed by the calculation framework. Acknowledgments The authors acknowledge funding from FEDER funds through the COMPETE 2020 Programme and National Funds through FCT (Portuguese Foundation for Science and Technology) under projects UID-B/05256/2020, UID-P/05256/2020. The authors thank Prof Jan Carmeliet and Dr. Aytac Kubilay from the Department of Mechanical and Process Engineering of ETH Zurich, and Dr. Clement Nevers from the Universit´ e de Sherbrooke (Canada) for their guidance and help is the use of uMF. References [1] B. Maronga, M. Gryschka, R. Heinze, F. Hoffmann, F. Kanani-S¨ uhring, M. Keck, K. Ketelsen, M. O. Letzel, M. S¨ uhring, and S. Raasch, “The parallelized large-eddy simulation model (palm) version 4.0 for atmospheric and oceanic flows: model formulation, recent developments, and future perspectives,” Geoscientific Model Development, vol. 8, no. 8, pp. 2515–2551, 2015. [2] urbanMicroclimateFoam, “An open-source solver for coupled physical processes modeling urban microclimate based on openfoam,” https://github.com/OpenFOAM-BuildingPhysics/urbanMicroclimateFoam, 2020, download available at https://github.com/OpenFOAM-BuildingPhysics/urbanMicroclimateFoam (accessed 6 May 2025). 151 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria BUBBLE FORMATION AND FLOW DYNAMICS INDUCED BY LIQUID SLOSHING KOHEI OKITA1, HIROKAZU SUGIZAKI2, TOMONOBU SHIOYASU3, SHU TAKAGI4 1Nihon University, [email protected] 2Komatsu Ltd., [email protected] 3Komatsu Ltd., t [email protected] 4the University of Tokyo, [email protected] Keywords: Liquid Sloshing, Bubble Formation, VOF model, Discrete Bubble model The hydraulic oil tank in construction and mining machines experiences acceleration due to the machine's forward and backward movement, as well as steering and turning. Inside the tank, hydraulic oil sloshes, and bubbles are formed. These bubbles in the hydraulic oil induces cavitation in the hydraulic pump, which results in material damage. Therefore, numerical simulations of bubble formation caused by liquid sloshing are essential for optimizing tank shapes to mitigate bubble formation. Liquid sloshing has been simulated using the volume of fluid (VOF) model to capture the gas-liquid interface. However, representing smaller-scale bubbles requires a large number of grids, even with the use of adaptive mesh refinement. In the present study, a hybrid model solver was developed based on atomizationFoam [1]. The solver utilizes a hybrid approach combining the VOF method and the discrete bubble model (DBM), where bubbles smaller than the grid cells are represented using a Lagrangian formulation with bubble translation motion taken into consideration. The solver was employed to simulate bubble formation induced by liquid sloshing, and its correlation with flow dynamics was analyzed. A closed tank with dimensions 144×144×34mm, filled with liquid to a height of ℎ=85mm, was driven by a sinusoidal oscillation at a frequency of 𝑓=2Hz and an amplitude of 𝐴𝑥=10mm. The liquid phase is water and the gas phase is air. Figure 1 shows time evolution of a gas-liquid interface and discrete bubbles. At 𝑡=840ms, a portion of the wavehead plunges into the interface. Then the gas phase follows and is entrained in the liquid, resulting in the formation of a discrete bubble at 𝑡=845ms. By 𝑡=920ms, VOF and discrete bubbles are formed near right sidewall because of the fragmentation of stretched gas phase. On the other hand, the downward interface flow near the left sidewall entrains the gas phase into the liquid phase, leading to bubble formation at 𝑡=975ms. In the presentation, simulation results will be compared with the corresponding experimental results for validation. Figure 1: Time evolution of a gas-liquid interface represented by an isosurface with 𝜶𝑳=𝟎.𝟓 and discrete bubbles visualized as white particles (diameters doubled for clarity). References [1] M. Heinrich and R. Schwarze, “3D-coupling of Volume-of-Fluid and Lagrangian particle tracking for spray atomization simulation in OpenFOAM”, SoftwareX, vol. 11, 100483, 2020. 152 The 20th OpenFOAM Workshop (OFW20), 30 June – 4 July 2025, Vienna, Austria A NOVEL UTILITY TO CALCULATE THE FORCE AND TORQUE FOR VISCOELASTIC MATERIALS MOHAMMADREZA AALI1, JORGE PEIXINHO2, J. MIGUEL NÓBREGA3 1 Institute of Polymer Processing and Digital Transformation (IPPD), Johannes Kepler University Linz, Altenberger Straße 69, 4040 Linz, Austria, [email protected] 2 Laboratoire PIMM, CNRS, Arts et Métiers Institute of Technology, Cnam, 151 Boulevard de l’Hôpital, Paris 75013, France, [email protected] 3Institute for Polymers and Composites – IPC, University of Minho, 4800-058 Guimarães, Portugal; [email protected] Keywords: Viscoelastic fluids, Force and torque calculation, OpenFOAM utility Viscoelastic fluids, which exhibit both viscous and elastic behavior, play a crucial role in a wide range of engineering and industrial processes, including polymer and food processing, additive manufacturing and biomedical applications. A comprehensive understanding of the mechanical interactions between these fluids and the boundaries of such flow systems is essential for accurately predicting performance, optimizing equipment and processes, and validating constitutive models. Specifically, the evaluation of forces and torques exerted by viscoelastic fluids on surfaces is critical for interpreting flow-induced stresses, assessing mechanical loads, and deriving rheological properties under complex flow conditions. This work presents the development and implementation of an innovative force and torque computation utility specifically designed for viscoelastic materials within the open-source computational fluid dynamics (CFD) framework OpenFOAM. Unlike existing force utilities, which are primarily intended for inelastic fluids and neglect the non-linear stress contributions characteristic of viscoelastic fluids, the proposed tool incorporates the complete extra stress tensor derived from user-specified viscoelastic constitutive models for two-phase flows. This enhancement facilitates the precise calculation of hydrodynamic forces and torques in simulations involving viscoelastic behavior. To verify the accuracy and applicability of the implementation, a benchmark case was conducted simulating a polymer melt undergoing extensional deformation. The setup utilized corresponds to a Sentmanat Extension Rheometer (SER), in which a parallelepiped polymer sample is elongated between two counter-rotating drums, and the torque was computed using the newly developed utility. The torque data obtained were employed to calculate the transient extensional viscosity of the polymer melt. The resulting extensional viscosity evolution, illustrated in Figure 1, demonstrates strong agreement between numeric predictions and semi-analytical results, thereby validating the accuracy of the implementation and confirming the physical consistency of the computed stress contributions. Figure 1: Evolution of the extensional viscosity of a viscoelastic material in an elongational rheometry test Acknowledgements The authors acknowledge funding from Project n° 49209UA - PHC PESSOA 2023 and FEDER funds through the COMPETE 2020 Programme and National Funds from FCT (Portuguese Foundation for Science and Technology) under projects UID-B/05256/2020 and UID-P/05256/2020. The authors also acknowledge the support of computational resources provided by the Leonardo supercomputer at CINECA, which was instrumental in performing the simulations presented in this work. References [1] M. L. Sentmanat, “Miniature universal testing platform: From extensional melt rheology to solid-state deformation behavior,” Rheol. Acta, vol. 43, pp. 657–669, 2004, doi: 10.1007/s00397-004-0405-4. 153