This dissertation presents the development, validation, and application of an integrated Computational Fluid Dynamics (CFD) solver, which is built upon the opensource OpenFOAM framework, to predict and ultimately reduce thermal nitric oxide (NO) formation in industrial rotary kilns without sacrificing process productivity. The work is structured into three parts: theoretical foundations, solver implementation and validation through increasingly complex case studies, and final conclusions and recommendations.

The theoretical part of the dissertation established a comprehensive mathematical framework for simulating the complex interactions within rotary kilns, focusing on turbulent combustion, conjugate heat transfer (CHT), and thermal NO formation. Each part of this multi-physics problem is highlighted with a discussion of one or more available models to resolve it, which is often a trade-off between computational speed and accuracy.... ... Read more

Is Julia suited for chemical process analyses? Yes.
This thesis investigates solving ordinary differential equations (ODEs) in programming language Julia, which it does firstly through a simple example of algae populations, and secondly through complex chemical simulations. A progressively more complex framework describing chemical combustion and diffusion is established and implemented. Problems of great stiffness and containing 200+ variables have been integrated in less than 20 seconds, showing the promise this new, open-source and accessible code language.

The thermochemical model that was developed could be improved upon: it has the problem of unexpectedly low flame temperatures because of heat lost to dissociative effects, the current solution for preventing variables overshooting into the negative, it is is ugly and there is ever more optimisation possible. What the program does is run, and it approximately simulates the kinetics and thermochemistry of combustion within acceptable runtimes, thus proving the technical feasibility but lacking the complete implementation of graph-based chemical analysis in Julia. ... Read more

Emulsions, which are oil-water mixtures, are ubiquitous in the food industry, and precise control of their flow and deformation (rheological properties like the storage modulus G'(ω) and the loss modulus G''(ω)) is critical. However, establishing a connection between the microscopic droplet characteristics of emulsions and their macroscopic viscoelastic behavior remains a major challenge. In this thesis, we propose a framework for estimating rheological model parameters of emulsions to address this issue. Our method combines comprehensive rheological testing of mayonnaise-type emulsions with two types of theoretical models: (1) frequency sweep viscoelastic models, including the classical Palierne model, Maxwell model, Kelvin-Voigt model, and combined models; (2) a energy minimization elasticity model (EEI), which accounts for entropic, electrostatic, and interfacial interactions between droplets. Model parameters were optimized via numerical algorithms to ensure that simulated storage and loss modulus spectra closely match experimental measurements. Results indicate that classical frequency sweep models alone cannot achieve optimal fitting, while a combined approach using Kelvin-Voigt and multi-mode Maxwell models achieves consistency with measured data. Additionally, the two adjustable parameters in the EEI model successfully capture droplet rearrangement and interfacial tension effects, accurately reproducing the low frequency elastic plateau. Notably, each emulsion formulation produces a distinct parameter set, highlighting how variations in droplet size, volume fraction, or surfactant chemistry translate into different rheological profiles. By quantitatively correlating microstructure with macroscopic rheology, this study provides insights into emulsion science and offers a theoretical foundation for designing emulsions with customized texture and stability. ... Read more

Magnetostatics play a crucial role in the detection and localisation of naval vessels. Also, minimising a vessel’s magnetic signature is essential to reduce the risk posed by naval mines, which often rely on magnetic detection. This research aims to improve the calculation of magnetic signatures using the Method of Moments (MoM) by implementing it in Julia, a high-performance programming language. A simplified version of TNO’s current MATLAB-based approach is implemented in Julia to establish a baseline for the accuracy and efficiency. Linear basis functions and automatic differentiation (AD) are incorporated into the methodology to explore potential improvements. These extended methods are compared to the baseline to evaluate their performance.
Results show that Julia can be of great value, since it significantly improves the assembly time of the interaction matrix. Point matching is not a suitable approach when using linear basis functions. The Galerkin method shows promising results, though its computational performance remains a significant drawback. Also, using AD shows potential to simplify the implementation of the MoM by eliminating the need for analytical integral expressions. However, AD disappoints in terms of computational performance. Moreover, the AD implementation relies on a mesh-dependent parameter.
... Read more

Combustion is and remains for the foreseeable future an essential chemical process. The physical and mathematical modelling of such processes can help to optimise the design and operation of combustion furnaces which is critical for fuel efficiency. Given the combined complexity and interaction of the gas flow dynamics and the reaction processes however, such modelling can become time consuming and demanding when it is comes to computational capacity. It is for this reason that alternative modelling and computational techniques are of interest. This paper serves as a thorough introduction to modelling this process using Chemical Reactor Networks (CRN) – which is known to be a relatively efficient model and computational approach.

It provides a step-by-step explanation of the CRN approach, as well as a hands-on implementation for a One-Step Mechanism, ie. a combustion process involving only a single stage oxidisation of the fuel. It also introduces the reader to the industry standard CHEMKIN format and the GRI 3.0 data-base, and investigates the possibility of incorporating the state-of-the-art GRI 3.0 database into Chemical Reactor Networks. Following on from this work it is recommended to validate CRN based results against experimental data and modelling results using different techniques such as Computational Fluid Dynamics to gauge both accuracy and computational efficiency. ... Read more

In this Bachelor project we researched the applicability of finite difference methods to solve the multi dimensional Westervelt equation. The Westervelt is an equation describing the propagation of a non-linear wave. The goal of this project was to solve the full Westervelt equation in higher dimensions. We used Julia for our implementations of the finite difference method. The choice of Julia allowed us to easily write high performing code. We were unsuccessful in modelling the full Westervelt equation since we were unable to implement attenuation into our models. We successfully implemented the non-linear wave equation without attenuation for a number of models, including a three dimensional model. The inability to implement attenuation into our models hampered the flexibility of our models. When the non-linearity created a shock front, the top of the wave would
experience numerical errors which quickly made the whole solution unphysical and ultimately unstable. An added implementation of attenuation would dampen the leading edge of the wave preventing it from becoming too steep. Additionally we made two comparisons with the linear wave equation to verify the correctness of our implementations.
We used a uniform grid to discretize the spatial dimensions. After which we used non-stiff FDM solvers with adaptive time stepping algorithms to solve for the remaining time dimension. We analyzed the use of stiff FDM solvers but these proved too computationally expensive to be viable, especially for our 3D model with 64.000 points.
... Read more

In agricultural studies it is often important to predict the performance of genetically different plants. To make sure predictions are done well, it is necessary to make sure they are not influenced by effects of the field on which they are planted. These field effects or spatial effects are in practice often quite complicated and can be due to a wide variety of reasons. To get a better view of these field effects a good mathematical model is desired. In this paper a model is presented which helps to find these field effects. This model tries to estimate the field effect by comparing data of the same plant on different positions of the field. Data is obtained in a finite amount of positions, which means that the model finds the field effect in a finite amount of positions as well. This field effect is found using a cross-validation technique obtained from Tikhonov regularization. The field effect in a finite amount of positions is extended to a field effect in every position of the field. To do this in a good way a kernel method is used, the advantage of which is that it does not depend on a mesh. This kernel method is here applied with a kernel function that is based on Gaussian distribution. This model is applied to several fields of crops to get a view of the performance of the model on real data. ... Read more

After decades of urban growth, mass transport, including buses, will play a significant role in our daily life. Therefore the requirements for buses and their bus door system are increasing. Ventura Systems is a company that is specialized in bus door systems and wishes to gain knowledge on their bus door system using mathematical modeling. This bachelor thesis provides the base for the mathematical models for modeling a bus door system. Multiple models are presented and one model is analyzed in more depth. ... Read more

The energy transition will significantly impact distribution transformers as they will have to deal with more load on them, which is also more variable due to renewable energy sources. However, currently, these transformers are often not monitored with sensors. Therefore, the Dutch network operator Stedin asked to investigate the possibilities of a digital twin for distribution transformers without many sensors. This thesis presents two ways to do so: the currently most used but more empirical loading guide and a more analytical method where we solve Maxwell's equations using the finite element method (FEM). Furthermore, the transition to more renewable energy sources and sources that draw instant power from the grid causes the current to get distorted. This distortion can be mathematically analyzed using harmonic functions, and we will consider the impact of these harmonics in both the loading guide and the FEM model.

The loading guide is a method written in slightly different ways in the IEEE and IEC standards that takes the load on a transformer and the ambient temperature around the transformer to determine the hot spot temperature inside the transformer. By saying that the hot spot temperature is the warmest point inside the transformer, we can determine the percentage of loss of life of a transformer for a particular loading pattern. However, the loading guide only considers one set of parameter values for distribution transformers, which can vary widely in rating, location and ventilation household, leading to an over-or underestimation of the temperature, which in both cases can lead to extra costs. Additionally, the impact of harmonics is an empirical addition to the loading guide and only considers the effect on temperature rise and not the losses.

Therefore, we consider the transformer in a finite element approach. We solve Maxwell's equation on a transformer cross-section with the FEM, resulting in a 2D model that can calculate the losses for a particular geometry. This model is made using COMSOL Multiphysics. Furthermore, we calculate the core and winding losses under a harmonic load, resulting in considerably higher losses than without harmonics.

As the FEM model is quick and straightforward to run, it can serve as a first step towards developing a digital twin of distribution transformer, giving a way to determine the losses analytically. With future development, the model can provide better insight into the temperature distribution in the transformer. ... Read more

Heat transfer by radiation is something you experience everyday: when walking outside in the sun or when cooking your diner. This thesis provides understandable models of the role that radiation plays in the transfer of heat. By studying literature we developed a mathematical basis for the heat transfer model, where we discuss diffusive, convective and radiative heat transfer, the weighted sum of grey gases and the chosen solving method: the discrete ordinates method. Using CONVERGE CFD and MATLAB software, we present several one- and threedimensional simulations of a furnace in which we consider different conditions. The overall conclusion is that radiation transports heat from places with high temperatures to their cooler surroundings, where also conditions as isolation, wall emissivity and gas mixture play an important role. ... Read more

Aluchemie produces anodes for the Aluminium industry and is the largest stand-alone anode factory in the world. Anode baking process The anode baking process is one of the crucial steps in the production of anodes for the aluminium industry. It improves the strength, conductivity of the anode and reduces reactivity during electrolysis. Efficient baking involves uniform heat distribution on the surface of the Anode. Hotspots appear close to the burner due to high local temperature gradients. Hotspots lead to an increase in unwanted NOX emissions. NOX emissions can be reduced by having an efficient burner design to create a wide temperature distribution and subsequently avoid hotspots. The temperature distribution in the furnace is highly dependent on the flow distribution, combustion and the heat transfer. Accordingly, a model needs to be developed that can model the above-mentioned phenomena and predict the NOX emissions.
State of the art review shows that there are many mathematical models available for the functioning and operation of ABF. But only a few of the models are tailored to estimate the emissions. Furthermore, these simulations are performed on the geometries with the simplest burner design, even though it is known that the burner configurations, significantly affect the emissions. Therefore, in the current study, a more sophisticated burner design is considered. Due to the complex nature of the burner design aspects of mesh generation are studied in detail and recommendations are made to improve the quality of the mesh. ... Read more

The chemical reactor network (CRN) approach is a practical tool for precisely predicting the species concentration in combustion processes with low computational cost. This work examines the capability of the emerging Julia programming language and its ecosystem in solving large CRNs. The packages DifferentialEquations.jl and ModelingToolkit.jl are employed to defining and solving stiff ordinary differential equations, for which the implicit time-integration methods Rodas5 and TRBDF2 with the GMRES linear solver are used. The graph structure of reactor networks is constructed by LightGraphs.jl and SimpleWeightedGraphs.jl. The differential equation solver and the graph data structure are connected via NetworkDynamics.jl. It is concluded that Julia is a competent tool for CRNs containing up to 1000 nodes each with 4 species. Julia is capable of simulating pollutant formation in large reactor networks with reasonable time and memory space. ... Read more

This essay shows a two dimensional implementation of the finite element method for the Westervelt equation. To do this the finite element method is first applied to the linear wave equation, then to non-linear diffusion and finally to the Westervelt equation. Both an element by element and a faster vectorized implementation are given for the finite element method. To verify the numerical solution two analytical solutions are used. The first is a one dimensional wave and the second a circularly symmetric wave.

We found that the two-dimensional implementation was successful in computing the Westervelt equation. The error of the solution scales with the mesh size with a power of around 1.7. It was also found that the time step used to compute the solution needs to be small enough for the implementation to converge. ... Read more

In the present era of environment-friendly and clean combustion systems, there is an increasing demand for fast and accurate tools for emission predictions. The best choice is the CFD-CRN method which is a combination of computational fluid dynamics (CFD) and chemical reactor network (CRN) for decoupled simulation of fluid flow and detailed chemical kinetics. This thesis describes an improved solver implementation for resolving a constructed CRN using only global resolution methods and debunks the notion of needing any form of sequential resolution. This research focuses on further improving the Python based computational tool AGNES, developed at the Delft University of Technology. AGNES can automatically cluster CFD cells into reactors, solve the network and visualise the results [?][?]. This project aims at boosting the performance by reducing computation time and selecting a PETSc [1] based global resolution approach using time integration and Newton’s method. Sandia flame D, a piloted methane-air jet flame (Re=22400) [2], is chosen as the test case. The CRN results for species concentration, mainly NO and CO, are validated with the experimental data and CFD simulation results to ensure no compromise on the accuracy. The solution time and convergence ratewere compared for pre-research AGNES (AGNES v1.1) and the current version (AGNES v1.2). Results show that opting for an entirely global resolution approach is computationally feasible at higher reactor count (more than 500 reactors) and proves superior to the pre-research version. The achieved speedup is around 13% and with smart Jacobian evaluation, this is risen up to 21%. The potential reason for the improved performance is identified as the capability of the global time integration method (in place of local sequential resolution approach) to provide sufficient convergence with an increasing number of reactors and increased complexity. Moreover, the solver can be further augmented by selecting an efficient way for Jacobian calculations either by using automatic differentiation or simplifying the current Python loop approach. ... Read more

This thesis covers the power load flow problem and three numerical methods which can be used as a solution. We firstly discuss the concepts from electrical engineering required to discribe the problem. This includes alternating current and voltage and admittances of various electrical elements. We then use these concepts to derive a set of equations that govern how power flows through a network. To solve these equations, we consider the Newton-Raphson, line-search and trust-region algorithms. For these methods, we also look at its convergence and the amount of computational power required. We conclude in the cases where the parameters of the electrical network are within reasonable bounds, the Newton-Raphson algorithm is sufficient in precision and requires the least computational power. Otherwise, one of the other two methods may be tried. ... Read more

This thesis produces a pre-characterization numerical model capable of handling and calculating electromagnetic fields within a rectangular reverberation chamber near its lowest usablefrequency at 200 MHz. As reverberation chambers strive to have high electric field uniformityto meet field uniformity standards, high electric fields, having the property of being naturallymore uniform, are out of the scope of interest. Additionally high frequencies require more computational memory and CPU time.The model is made using a finite element method based modelling software called Comsol Multiphysics. The modelled reverberation chamber consisting of an antenna, a reflectiveshielded chamber and a Z-fold mode-stirrer is gradually build up. This means that first ananalysis of only the antennas will be made, thereafter the antennas will be put into a reflec-tive shielded chamber environment and finally the antennas are put into a reflective shieldedchamber environment with a z-fold mode stirrer. Furthermore, an extra situation will be considered in which a dielectric object will be added to a shielded chamber environment excited by an antenna without mode stirrer. The effect of an added dielectric object will be studied because of physical interest and completeness and not as added intermediate step to build a reverberation chamber environment. As all situations have similar difficulties in modelling and measuring, the gradual development of the reverberation chamber will allow for the best error analysis of the model. The model replicates the setup of a real reverberation chamber located at Comtest, Zoeterwoude. The reliability and accuracy of the model is studied by comparing the modelled electric field to the measured one. It was found that all models showed a good resemblance between simulated and measured electric fields above 60 MHz except the most complex reverberation chamber model. The simulated field uniformity expressed as standard deviationis twice as high as measurements suggest. The model can therefore only be used as worst-case scenario prediction for the field uniformity. Two kinds of antennas are used during the modelling and measuring phase, a 3104c biconical antenna used in the frequency range 25-200 MHz and the 3146a log-periodic dipole array antenna in the frequency range 200-1000 MHz. The electromagnetic radiation pattern in the far field domain was modelled and corresponded to the expected omni-directional and directional field pattern respectively. Next the antennas were placed in a highly reflective shielded chamber with dimensions (4.05 m×2.55 m×2.925 m). By taking data in specific slices from the model and comparing these to the measurements at the same points it was concluded that the model does not ac-curately predict the electric fields below 60 MHz. Above 60 MHz the model does predict the general electric field pattern in the chamber. It however does not predict local maxima or minima of the electric field.The same comparison was done for the situation in which a dielectric object was added to the setup. The dielectric object was chosen to be a container filled with water (0.27 m×0.565 m×0.269 m). This container is placed in the formerly empty shielded chamber to change the inner electric field. Water was used due to its favourable properties for reflectivity. The same behaviour of the model was observed with this added dielectric object. Finally the electric field in the reverberation chamber at Comtest (5.03 m×3.97 m×2.85 m) was modelled and solved using a GMRES algorithm with geometric multigrid preconditioning.The preconditioning in the GMRES allows for faster convergence. The field uniformity wascomputed for both the model and the measurements as outlined by IEC 61000-4-21. This wasdone for 4 and 12 stirrer rotation positions. Both showed that the model had a less uniform field compared to the Comtest reverberation chamber. The Comtest reverberation chamber complied with the electromagnetic compatibility requirements for measurements using 12 stirrer positions whereas the model did not. However the 4 stirrer position model, which made use of perfect electric conductor boundary conditions, showed a maximum of 13% increase in field uniformity when steel walls were used instead. To show a glimpse of an innovation that put the Delft University of Technology on the map, the shifted laplacian preconditioner is briefly discussed. As an intermediate step in solving a problem with little to no damping, a complex preconditioning matrix is used. It is shown that for an increasing imaginary shift in the Helmholtz problem, expressed as an increasing electrical conductivity σ, the number of iterations needed to reach a relative tolerance smaller than 0.01 decreases. At last, the effect of a shifted laplacian contribution to a multigrid preconditioning on the convergence speed is studied for an non-damped pressure acoustic Helmholtz problem. It is shown that the added contribution slows convergence in a simple geometry, while a more complex geometry cannot be solved without this contribution. ... Read more

A high-resolution hydrological solver that solves on time-adaptive grids is created. The solver uses h-refinement strategies which adapts the resolution of the grid using a posteriori measure. It is based on an open-source scientific computing library named PETSc. With the use of highly scalable time-stepping solvers provided by PETSc, an optimised numerical model is obtained. Validation tests are conducted using dam-break scenarios and flow over humps. Strong scaling tests are conducted to test the scalability of the solver. The model is found to be robust and could thus provide accurate flood-forecasting information. ... Read more

In this work, the parareal algorithm is analysed and executed on the model for combustion of methane. The parareal algorithm is designed to generate an approximation to an initial value problem faster than a serial numerical time-integration method by using two propagators, the coarse propagator and the fine propagator. With the use of two different propagators, some computations can be carried out in parallel, which leads to a faster method. In this research, the parareal algorithm is executed on the two-step mechanism for combustion of methane. The temperature rise due to the combustion is assumed to be zero. The model is implemented in Python and with the use of the librarymultiprocessing, computations are executed in parallel. Different time-integration methods are implemented that can be used in the coarse and fine propagator. In this research, we will focus on the case that the both propagators use the same time-integration method. One can distinguish the propagators by using a different time-step for a chosen time-integration method. With the use of the absolute error the accuracy can be examined. Because the analytic solution to the problem for combustion of methane is unknown, a time-integration method, from which we know that it gives a small absolute error, is used as representation of the analytic solution. The parareal algorithm executed on the model for combustion of methane gives an accurate result for the right choice of propagators. However, for this choice of propagators, the parareal algorithm does not result in a significant speedup compared to the fine propagator in serial, assuming that we have enough processors available. This is because the running time of the propagators do not differ much. To generate a better speedup, two different time-integration methods can be considered for the propagators. Moreover, the model for combustion of methane can be divided into more than two partial reaction and then themulti-level parallelization [1] can be examined. ... Read more

Decreasing NOx emissions is becoming increasingly important as it has many lifechanging implications on humanity and nature. In this work several models were constructed to simulate the combustion of methane in a cement rotary kiln, and it was investigated whether oxygen-enhanced combustion leads to less NO_x production for the obtained models. This was done using the Cantera software package with Python. Starting off with a single zero-dimensional reactor equipped with a reduced one-step global reaction mechanism, which was expanded to include two-step and four-step reaction mechanisms. Combustion inside this homogeneous reactor was simulated for various stoichiometric conditions for an initial temperature of 1000K. Subsequently, a one-dimensional model of chained reactors to simulate flow was considered, equipped with both the two-step and four-step mechanisms. This was realized using the scalar convection-diffusion equation to compute the flow throughout the reactor. All aforementioned models were adjusted to replace air with oxygen-enhanced air, containing higher levels of oxygen for every iteration. The simulated
temperature evolution was examined using the exponential relationship of temperature and thermal NO_x. ... Read more

Anode baking is a huge part of the Aluminum industry. The quality and characteristics of the anode blocks will be formed in the anode baking furnaces. The process itself is very expensive, time consuming and ecologically unfriendly. From the other hand, anode baking process can use its own waste as a secondary fuel source. Here, by waste we mean volatile gases from the baking process. In the factories, raw anode and baked anode are weighted before and after baking process. The difference in the weight counted as a released gas. In this work we will try to build a model and solve it that can find volume of the volatiles. Different parameters of the anode material will be included into the model, and analyzed to find insights how it can be used for future investigations. The results of this work will be compared with experimental results. Finally, recommendations for further works will be given. ... Read more