Circular Image

S. Hickel

info

Please Note

15 records found

The aerospace industry is under increasing pressure to reduce emissions and transition toward a more sustainable future. With recent technological advances and stronger emphasis on environmental protection, liquid hydrogen has gained attention as an alternative fuel for civil aviation due to its a vorable thermo chemical properties, including low minimum ignition energy, wide flammability limits, high energy content, and potential as a zero-carbon alternative to fossil fuels. However, operating hydrogen under such conditions presents challenges, including potentially higher engine-out NOx emissions due to higher flame temperatures, onboard storage difficulties, higher laminar flame speeds, and an adia batic stoichiometric flame temperature higher than that of natural gas. To mitigate NOx formation, water injection has emerged as a possible method due to water’shigh specific heat capacity and latent heat of evaporation absorb heat and lower combustion temperatures, with liquid water injection generally prov ing more effective than steam. Additionally, introducing flame strain has emerged as another potential strategy as highly strained flames have shown promising results for reducing NOx emissions. This thesis aims to qualitatively and quantitatively analyze the effect of water injection on highly strained premixed laminar hydrogen flames through the use of computational methods and in this way infer the possibility of combining these two methods to reduce NOx emissions in hydrogen combustion. In order to do that, Direct Numerical Simulations (DNS) of highly strained hydrogen flames with water injection with different setups were performed to conducted a parametric analysis of the effect of spray injection velocity, droplet diameter, and strain rate on the flame structure and emissions of NOx related species. The results revealed that for the baseline case with water injection a sharp reduction in the presence of key flame radicals and reductions in hydrogen reactivity and domain temperature, which lead to reductions of NNH, N2O, and NO emissions. These effects are enhanced with increasing water injection velocity likely due to the increasing momentum resulting in the evaporation of the droplet occurring closer to the flame front. Regarding the effect of increasing droplet diameter, it was verified that increasing diameters are also associated with larger reductions in NO emissions and radical presence, likely due to the higher droplet volume requiring more energy to evaporate and therefore evaporation occurring closer to the flame front. In a computational setup with a flame with higher strain rate than the baseline case, when injected with water at similar water loading, the reductions in radical compositions, hydrogen rate of production and NO emissions are more significant in the case with higher bulk flame strain rate. From these results, it can be concluded that increasing droplet diameter and water injection velocity has positive effects on reducing emissions of NOx related species. Furthermore, it can also be concluded that higher strain rates enhance the effects of water injection, presenting sharper reductions in key radicals and emissions of NOx related species. ...

Reconstructing Geometric Properties Using Modal Participations

Master thesis (2026) - R. de Vroomen, T.J.C. van Terwisga, L.P. Lagendijk, H.C.J. Wijngaarden, H.C. Neatby, A.H. van Zuijlen, S. Hickel
Propellers are often the sole form of ship propulsion, making their design uniquely important in the vessel's operation. Despite many advantages, they suffer from issues like cavitation, vibrations and underwater radiated noise. Flexible composite materials have been proposed as a way to address these issues and increase the overall efficiency envelope, achieved by a passive pitch distribution reduction in the vessel's wake. Due to an increase in variables, optimization is increasingly dependent on the accuracy and validity of numeric simulations.

The objective of this report is to study experimental deformations of flexible propellersin terms of their underlying parametric properties. These are compared to simulations using the unsteady Reynolds averaged Navier-Stokes equations for fluid calculations, coupled with a finite element model for the structural calculations.

A custom preprocessor was developed and used to reconstruct the full parametric blade geometry from digital image correlation measurements. The first 10 mode shapes of the undeformed mesh were found using modal decomposition, which was fed as input to a flexible point cloud registration. The participation factors for each mode shape were optimized to minimize an error function, the symmetric chamfer distance. The deformed mesh was processed with a software library, PropArt, to recover the underlying geometry parameters.

It was found that the hydrodynamic pitch β was the best predictor of rake, skew, and camber deformations, with increased inertial forces increasing the magnitude of the deformation. Pitch deformations are largely predicted by the blades' inertial loading, and is mostly independent of β. The blade thickness and chord length deformations were negligible. Simulations tended to over predict the deformations of the propeller, especially at low β, being exacerbated by increased inertial loading. At high β, the deformations were small and well predicted.

The relation between pitch deformation and inertial loading suggests an optimal material stiffness per operating condition exists, to be tuned for an appropriate amount of pitch reduction. The lack of correlation to β also poses a challenge in the design as the wake peak is where the change in pitch is desired. The increase in camber also leads to an increase in thrust at low β, which is contrary to the desired wake peak thrust reduction. These relations should be investigated further, for a broader range of material stiffnesses to see if this holds. Finally, the difference between experimental results and simulations suggests a systematic error due to the error scaling with the inertial forces, which should be further investigated.
...
Master thesis (2026) - J.P. Kok, I. Langella, A. Gangoli Rao, S. Hickel
To reduce humanity’s dependency on fossil fuels, hydrogen is increasingly investigated as a promising energy carrier. When converting the chemical energy of hydrogen to the desired form, hydrogen reacts with oxygen to form water. This makes hydrogen very promising for achieving net-zero carbon emissions. For the aviation sector, hydrogen combustion has been suggested as a possible solution.
As with any technology, challenges remain. Whilst combustion modelling using Computational Fluid Dynamics (CFD) is not a solved problem to begin with, hydrogen combustion is made more complex by differential diffusion. This occurs when different species do not diffuse at the same rate. Hydrogen tends to diffuse faster than other species, causing local changes in mixture fraction, in addition to super-adiabatic temperatures and super-equilibrium product concentrations. This effect is stronger in strained or curved flames. Flame strain and curvature also influences the local reaction rate. These phenomena occur on a very small length and time scale. To simulate this directly, very fine meshes and very precise transport equations for each species are required. This can be prohibitively expensive for large scale designs or for design iteration.
In cases where simulating this behaviour directly is not possible, an alternative model is required. In this paper, five variants of such a model are assessed. These models are based on the assumption that the flame can be approximated as an ensemble of 1-dimensional flamelets. These flamelets can be analysed before the CFD simulations are performed. The results are then parametrized by a number of control variables, using a presumed filtered probability density function which aims to include sub-grid scale effects. By changing the flamelet conditions and control variable definition, different manifolds can be generated. In this paper, a number of these Flamelet-Generated Manifolds (FGMs) are compared to a higher fidelity model (using an Eulerian Stochastic Fields (ESF) approach). These FGMs use the mixture fraction, the progress variable and their sub-grid variances as control variables. Five different FGMs are tested, characterized by the progress variable definition (water vs hydrogen mass fraction) and flamelet strain rate (0 (unstrained) vs 3000 vs 6000 vs 13000 s^-1.
The performance of these models is assessed by means of a Large Eddy Simulation (LES) of a combustor with a bluff-body, using the open source CFD software OpenFOAM. The bluff-body causes a recirculation zone in its wake. Additionally, it causes a strain on the flame front, where differential diffusion is expected to cause super-adiabatic temperatures and super-equilibrium mixture fractions.
The results show that the FGMs give a good prediction of the conditionally averaged reaction rates. They are also capable of qualitatively predicting the increase in mixture fraction, super-adiabatic temperature and super-equilibrium product mass fractions. They can therefore be used in a strained combustor setting with lean, premixed hydrogen. However, there are challenges regarding the prediction of temperature, especially for FGMs using a hydrogen based progress variable. This should be investigated in the future. Furthermore, the FGMs predict the reaction rate well despite over-predicting the local effects of strain and differential diffusion. ...
The present Master's thesis investigated the low-frequency unsteadiness characteristic of highly separated transitional oblique shock wave-boundary layer interactions (hereinafter "OSBLIs"). Such phenomena, encountered notably in engine components operating at transitional Reynolds numbers, are relevant due to their impact on aerodynamic efficiency, structural integrity, and system reliability. The primary research question examined how variations in different parameters such as Mach number, Reynolds number and inviscid pressure jump influence a low-frequency shock oscillation mechanism which was previously identified in literature. To this end, experimental studies were conducted in the TST-27 transonic-supersonic wind tunnel at TU Delft, using high-speed and spark-light Schlieren visualizations capture the relevant flow phenomena. These recordings were processed using digital and spectral analysis.

The findings from this study revealed that the investigated transitional OSBLIs exhibited low-frequency shock oscillations strongly correlated with the periodic formation and disappearance of a Mach stem, which was denoted as the "dual domain" phenomenon. Through carefully chosen variations in Mach number and Reynolds number, it was shown that slight adjustments significantly impacted both the presence of the dual domain and the characteristics of the shock oscillations. Moreover, the Reynolds number regime identified as transitional for the natural flat plate boundary layer in previous research was validated.

Another aspect of this thesis involved the implementation of passive flow control techniques, specifically the introduction of thin two-dimensional steps (in height increments of 60 microns) designed to artificially trip the boundary layer. The experimental results demonstrated that even minimal boundary layer tripping significantly dampened the shock oscillations and modified the interaction dynamics of cases where the oscillation mechanism had previously clearly been identified. The frequency analysis confirmed this, as none of the oscillation peaks which had previously been identified were observed when tripping the boundary layer.

A non-dimensional analysis of the dominant oscillation frequencies indicated a consistent Strouhal number convergence around St = 0.33, particularly in cases where the "dual-domain" behavior and high oscillation amplitudes were observed. An increase in the Reynolds number consistently resulted in reduced laminar separation amplitudes and increased oscillation frequencies, which aligned with the theoretical expectations of accelerated boundary-layer transition dynamics.

In conclusion, this study was successful in identifying the main parameters that cause unsteadiness in transitional OSBLIs. It confirmed the existing transitional Reynolds number ranges which had previously been analyzed in the context of weak OSBLIs and the natural boundary layer of the flat plate which was used, and showed that even simple flow control methods with low 2D step heights can effectively reduce shock oscillations. Additionally, a meaningful non-dimensional scaling was done, which can aid in further research and comparison of the phenomena which were investigated in the present thesis. These findings provide a solid basis for future studies aiming to apply this knowledge to more general cases and improve the prediction and design of aerospace components affected by transitional shock-induced boundary-layer interactions. ...
The recent growth of the size of wind farms highlights the need for a deeper understanding of the mesoscale phenomena in a stably stratified atmosphere, such as atmospheric gravity waves, as the effect on the power generation can be significant. This thesis is a numerical study of the effect of wind farm layout on atmospheric gravity wave excitation and the resulting feedback on wind farm performance. Wind farms in this study have varying power density, streamwise or spanwise turbine spacing, aspect ratio, hub height, rotor diameter, shape, or orientation with respect to the freestream, and are situated in a conventionally neutral boundary layer in offshore conditions. Moreover, the wind farms can be horizontally or vertically staggered. To investigate the atmospheric gravity wave excitation, the effect of wind farm layout on the Froude number and inversion Froude number governing the internal and interfacial waves respectively is studied using several high-fidelity Large Eddy Simulations. Then, AGW wavelength and wind farm efficiency are parametrically studied using a fast reduced-order model. Specifically, the non-local efficiency is considered, which is a measure of the global blockage effect induced by the atmospheric gravity waves. It is found that the length scale used in the Froude number must be adjusted to the wind farm, and is dependent on the turbine spacing and farm shape. The inversion Froude number is based on the phase speed of the interfacial waves. It is suggested that the phase speed must be based on shallow-water theory and deep-water theory for small and large turbine spacings respectively. In other words, for small turbine spacings the wind farm acts as an entity, while for large turbine spacings, the farm acts as a collection of individual turbines. Finally, the streamwise and spanwise turbine spacing (and consequently the power density), and the aspect ratio primarily govern the non-local efficiency. ...
An impinging Shock Wave-Turbulent Boundary Layer Interaction (SWTBLI) at Mach 2 was investigated experimentally while implementing two-dimensional Shock Control Bumps (SCBs). The aim was to investigate the changes in the unsteady dynamics while changing the bump ramp angle, tail angle, spanwise shape, and the impinging shock location. Schlieren and oil flow visualisations were used to identify these changes. Unsteadiness was quantified through a spectral analysis based on Welch’s method, and a Canny-based edge-detection algorithm was developed to track the impingement and separation location in the Schlieren images.

An uncontrolled SWTBLI could successfully be generated, and the implementation of the baseline bump showed a replacement of the unsteady separation shock by a steady compression ramp shock originating from the leading edge. Spectral analysis confirmed this behaviour since the characteristic low-frequencies of the uncontrolled separation shock seemed to be removed for this compression ramp shock. An increase in the bump ramp angle showed the progressive generation of a separation shock upstream of the bump with spectral content trending towards the low frequencies of the uncontrolled interaction. On the contrary, no alterations were observed for an increase in tail angle. An upstream impingement produced similar behaviour as the increase in ramp angle; a separation shock was generated upstream of the bump with increased low-frequency spectral content. The downstream impingement, however, did not show any alteration in the same region. Major factors influencing the unsteady dynamics are identified as the impingement location and the ramp shape of the bump.

The developed edge-detection algorithm proved unsuccessful in quantifying the unsteadiness although it could detect the impinging and reflection shock of the interactions. Spatial standard deviation distributions of the interaction revealed increased deviation values in the impinging shock suggesting that the impinging shock fluctuates. Rather, it is suspected to be a form of noise intrinsic to the experimental technique. Therefore, this affects the detection of the edges and the calculated impingement and separation location. Future research is suggested to improve the noise mitigation method in the algorithm. Additionally, the benefits associated with the 2D-SCB would be most noticeable in an integrated approach with an industrial application. Finally, numerical simulations and/or different experimental techniques are suggested for future research to obtain a better quantification of the interaction and unsteady dynamics.
...

Application to the TU Delft LEI V3 kite as a case study

Master thesis (2022) - S. Sen, M. Kotsonis, S. Hickel, J. Casacuberta Puig
The flow of air over a swept wing initially starts in a smooth laminar state (referred to as base flow) and entrains disturbances which subsequently grow and transition the flow to a chaotic, turbulent state. Active efforts have been made to study and control these disturbances, which manifest as stationary crossflow modes. Excrescences in the form of a forward-facing step (FFS) impose a base flow deformation in the form of a rapid near-wall pressure change, flow separation, and strong upwash. This modifies the behaviour of stationary crossflow instability and subsequently leads to an upstream or downstream shift of transition location depending on FFS height. The mechanisms responsible for the above behavioural modification are unknown, motivating the current thesis. The evolution of primary stationary crossflow instability, in its linear growth phase, is studied through a spanwise invariant, synthetic, and idealized rapid base flow deformation imposed by changing the near-wall pressure distribution of a clean swept flat plate (possessing a favourable pressure gradient) via Gaussian-like pressure variations. An energy balance framework is developed that identifies the production term's behaviour as the differentiator between regions of perturbation growth and decay. The behaviour of the production term is described by two competing mechanisms, the first controlled by wall-tangential base flow shear and the second controlled by wall-tangential base flow acceleration or deceleration. The balance of these mechanisms shows that perturbations grow faster than the clean case in regions of wall-tangential base flow deceleration and slower than the clean case in regions of wall-tangential base flow acceleration. Perturbations are even found to attenuate in some cases when a region of wall-tangential base flow acceleration follows a region of wall-tangential base flow deceleration. The modification of energy transfer mechanisms brings into question whether the initial modal stationary crossflow mode deviates from modal character on interacting with the base flow deformation. The Orr mechanism is shown to identify differences in perturbation behaviour from local modal character. However, criteria from the literature hint towards an absence of any non-modal effects. Finally, the extent to which a synthetic idealized base flow deformation mimics the effects of an FFS on deforming the base flow and changing trends of stationary crossflow instability evolution is tested to show the applicability of methods developed in the thesis to instances of natural base flow deformation. The progress in understanding mechanisms by which a deformed base flow affects the linear phase of primary stationary crossflow instability growth leads to suggestions on devices that can be tested to delay this phase of instability growth. These devices could potentially also delay subsequent stages of instability growth and hopefully lead to the development of novel transition delay techniques. ...
 Distributed surface roughness elements characterise Thermal Protection Systems (TPS) typical of supersonic and hypersonic flows. The presence of these distributed roughness elements cause an increase in drag and heat transfer.  As opposed to incompressible flow over roughness elements, there are very few experimental and numerical studies on supersonic flow over roughness. The most fundamental computational technique wherein, all scales of turbulence is resolved is Direct Numerical Simulation (DNS). The cost of performing DNS of fully resolved roughness is even higher than canonical DNS because of the refined mesh needed to solve the roughness elements.  To overcome this limitation, the current thesis aims at exploring low cost alternatives to DNS of fully resolved roughness for studying the effect of drag and heat transfer in supersonic flow over rough walls. DNS results from fully resolved full channel cube roughness for Mach 2 and Mach 4 at friction Reynolds number Reτ = 500,1000  are analyzed. The results from the low cost models are compared against the fully resolved roughness simulated using full channels. Three lower-cost alternative, namely DNS of minimal channel flow of fully resolved roughness, DNS of modelled roughness and resolved RANS are considered. As for the DNS of minimal channel flow, it is found that the velocity shift ΔU+ is predicted accurately and therefore the added drag. However, it cannot be used to predict the temperature field because of lack of outer layer similarity for the thermodynamic statistics. As for the modeled roughness, an extension of the model by Busse and Sandham originally developed for incompressible flows is considered. In this case the roughness geometry is substituted by the additional drag and heat transfer that it induces on the flow, which take the form of source terms in the momentum and energy equations. We perform 17 DNS simulations with modeled roughness and compare the results to the fully resolved simulation. We find that the parametric forcing method is able to predict the velocity shift with good accuracy, although recovering the equivalent roughness height from the model parameters can only be done a posteriori. The model is able to qualitatively reproduce the temperature field, but thermodynamic statistics are inaccurate when compared to DNS of the fully resolved geometry. The final computational technique is RANS. In real case applications, RANS require the use of wall functions, and in the case of rough walls knowledge of the equivalent roughness height ks+ is necessary. We attempt to see if RANS of fully resolved roughness can be used to estimate the velocity shift ΔU+ and therefore  ks+ by limiting ourself to the  linear Spalart-Allmaras (SA) model. It is found to be inaccurate in computing the mean velocity profile at  ks+≈ 40 with improvements in accuracy observed for ks+≈ 80 when compared with the results from DNS for cube roughness element. However, the accuracy is still low to be used for estimating ks+...
Large Eddy Simulations (LES) of a novel type of wing/body junction called the anti-fairing are performed in the current thesis to study the complex turbulent flow physics involved in the junction area and also to obtain a clear understanding of the drag reduction capabilities of the anti-fairing. In regards to that, two separate LES are performed: one for the baseline case with a Rood wing/flat plate combination and another with the Rood wing/anti-fairing combination. A detailed comparative study is performed between the two cases to observe important differences in junction flow characteristics. Both the simulations are performed on a 25 million immersed boundary Cartesian mesh by solving the incompressible Navier-Stokes equations using the in-house finite volume LES solver called INCA. Results from the LES study confirms the existence of the propulsive pressure mechanism of drag reduction for the anti-fairing case, previously proposed by Belligoli et al. However, the results also show that there exists a secondary drag reduction mechanism caused by a combination of increase in approach boundary layer momentum thickness and dampening of the turbulence associated with the horseshoe vortex (HSV) upstream of the wing. This secondary mechanism has been found to be caused by the convex dent present at the start of the anti-fairing geometry. The total drag reduction for the anti-fairing case comes out to be 1.8%. A new parameter called junction drag is defined which accounts for the drag only due to the presence of a junction. The reduction in junction drag obtained for the anti-fairing case is about 6.8%. Apart from the LES analysis, a RANS analysis has also been performed to further investigate the drag reduction capabilities of anti-fairing for different approach boundary layer thicknesses and anti-fairing depths. All the RANS analysis have been performed on a 5 million body-fitted mesh by solving the incompressible Navier-Stokes using the open source finite-volume solver OpenFOAM. Results from the RANS analysis indicate that there exists an optimum depth for the anti-fairing which corresponds to the least drag. Furthermore, it is found that the effect of approach boundary layer thickness is mostly on changing the base drag of the case where no anti-fairing is present, rather than actually affecting the performance of the anti-fairing at different depths. ...

Solution of the vector Laplace and the Stokes’ equation

Being able to solve numerically partial differential equations is fundamental for engineers to evaluate, optimize and improve industrial equipments. The framework of mimetic finite element methods allows engineers to find solutions characterized by strong conservation properties: this may result in a pointwise divergence free-flow field. However, sometimes, the computation of the solution of partial differential equations is time consuming, as a results, to reduce the computational time, engineers and mathematicians have developed hybrid methods.
The objective of this thesis is the development of a hybrid mixed finite element formulation of the vector Laplace equation without spurious modes. Discontinuous elements permit an higher degree of parallelism, and, at the end, a lower computational time. Lagrange multipliers are used to impose continuity between discontinuous elements. These turn out to be not only mathematical features but they are connected to the physical variables of the problem. Furthermore, it has been found that the usage of a new Lagrange multiplier, on the intersection of 4 or more elements, removes the spurious modes. Therefore, the associated system of equation is non-singular. The usage of the hybrid finite element methods reduces the computational time while maintaining the pointwise divergence constraint and the optimal convergence rate of all variables.
At the end, the mixed hybrid formulation is modified to solve the Stokes equations. Lagrange multipliers are used as boundary conditions. Solution of the lid-driven Stokes flow is shown.
...
Master thesis (2018) - Jeroen Kunnen, Marc Gerritsma, Stefan Hickel, Matthias Möller
In 2011 a relatively new type of numerical scheme has been introduced: Active Flux schemes. In this type of scheme an extra degree of freedom is added to the cell interfaces of a regular finite volume grid. This enables the use of non-conservative update methods for these additional variables, as conservation is automatically adhered to by the cell integral values. This increases the order of accuracy of the scheme, while allowing a broader range of update methods. This thesis proposes a new update method based on minimizing the truncation error of a Taylor series expansion. This way, a linear update scheme can be created for each unique stencil. An adaptive mesh refinement algorithm is implemented to conform the mesh to local high-frequency phenomena such as shock waves. A high-resolution simulation shows that the adaptive method reaches error levels of a uniform mesh while using ~9.6 times less computational cells. ...
Master thesis (2018) - Jacob Butler, Richard Dwight, Matteo Pini, Stefan Hickel
The project concerns uncertainty reduction of parameters of a thermodynamic equation of state for a dense gas, using Bayesian inference. The dense gas considered is D6 siloxane and the equation of state used is the polytropic van der Waals equation. The shock tube data comes from the flexible asymmetric shock tube (FAST) experiment. This is modeled using the quasi-one-dimensional Euler equations with a source term that depends on time. A surrogate model based on sparse grids and a sensitivity analysis using Sobol' indices are both applied. The Markov chain Monte Carlo technique is applied to sample from the posterior probability distribution on the chosen parameters of the computer model. The results indicated that some of the thermodynamic parameters were identified, but that their mean values showed a disagreement with the true values in the literature. ...
Master thesis (2018) - Sanjeet Desai, Okko Boelens, Marc Gerritsma, Leon Rietveld, Stefan Hickel, Mirjam Snellen
Feadship yachts are designed for the leisure and cruising across the oceans. These luxury yachts are mostly powered by diesel engines or in some cases, a diesel-hybrid system. To prevent the inconvenience and the discomfort arising from the diesel exhaust gases for passengers, Feadship yachts are equipped with an underwater exhaust outlet. These underwater exhaust outlets are located on the side of the hull close to the dynamic waterline. It consists of an external appendage called "scoop" which creates a low pressure region for the exhaust outlet. During recent sea trials of the Feadship yachts, undesirable variations in the exhaust back-pressure were observed at the underwater outlet. These undesirable variations led to a situation with either too high back-pressure or too low back-pressure. An excessive back-pressure will increase the fuel consumption and will damage the diesel engine. Contrary, an extremely low back-pressure will give a visible exhaust flow above water thereby discolouring the hull and contaminating the deck with exhaust gases and steam.

An ideal scoop design would substantially reduce the above described problems. To investigate the optimal design for a scoop, a numerical method will be used. The method applicable in this study will be the Mutiphase Flow models from the commercial Computational Fluid Dynamics (CFD) software called Star CCM+. A multiphase fluid interaction between the exhaust gases and sea water will be examined to find the physical phenomenon affecting the back-pressure at the underwater outlet. In return, this phenomenon will be useful for a thorough analysis of the different scoop designs and how this design could impact the back-pressure. Furthermore, the validation of the numerical method will be carried out against the data procured from sea trials of the yachts with current scoop design.

To conclude, design recommendations for an optimal scoop geometry will be provided such that it can reduce the excessive back pressure, have a low resistance and prevent the discolouring of the hull. ...
Master thesis (2017) - Floris van de Beek, Marios Kotsonis, Pénélope Leyland, Stefan Hickel, Steven Hulshoff, Valeria Gentile
A robust and flexible numerical framework was developed for modelling of dielectric-barrier discharge applicable to plasma actuated flow control. In OpenFOAM, a solver was built which solves convection-diffusion equations for ionized species coupled to Gauss’ law for electrostatics and includes a chemical model to account for the various chemical reactions taking place in non-equilibrium cold plasma. The framework relies on a predictor-corrector method and was verified and qualitatively validated in 1D using a volumetric cathode-fall glow discharge case. Quantitative validation in 1D and extension to 2D was recommended for future work as steps towards increased capability of the solver. ...