HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 1–83 (2019) hjic.mk.uni-pannon.hu Contents Investigations Into Succinic Acid Fermentation ÁRON NÉMETH 1–4 Production of the Enzyme Cyclodextrin Glycosyltransferase Using Different Fermentation Techniques RÉKA CZINKÓCZKY AND ÁRON NÉMETH 5–10 Analysis of Crude Oil in Terms of Fouling and Corrosion RICHÁRD KATONA, ROLAND LÓCSKAI, GERGŐ BÁTOR, ANTAL KRÓJER, AND TIBOR KOVÁCS 11–15 Immobilization of Laccase in Alginate Beads BÉLA VARGA, MÓNIKA MEICZINGER, AND VIOLA SOMOGYI 17–23 Multistage Gravity Beneficiation of Rutile in a Tar-Free Sand Residue RUKAYAT AKANDE AND ABRAHAM ADELEKE 25–30 GPU-Accelerated Simulation of a Rotary Valve by the Discrete Element Method BALÁZS FÜVESI AND ZSOLT ULBERT 31–42 Residence Time Distribution-Based Analysis of an Industrial-Scale Biogas Fermenter ANDRÁS TANKOVICS, DÁVID TAKÁCS, JUDIT SZENDEFY, BÉLA CSUKÁS, AND MÓNIKA VARGA 43–51 Programmable Process Structure Based Analysis of Hydrogen Supply Chains MÓNIKA VARGA, CATHERINE AZZARO-PANTEL, JOSÉ MANUEL FLORES-PEREZ, AND BÉLA CSUKÁS 53–61 Investigation of Mixing in Tanks of a Special Geometry BÁLINT LEVENTE TARCSAY, ATTILA EGEDY, JANKA BOBEK, AND DÓRA RIPPEL-PETHŐ 63–70 Assessment of the Ecotoxicity of Nanoplastics NÓRA KOVÁTS, BETTINA ECK-VARANKA, ZSÓFIA BÉKÉSSY, DORINA DIÓSI, KATALIN HUBAI, AND JÁNOS KORPONAI 71–75 Biocomponent-Based Quality Improvement Opportunities of Binders for Road Construction GRÉTA TÁLOSI, PÉTER GERGÓ, AND ANDRÁS HOLLÓ 77–83 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 1–4 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-13 INVESTIGATIONS INTO SUCCINIC ACID FERMENTATION ÁRON NÉMETH *1 1Department of Applied Biotechnology and Food Science, Budapest University of Technology and Economics, Műegyetem rkp. 3, Budapest, 1111, HUNGARY Succinic acid (SA) is an important chemical intermediate from which fine chemicals ( e.g. detergents), additives (for pharmaceuticals, food (taste), plant growth stimulants) as well as other important intermediates (maleic anhydride, suc- cinimide, 2-pyrrolidinone, dimethyl succinate) can be manufactured. Since SA is involved in the central metabolism of cells (in the tricarboxylic acid (TCA) cycle), it is a key player in the biochemistry of life, which has the potential of biotech- nological production. Since SA is formed in the “middle” of the TCA cycle it can be formed by both CO2 production and fixation. The significance of the latter is that the amount of the product can be controlled by the availability of CO2, since stoichiometrically one molecule of CO2 is fixed by one molecule of SA. In our studies of compositions of Actinobacillus succinogenes media, the role and effect of pH regulator compounds as well as the effect of an inert atmosphere were investigated in terms of the yield. Furthermore, in fermentation experiments, the application of higher sugar concentra- tions was also studied. On the basis of different fermentations, a neural network for modelling and describing how factors influence SA production was established. Keywords: succinic acid, Actinobacillus succinogenes, neural network, modelling 1. Introduction Succinic acid (SA), as an intermediate of the tricar- boxylic acid (TCA) cycle, plays an essential role in the metabolism of microorganisms. SA can be produced by many anaerobes or facultative anaerobes as a metabolic product, thus can be used as an important platform chem- ical, a precursor of many pharmaceuticals, feed additives, green solvents or biodegradable polymers. SA itself is a colourless, odourless and crystal-forming compound. Since this metabolite is bifunctional (pKa1 = 4.21, pKa2 = 5.72 [1]), it is very reactive so has many potential applications, e.g. it plays important roles in the synthesis of γ-butyrolactone, maleic anhydryde, succinimide, 1,4- butanediol, dimethyl succinate, succinonitrile and 1,4- diaminobutane. It is industrially produced, mainly syn- thetically, in a complex way from maleic anhydride found in crude oil, which is both economically and environ- mentally unfavourable [2]. Therefore, its biotechnologi- cal production is a current research topic to find an alter- native method to avoid the above-mentioned side effects [3]. Furthermore, a great advantage of the microbial pro- duction of SA is that one of the initial biochemical re- actions is the carboxylation of phosphoenolpyruvate [4] which is regulated in the case of anaerobic bacteria by the availability of CO2 [2], thus the elevation of the CO2 con- centration can shift product portfolio from formate and ethanol towards SA (Fig. 1). *Correspondence: naron@f-labor.mkt.bme.hu Despite this fact, the process can help to de- crease the CO2 emissions of the human population [5]. Both fungi and bacteria can be found among the SA- producing microorganisms, but their ability to produce SA differs significantly: fungi ∼45 g/L (Aspergillus niger, rec. Yarrowia lipolytica), Gram-negative bacte- ria (wild-type Actinobacillus succinogenes 98.7 g/L, Mannheimia succiniciproducens 90 g/L) and Gram- positive bacteria (Clostridium thermosuccinogenes 82.5 g/L, rec. Corynebacterium glutamicum 145 g/L) [5]. The most frequently applied strains in the industry are Aspergillus niger, Aspergillus fumigatus, Byssochlamys nivea, Lentinus degener, Paecilomyces variotii, Peni- cillium viniferum, Saccharomyces cerevisiae and Acti- nobacillus succinogenes [4]. The latter bacterium is one of the most prominent producer strains that is isolated from bovine rumen and has been identified as a member of the genus Pas- teurella, which is a facultative anaerobic, non-motile, Gram-negative pleomorphic, rod-shaped bacterium [2]. It has great potential in terms of SA production, because of its higher yield and wide range of applicable sub- strates, e.g. glucose, cellobiose, maltose, lactose, saccha- rose, fructose and sorbitol. Furthermore, it can tolerate high initial concentrations of glucose, therefore, is suit- able for simple batch fermentation instead of the more complex and costly fed-batch culture technique [5]. The most widely applied strains are Actinobacillus succino- https://doi.org/10.33927/hjic-2019-13 mailto:naron@f-labor.mkt.bme.hu 2 NÉMETH Figure 1: SA metabolism with CO2 regulation genes 130Z and its variants (FZ6, 9, 21, 45 and 53) [2], which can tolerate both high glucose and SA concentra- tions as well as achieve high SA yields [5]. The aim of this paper is to compare fermentations of Actinobacillus succinogenes under different experimental conditions, then – due to the many influential factors – to set up a neural network-based model which can be used to predict high SA titres. 2. Materials and Methods 2.1 Cultivation of the bacteria Actinobacillus succinogenes 130Z (DSM22257) was cul- tivated in 10 ml of tryptic soy broth (TSB) (Sigma) in an impedimetric BacTrac 4100 (SY-LAB, Austria) anaero- bic cell. For a 1 L fermentation with a working volume of 0.8 L, an AS medium was applied according to Liu et al. [5] as follows: 62.5 g/L total sugar including 44.9 g/L saccharose, 9.8 g/L glucose and 7.2 g/L fructose sup- plemented with 15 g/L yeast extract, 1.5 g/L NaHPO4, 1g/L Na2HPO4, 1 g/L NaCl, 0.2 g/L MgCl2 and 0.2 g/L CaCl2. During preliminary tests, a temperature of 37◦C was not successful for SA production, therefore, 34◦C was applied and the pH regulated by 3M Na2CO3. Af- ter each fermentation, 200 ml of broth remained in the reactor and 600 ml was extracted, while 600 ml of fresh AS media was introduced. The 5 different fermentations are compared in the Results section. An innovative solution was to apply a CO2 economi- cal gas supply through the oxygen enrichment system of a Biostat Q DCU3 fermentor system in the absence of air and oxygen. The gas mix oxygen enrichment regulator was set at 4 %, therefore, periodically a small amount of CO2 was introduced into the fermentation broth. 2.2 Analysis During fermentations, samples were taken periodically and the optical density (OD) determined by a spectropho- tometer (Ultrospec Plus, Pharmacia LKB) at a wave- length of 600 nm (OD600) against the supernatant of a centrifuged sample by applying the same dilution factor (5×) as in the case of the samples. The cell dry weight (CDW) was obtained by a multiplication factor of 2 from OD600. Substrate consumption and product as well as by-product formation were detected by the Waters Breeze HPLC System by applying 5 mM H2SO4 in deionized water at a flow rate of 0.5 ml/min through a BioRad Aminex HPX87H column at 65◦C in a refractive index (RI) detector at 40◦C. 2.3 Neural networking For model building and evaluation, Neural Designer v2.9.5 was used by applying the following 4 steps: 1. Fermentation data were combined into a single MS Excel spreadsheet and exported to a tab-delimited text file, which could be imported into the modelling software. 9 variables, i.e. 7 inputs (time, lactic acid, acetic acid, propionic acid, glycerol, ethanol, total sugar) and 2 outputs (succinic acid, CDW), were applied to 58 fermentation samples from which 36 were used for training the network, 11 were selected and 11 were used for testing the behaviour of the model. 2. From among the many options, a model was defined (Fig. 2) by automatic scaling, without any principle component, with 2 layers and 3 neurons/hidden lay- ers that exhibit a logistic activation function in the absence of a bounding layer. In terms of a training strategy, the normalized mean squared error method was selected using a Quasi-Newton algorithm and a maximum of 1000 iterations. Incremental order was chosen as the order selection algorithm together with the growing inputs. 3. Model fit, i.e. performance training, was conducted. 4. Output of the model: impact figures of factors were determined (the other parameters were fixed), model equations obtained and predictions made by imple- menting input data into an input data matrix. 3. Results and Discussion Since Actinobacillus succinogenes is a facultative anaer- obic microorganism, the first fermentation experiment (Fig. 3A) was started in the absence of any specific at- mosphere. However, it ran very slowly, therefore, around 48 h (denoted by a red arrow) of continuous 5 % CO2 enrich- ment was applied via a zero flow rate gas inlet. The exper- iment confirmed that the application of CO2 is essential to form SA, therefore, finally 6 g/L SA was achieved and the model fitted very well for both CDW and SA mea- surements. Hungarian Journal of Industry and Chemistry INVESTIGATIONS INTO SUCCINIC ACID FERMENTATION 3 Figure 2: The used neural network for describing SA fer- mentations Therefore, the next fermentation (Fig. 3B) was con- ducted under 5 % CO2 enrichment of zero flow rate gas inlet and resulted in the highest 16 g/L SA (besides 11 g/L residual sugar) corresponding to a yield of 37 %. To increase the economic feasibility, the amount of CO2 was reduced by 50 % to 2.5 % in the following experi- ment (Fig. 3C). To avoid a fall in the SA concentration, detected by-products (such as lactic acid 3.5 g/L, acetic acid 5.7 g/L, propionic acid 1.7 g/L and glycerol 1.3 g/L) should be repelled as a result of the addition of 20 g/L calcium lactate. The model once again fitted well to the measured CDW and SA values. Unfortunately, the carbon flux shifted towards propionic acid formation, therefore, the subsequent experiment (Fig. 3D) was both supplemented with 20 g/L lactic acid and 11 g/L propionic acid, but this resulted in a low level of SA and a high level of lactic acid (52 g/L). The model yet again fitted well to the measured CSW and SA values. Independent fermentation results were also checked by the model (Fig. 3E) and resulted in very good fits. It can be concluded that the artificial neural network model constructed described well the SA fermentations, which is in line with the results of Li et al. [7], namely that ar- tificial neural network models can describe succinic acid fermentation better than response surface methodology. The presented results revealed that for SA fermenta- tion, 5% CO2 enrichment is essential and cannot be fully or partially replaced with the addition of lactic acid or propionic acid. After validating the model by conducting 5 different fermentations, it was used to optimize influential factors via impact figures (directional output plots) (Fig. 4). These trends show that all the presented factors cor- relate with SA concentration, i.e. any increment in their amounts resulted in an increment in SA, with the excep- tion of propionic acid which exhibited a negative corre- lation. These suggest that the addition of propionic acid can decrease the concentration of SA obtained and the addition of lactic acid can increase it. (A) (B) (C) (D) (E) Figure 3: Actinobacillus succinogenes fermentations with different degrees of CO2 enrichment in combination with lactic acid (LA) and propionic acid (PA). Dots indicate the measured values and lines indicate the model prediction. 47(2) pp. 1–4 (2019) 4 NÉMETH Figure 4: Factors impact plots: PrOH – propionic acid (g/L), Total sugar (g/L), LA – lactic acid (g/L) and CO2 – carbon dioxide (%) 4. Conclusion Five Actinobacillus succinogenes fermentations were run under different conditions: semi-anaerobic method, con- trolled introduction of CO2, and introduction of CO2 combined with either the addition of lactic acid or both lactic acid and propionic acid. All together 58 samples were taken and analysed, the results of which were en- tered into Neural Designer modelling software and used for training and testing the model. While the fermenta- tions resulted in very different final SA concentrations, the established model fitted well to all of the fermen- tations, even to the one which was not used for model building, testing and validation. While economical CO2 enrichment was successfully applied and resulted in the highest SA yield (37 %), the addition of lactic acid and propionic acid was not successful in terms of SA concen- tration. Acknowledgement This research was conducted within the framework of and in cooperation with PROGRESSIO Engineering Bureau Ltd., MÉL Biotech K+F Kft. and Budapest University of Technology and Economics. REFERENCES [1] Tonova, K.: State-of-the-art recovery of fermenta- tive organic acids by ionic liquids: An overview Hung. J. Ind. Chem. 2017 45(2), 41–44 DOI: 10.1515/hjic-2017-0019 [2] Song, H.; Lee, S. Y.: Production of succinic acid by bacterial fermentation Enzyme Microb. Tech. 2006 39(3), 352–361 DOI: 10.1016/j.enzmictec.2005.11.043 [3] Sauer, M.; Porro, D.; Mattanovich, D.; Branduardi, P.: Microbial production of organic acids: Expand- ing the markets Trends Biotechnol. 2008 26(2), 100–108 DOI: 10.1016/j.tibtech.2007.11.006 [4] Zeikus, J. G.; Jain,·M. K.; Elankovan, P.: Biotech- nology of succinic acid production and markets for derived industrial products Appl. Microbiol Biot. 1999 51(5), 545–552 DOI: 10.1007/s002530051431 [5] Liu, Y.-P.; Zheng, P.; Sun, Z.-H.; Ni, Y.; Dong, J.- J.; Zhu, L.-L.: Economical succinic acid production from cane molasses by Actinobacillus succinogenes Bioresource Technol. 2008 99(6), 1736–1742 DOI: 10.1016/j.biortech.2007.03.044 [6] Mika, L.T.; Cséfalvay, E.; Németh, Á.: Catalytic conversion of carbohydrates to initial chemicals: Chemistry and sustainability Chem. Rev. 2018 118(2), 505–613 DOI: 10.1021/acs.chemrev.7b00395 [7] Li, X.; Jiang, S.; Pan, L.; Wei, Z.: Optimization for the bioconversion of succinic acid based on re- sponse surface methodology and back-propagation artificial neural network in: Wang, H; Low, K. S.; Wei, K.; Sun, J. (eds.): Fifth International Con- ference on Natural Computation 2009 3, 392– 398 (IEEE Computer Society, Los Alamitos, USA) ISBN: 978-0-7695-3736-8 DOI: 10.1109/ICNC.2009.20 Hungarian Journal of Industry and Chemistry https://doi.org/10.1515/hjic-2017-0019 https://doi.org/10.1515/hjic-2017-0019 https://doi.org/10.1016/j.enzmictec.2005.11.043 https://doi.org/10.1016/j.tibtech.2007.11.006 https://doi.org/10.1007/s002530051431 https://doi.org/10.1016/j.biortech.2007.03.044 https://doi.org/10.1016/j.biortech.2007.03.044 https://doi.org/10.1021/acs.chemrev.7b00395 https://doi.org/10.1109/ICNC.2009.20 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 5–10 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-14 PRODUCTION OF THE ENZYME CYCLODEXTRIN GLYCOSYLTRANS- FERASE USING DIFFERENT FERMENTATION TECHNIQUES RÉKA CZINKÓCZKY1 AND ÁRON NÉMETH *1 1Department of Applied Biotechnology and Food Science, Budapest University of Technology and Economics, Műegyetem rkp. 3, Budapest, 1111, HUNGARY Cyclodextrins produced by cyclodextrin glycosyltransferase (CGTase) are widely used in the pharmaceutical industry to improve the solubility of drug substances as well as protect them against oxidation. The use of this enzyme in the cos- metics industry is also significant. CGTase is an enzyme that belongs to the α-amylase family, which is part of the group of non-Leloir glycosyltransferases. Enzyme-catalysed transglycosylation reactions may involve cyclization, coupling and disproportionation processes. The enzyme CGTase is mostly used to produce cyclodextrins (CDs). CGTase can produce α-, β- and γ-CDs during transglycosylation reactions, depending on the number of glucopyranose units involved (6, 7 or 8). The enzyme CGTase can also be used for enzymatic bioconversion, e.g., in the development of alternative sweeten- ers, where the bitter aftertaste of the product is reduced during the enzymatic bioconversion of steviol glycosides, thereby obtaining an even sweeter and more advantageous material. In our research, the enzyme CGTase was produced using different fermentation techniques to compare the activity and amount of CGTase produced by each process and optimize the subsequently planned scale-up. In our studies, the strain DSM 13 of Bacillus licheniformis was used, which produced CGTase extracellularly. During the experiments the batch, fed-batch and semi-continuous fermentation techniques were compared in terms of enzymatic production. All cultivation processes were carried out in a desktop lab scale fermenter. Keywords: Cyclodextrin glycosyltransferase, fermentation, cyclodextrins, Bacillus licheniformis 1. Introduction Cyclodextrin glycosyltransferase (CGTase, EC 2.4.1.19) is a starch-degrading enzyme, which is a member of the α-amylase family. The formal name of CGTase is [1,4-α- D-glucan 4-α-D-(1,4-α-glucano)-transferase(cyclizing). Kuriki et al. [1] reported that CGTase has the same four highly conserved regions as the α-amylases. CGTase catalyses four kinds of transglycosylation reactions (Fig. 1): cyclization, coupling, disproportionation and hydrol- ysis. These reactions are all transglycosylations, in which cyclization is intramolecular, coupling and disproportion- ation are intermolecular, and hydrolysis is the conversion of sugar to H2O [2]. The formation of cyclodextrin (CD) by the enzyme CGTase is an intermolecular transglyco- sylation reaction [3]. Many microorganisms are capable of producing CG- Tase, e.g., Bacillus macerans [4, 5], Bacillus amyloliq- uefaciens [6], Bacillus clarkii [7], Bacillus megaterium [8], Bacillus subtilis [9], Bacillus licheniformis [10, 11], Bacillus firmus [12,13], Bacillus circulans [14,15], Bacil- lus ohbensis [16, 17], Geobacillus stearothermophilus [18], Thermoanaerobacter sp. [19], Klebsiella pneumo- niae, and Klebsiella oxytoca [20]. *Correspondence: naron@f-labor.mkt.bme.hu Figure 1: CGTase-catalysed reactions: A: cyclization, B: coupling, C: disproportionation, D: hydrolysis https://doi.org/10.33927/hjic-2019-14 mailto:naron@f-labor.mkt.bme.hu 6 CZINKÓCZKY AND NÉMETH Figure 2: General structure of cyclodextrins (n: number of glucopyranose units, n = 6 α-CD, n = 7 β-CD, n = 8 γ-CD) (Figure adapted from: https://commons.wikimedia.org/wiki/File:Cyclodextrin.svg CC BY-SA 3.0) (08.07.2019) ) The molecular weight of CGTases may vary from 60 to 110 kDa, typically its proteins have a mass of 75 kDa [21]. The most important demand of metal ions for them is Ca2+, which protects the protein against heat denatura- tion. Most CGTases are strongly inhibited by Zn2+, Cu2+ and Fe2+ [22]. Cyclodextrins, produced from starch or its deriva- tives via enzymatic conversion, proceed through an in- tramolecular transglycosylation reaction using CGTases and to a lesser extent α-amylases [3]. They are cyclic oligosaccharides composed of α-1,4-glycosidic-linked glucosyl residues [23]. Three different types of cyclodex- trins exist and are characterised according to the num- ber of glucosyl residues in the molecule: α-, β- and γ- cyclodextrins consist of 6, 7 and 8 glucose units, respec- tively (Fig. 2). Cyclodextrins are cyclic molecules with a hydrophilic exterior and a hydrophobic cavity that en- ables them to form specific inclusion complexes with small hydrophobic molecules [24]. Cyclodextrins are chi- ral non-reducing oligosaccharides. Glucose is the decom- position product of all cyclodextrins in acidic solutions. The rate of hydrolysis follows the order of γ > β > α. Under acidic conditions, cyclodextrins are more slowly hydrolyzed than maltooligosaccharides. The gly- cosidic bonds in the cyclodextrins can be hydrolyzed by α-amylase, but β-amylase is unable to perform this hy- drolysis. The rate of enzymatic hydrolysis is the fastest for γ-CD, followed by β- then α-CD. All CDs are very stable and soluble in alkaline solutions at high pH. CDs are more resistant to acid or alkaline degradation than starch. CDs do not even degradate at temperatures as high as that of caramelization (> 200 ◦C, sterilization) under both dry or aqueous conditions of between pH 2 and 12. They are also stable up to 250 ◦C under an inert atmo- sphere of, for example, nitrogen [20, 25, 26]. The widespread use of cyclodextrins is due to their specific structure. Since each guest molecule is uniquely surrounded by the CD (or one of its derivatives), it is microencapsulated from a molecular microscopic point of view. This can result in beneficial changes to the chemical and physical properties of guest molecules, e.g., light- or oxygen-sensitive materials can be stabilized; very volatile substances fixed; the chemical reactivity of molecules modified; the solubility of materials improved; changes between phases achieved from liquid substances into powders; degradation of microorganisms avoided; bad smells and tastes masked; and pigments or colors of materials coated. As a result of these characteristics, CDs (and their derivatives) can be used in analytical chemistry, agriculture, the pharmaceutical as well as food industries and other masking areas. CGTase can be used for the transglycosylation of stevioside to rebaudioside through which the edulcorant quality can also be improved by increasing the substitution of steviol glycoside with the help of cornstarch hydrolyzate and CGTase [27–30]. 2. Materials and Methods 2.1 Cultivation of the bacteria The applied bacterial strain was Bacillus licheniformis B.01470 (DSM 13) purchased from the National Col- lection of Agricultural and Industrial Microorganisms in Hungary. Bacillus licheniformis is a Gram-positive, rod- shaped, endospore forming, facultatively anaerobic bac- teria. Nutrient agar was used to maintain the bacterium in Petri dishes [31]. In our research, three types of fermenta- tion techniques for the production of CGTase were com- pared: batch, fed-batch, and semi-continuous fermenta- tion techniques. All cultivation processes were carried out in a benchtop lab scale fermenter (Fig. 3) (Biostat Q, B. Braun Biotech International, Germany). In the fed-batch fermentation, after 24 hours 15 v/v % of the medium was fed into the bioreactor. During the semi-continuous fer- mentation at the end of each cycle, 80 % of the broth was replaced by fresh media. For the experiments in the bioreactor, Horikoshi II medium was used for the cultivation of bacteria contain- ing 1.0 % soluble starch, 0.5 % peptone, 0.5 % yeast ex- tract, 0.1 % K2HPO4, 0.02 % MgSO4 •7 H2O, and 1.0 % Na2CO3 (all concentrations are given in w/v in distilled water) [32]. Hungarian Journal of Industry and Chemistry https://creativecommons.org/licenses/by-sa/3.0/deed.en PRODUCTION OF THE ENZYME CYCLODEXTRIN GLYCOSYLTRANSFERASE 7 Figure 3: The bioreactor used in the experiments 2.2 The modelling of microbial growth In order to monitor the growth of bacterial cells, samples were taken during fermentations and the optical density (OD) measured at 600 nm. Microbial growth is described by µx = 1 x dx dt , (1) where µx is the specific growth rate of the microbe. It was evaluated through fitting the generalized logistic function Z = Zmax 1 + exp(a+ bt+ ct2 + dt3) (2) to the measured cell dry weight (CDW) values (calculated from at OD600). To fit the curve, SigmaPlot Version 12.0 software was applied. If the coefficient of determination (R2) was not high enough, the last two members of the generalized logistic function were omitted resulting in Z = Zmax 1 + exp(a+ bt) , (3) Figure 4: Colorimetric analysis of CGTase activity which also corresponds to the modified Monod model. The derivative of the fitted function is dZ dt = −Z ( 1− Z Zmax ) dH dt (4) where dH dt = b+ 2ct+ 3dt2 (5) is the derivative of the internal function. The auxiliary variable Z in Eq. 2, 3, and 4 was x (biomass in g/L), S (substrate in g/L), and Pi (product in g/L), respectively, while Zmax was xmax, S0 and Pi,max, respectively. If the fit was successful, then, by using determined constants of the model, the velocities and specific growth rates could be calculated by derivation from Eq. 4. 2.3 Measurement of CGTase activity During the fermentations, samples were regularly ex- tracted into Eppendorf tubes, which were centrifuged at 12,000 rpm for 6 minutes, then the cell-free supernatants were used to determine the enzyme activity. The measurement of extracellular CGTase activity was adapted from the method of Kaneko et al. [33] with slight modifications (with a reduced concentration of phenolphtalein). The colorimetric reaction (Fig. 4) was measured by a spectrophotometer at 550 nm. The experiments were conducted in 15 ml centrifuge tubes in a water bath at 40 ◦C. First, 4.5 ml 50 mM Tris- HCl buffer (pH = 9) was added containing a 1 % (w/v) water-soluble starch suspension, then 0.5 ml of cell-free supernatant containing the extracellular CGTase enzyme was introduced and homogenized thoroughly with a vor- tex mixer. Then four 0.5 ml samples were taken from each tube which were boiled for 5 minutes to inactivate the en- zyme. The boiled samples were transferred into 2 ml cu- vettes that contained a staining solution (1.2 ml 0.06 mM phenolphthalein in 0.5 M Na2CO3 solution). Four ab- sorbances at 550 nm of a given assay were plotted against time and the gradient (mmol/min) converted into enzyme activity with the help of a molar extinction coefficient (32, 263 M−1cm−1) resulting in the CGTase activity in unit/ml supernatant. 47(2) pp. 5–10 (2019) 8 CZINKÓCZKY AND NÉMETH Figure 5: Microbial growth during the batch fermentation Figure 6: Microbial growth during the fed-batch fermen- tation 3. Results and Discussion 3.1 Batch fermentation Fig. 5 shows the microbial growth during the batch fer- mentation and also represents the changes in the specific growth rate. The maximum value of the specific growth rate was 0.53 1/h. At the end of the fermentation, the fi- nal activity of CGTase was 0.3 U/ml and the productivity was 11.8 mU/(ml h). 3.2 Fed-batch fermentation Fig. 6 represents microbial growth during the fed-batch fermentation. Unfortunately, due to poorly scheduled sampling, it was not possible to adjust the generalized lo- gistic equation, therefore, it was impossible to calculate the specific growth rate. A fresh medium of 15 % was in- jected after 24 hours. The final enzyme activity was 0.5 U/ml and the enzyme productivity was 12.3 mU/(ml h). 3.3 Semi-continuous fermentation The semi-continuous fermentation is shown in Fig. 7, which consisted of 3 cycles. The highest value of the maximum growth rate was during the first cycle (0.5 1/h). As the fermentation progressed, the maximum specific growth rate decreased. Figure 7: Optical density and specific growth rate changes during the semi-continuous fermentation Figure 8: CGTase activities at the end of each cycle Fig. 8 shows that the enzyme activity of CGTase in- creased as the fermentation progressed. 3.4 Comparison of the different fermentation techniques Table 1 summarizes the enzyme activities and productiv- ities achieved by each fermentation technique. The maxi- mum specific growth rates reached in the batch and semi- continuous fermentations were approximately the same, which is characteristic of when the microorganism can multiply. The enzyme activities at the end of the fermentations Table 1: Comparison between the results of the different fermentation techniques Type µmax [1/h] Final enzyme activity [U/ml] Productivity [mU/(ml h)] Batch 0.53 0.3 11.8 Fed-batch n.d. 0.5 12.3 Semi- continuous 0.50 2.4 29.95± 0.3 Hungarian Journal of Industry and Chemistry PRODUCTION OF THE ENZYME CYCLODEXTRIN GLYCOSYLTRANSFERASE 9 rose as the complexity of the fermentation technique in- creased. While the fed-batch fermentation elongated the declining phase of the microbial growth cycle, in the semi-continuous fermentation technique an attempt was made to operate in the exponential growth phase. It is as- sumed that this difference caused the higher activity and productivity in the case of the semi-continuous fermenta- tion. 4. Conclusion In our experiments, the effect of the fermentation tech- nique on the activity of the produced enzyme CGTase was investigated. There was no significant difference between the activities of the produced CGTase and productivities of the systems when the batch and fed-batch fermenta- tions were compared. In contrast, bacteria produced a much more active enzyme during the semi-continuous fermentation, moreover, the productivity of this system was also significantly higher than that of the other two fermentation techniques. From the results, it can be assumed that the microbes produce the enzyme during the exponential growth phase, since no significant difference was observed between the batch and fed-batch fermentations. Meanwhile, a repeated exponential growth phase resulted in a much higher activity and productivity. This suggests that CG- Tase production follows growth associated-type product formation. Acknowledgement The research was supported by the Gedeon Richter’s Talentum Foundation, founded by Gedeon Richter Plc. (Gedeon Richter PhD fellowship). REFERENCES [1] Kuriki, T.; Imanaka, T.: The concept of the α- amylase family: Structural similarity and common catalytic mechanism, J. Biosci. Bioeng. 1999 87(5), 557–565 DOI: 10.1016/S1389-1723(99)80114-5 [2] Kobayashi, S.: Cyclodextrin producing enzyme (CGTase) in Park, K.-H.; Robyt, J. F.; Choi, Y.- D. (eds.): Enzymes for carbohydrate engineering (Progress in Biotechnology) 1996 12 23–41 (Else- vier, Amsterdam, The Netherlands) ISBN: 978-0-444- 82408-0 DOI: 10.1016/S0921-0423(96)80360-1 [3] van der Veen, B. A.; van Alebeek, G.-J. W. M.; Uitdehaag, J. C. M.; Dijkstra, B. W.; Dijkhuizen, L.: The three transglycosylation reactions catalyzed by cyclodextrin glycosyltransferase from Bacillus circulans (strain 251) proceed via different kinetic mechanisms Eur. J. Biochem. 2000 267(3) 658–665 DOI: 10.1046/j.1432-1327.2000.01031.x [4] Isao, K.; Yoshida, N.: Method of producing β- cyclodextrin 1996 US 5556775, United States Patent and Trademark Office https://patents. google.com/patent/US5556775A/en [5] Shieh, W. J.; Hedges, A. R.: Process for producing alpha-cyclodextrin using cyclomaltodextrin glucan- otransferase in presence of cyclohexane 1994 US 5326701, United States Patent and Trademark Of- fice https://patents.google.com/patent/US5326701A/ en [6] Yu, E. K. C.; Aoki, H.; Misawa, M.: Specific alpha-cyclodextrin production by a novel ther- mostable cyclodextrin glycosyltransferase Appl. Microbiol. Biotechnol. 1988 28(4–5) 377–379 DOI: 10.1007/BF00268199 [7] Wu, D.; Chen, S.; Wang, N.; Chen, J.; Wu, J.: Gamma-cyclodextrin production using cyclodex- trin glycosyltransferase from Bacillus clarkii 7364 Appl. Biochem. Biotechnol. 2012 167(7) 1954–1962 DOI: 10.1007/s12010-012-9741-5 [8] Pishtiyski, I.; Popova, V.; Zhekova, B; Characteri- zation of cyclodextrin glucanotransferase produced by Bacillus megaterium Appl. Biochem. Biotechnol. 2008 144 263–272 DOI: 10.1007/s12010-007-8009-y [9] Cheirsilp, B.; Kitcha, S.; Maneerat, S.: Kinetic characteristics of β-cyclodextrin production by cy- clodextrin glycosyltransferase from newly isolated Bacillus sp. C26 Electron. J. Biotechnol. 2010 13(4) DOI: 10.2225/vol13-issue4-fulltext-6 [10] Bonilha, P. R. M.; Menocci, V.; Goulart, A. J.; Polizeli, M. L. T. M.; Monti, R.: Cyclodextrin gly- cosyltransferase from Bacillus licheniformis: Opti- mization of production and its properties Braz. J. Microbiol. 2006 37(3) 317–323 DOI: 10.1590/S1517- 83822006000300022 [11] Thombre, R. S.; Kanekar, P. P.: Synthesis of β- Cyclodextrin by Cyclodextrin glycosyl transferase produced by Bacillus licheniformis MCM–B1010 J. Microbiol. Biotech. Res. 2013 3(1) 57–60 [12] Matioli, G.; Zanin, G. M.; De Moraes, F. F.: In- fluence of substrate and product concentrations on the production of cyclodextrins by cgtase of Bacil- lus firmus, strain No. 37 Appl. Biochem. Biotechnol. 2002 98(1–9) 947–961 DOI: 10.1385/ABAB:98-100:1- 9:947 [13] Gawande, B. N.; Goel, A.; Patkar, A. Y.; Nene, S. N.: Purifcation and properties of a novel raw starch degrading cyclomaltodextrin glucanotransferase from Bacillus firmus Appl. Microbiol. Biotechnol. 1999 51 504–509 DOI: 10.1007/s002530051424 [14] Pinto, F. S. T.; Flôres, S. H.; Ayub, M. A. Z.; Hertz, P. F.: Production of cyclodextrin glycosyltransferase by alkaliphilic Bacillus circulans in submerged and solid-state cultivation Bioprocess Biosyst. Eng. 2007 30(5) 377–382 DOI: 10.1007/s00449-007-0134-z [15] Iyer, J. L.; Shetty, P.; Pai, J. S.: Immobilisa- tion of cyclodextrin glucanotransferase from Bacil- lus circulans ATCC 21783 on purified seasand J. Ind. Microbiol. Biotechnol. 2003 30(1) 47–51 DOI: 10.1007/s10295-002-0009-x 47(2) pp. 5–10 (2019) https://doi.org/10.1016/S1389-1723(99)80114-5 https://doi.org/10.1016/S0921-0423(96)80360-1 https://doi.org/10.1046/j.1432-1327.2000.01031.x https://patents.google.com/patent/US5556775A/en https://patents.google.com/patent/US5556775A/en https://patents.google.com/patent/US5326701A/en https://patents.google.com/patent/US5326701A/en https://doi.org/10.1007/BF00268199 https://doi.org/10.1007/BF00268199 https://doi.org/10.1007/s12010-012-9741-5 https://doi.org/10.1007/s12010-007-8009-y https://doi.org/10.2225/vol13-issue4-fulltext-6 https://doi.org/10.1590/S1517-83822006000300022 https://doi.org/10.1590/S1517-83822006000300022 https://doi.org/10.1385/ABAB:98-100:1-9:947 https://doi.org/10.1385/ABAB:98-100:1-9:947 https://doi.org/10.1007/s002530051424 https://doi.org/10.1007/s002530051424 https://doi.org/10.1007/s00449-007-0134-z https://doi.org/10.1007/s10295-002-0009-x https://doi.org/10.1007/s10295-002-0009-x 10 CZINKÓCZKY AND NÉMETH [16] Nishida, T.; Nakamura, A.; Masaki, H.; Uozumi, T.: Regulation of cyclodextrin glucanotransferase synthesis in Bacillus ohbensis FEMS Microbiol. Lett. 1997 149(2) 221–226 DOI: 10.1111/j.1574- 6968.1997.tb10332.x [17] Cami, P. H.; Majou, D. B.: Process for production of cyclodextrins 1994 US 5376537, United States Patent and Trademark Office https://patents. google.com/patent/US5376537A/en [18] Shiosaka, M.: Heat stable cyclodextrin glycosyl- transferase 1976 US 3988206, United States Patent and Trademark Office https://patents.google.com/ patent/US3988206A/en [19] Martín, M. T.; Plou, F. J.; Alcade, M.; Balles- teros, A.: Immobilization on Eupergit C of cy- clodextrin glucosyltransferase ( CGTase ) and prop- erties of the immobilized biocatalyst J. Mol. Catal. B Enzym. 2003 21(4-6) 299–308 DOI: 10.1016/S1381- 1177(02)00264-3 [20] Wimmer, T.: Cyclodextrins in Ullmann’s Ency- clopedia of Industrial Chemistry 2012 11 DOI: 10.1002/14356007.e08_e02 [21] Uitdehaag, J. C. M.; Kalk, K. H.; van der Veen, B. A.; Dijkhuizen, L.; Dijkstra, B. W.: The cycliza- tion mechanism of cyclodextrin glycosyltransferase (CGTase) as revealed by a γ-cyclodextrin-CGTase complex at 1.8-Å resolution J. Biol. Chem. 1999 274(49) 34868–34876 DOI: 10.1074/jbc.274.49.34868 [22] Thangadurai, D.; Sangeetha, J. (eds.): Industrial Biotechnology, Sustainable Production and Biore- source Utilization 2017 (Apple Academic Press, Oakville, Canada), ISBN: 978-1-771-88262-0 [23] Biwer, A.; Antranikian, G.; Heinzle, E.: Enzy- matic production of cyclodextrins Appl. Microbiol. Biotechnol. 2002 59(6) 609–617 DOI: 10.1007/s00253- 002-1057-x [24] Szejtli, J.: Introduction and general overview of cy- clodextrin chemistry Chem. Rev. 1998 98(5) 1743– 1753 DOI: 10.1021/cr970022c [25] Saenger, W.: Cyclodextrin inclusion compounds in research and industry Angew. Chem. Int. Ed. Engl. 1980 19(5) 344–362 DOI: 10.1002/anie.198003441 [26] Szejtli, J.: The cyclodextrins and their applications in biotechnology Carbohydr. Polym. 1990 12(4) 375–392 DOI: 10.1016/0144-8617(90)90088-A [27] Del Valle, E. M. M.: Cyclodextrins and their uses: a review Process Biochem. 2004 39(9) 1033–1046 DOI: 10.1016/S0032-9592(03)00258-9 [28] Singh, M.; Sharma, R.; Banerjee, U. C.: Biotech- nological applications of cyclodextrins Biotech- nol. Adv. 2002 20(5–6) 341–359 DOI: 10.1016/S0734- 9750(02)00020-4 [29] Li, S.; Li, W.; Xiao, Q. Y.; Xia, Y.: Transglycosyla- tion of stevioside to improve the edulcorant quality by lower substitution using cornstarch hydrolyzate and CGTase Food Chem. 2013 138(2–3) 2064–2069 DOI: 10.1016/j.foodchem.2012.10.124 [30] Kochikyan, V. T.; Markosyan, A. A.; Abelyan, L. A.; Balayan A. M.; Abelyan, V. A.: Combined en- zymatic modification of stevioside and rebaudioside A Appl. Biochem. Microbiol. 2006 42(1) 37–43 DOI: 10.1134/S0003683806010030 [31] DSMZ Nutrient agar (n.d.) Retrieved from https://www.dsmz.de/microorganisms/medium/pdf/ DSMZ_Medium1.pdf, date accessed: 2019.10 [32] Park, C. S.; Park, K. H.; Kim, S. H.: A rapid screen- ing method for alkaline β-cyclodextrin glucan- otransferase using phenolphthalein-methyl orange- containing solid medium Agric. Biol. Chem. 1989 53(4) 1167–1169 DOI: 10.1080/00021369.1989.10869443 [33] Kaneko, T.; Kato, T.; Nakamura, N.; Horikoshi, K.: Spectrophotmetric determination of cyclization ac- tivity of β-cyclodextrin-forming cyclomaltodextrin glucanotransferase J. Jpn. Soc. Starch Sci. 1987 34(1) 45–48 DOI: 10.5458/jag1972.34.45 Hungarian Journal of Industry and Chemistry https://doi.org/10.1111/j.1574-6968.1997.tb10332.x https://doi.org/10.1111/j.1574-6968.1997.tb10332.x https://patents.google.com/patent/US5376537A/en https://patents.google.com/patent/US5376537A/en https://patents.google.com/patent/US3988206A/en https://patents.google.com/patent/US3988206A/en https://doi.org/10.1016/S1381-1177(02)00264-3 https://doi.org/10.1016/S1381-1177(02)00264-3 https://doi.org/10.1002/14356007.e08_e02 https://doi.org/10.1002/14356007.e08_e02 https://doi.org/10.1074/jbc.274.49.34868 https://doi.org/10.1007/s00253-002-1057-x https://doi.org/10.1007/s00253-002-1057-x https://doi.org/10.1021/cr970022c https://doi.org/10.1002/anie.198003441 https://doi.org/10.1016/0144-8617(90)90088-A https://doi.org/10.1016/S0032-9592(03)00258-9 https://doi.org/10.1016/S0734-9750(02)00020-4 https://doi.org/10.1016/S0734-9750(02)00020-4 https://doi.org/10.1016/j.foodchem.2012.10.124 https://doi.org/10.1134/S0003683806010030 https://doi.org/10.1134/S0003683806010030 https://www.dsmz.de/microorganisms/medium/pdf/DSMZ_Medium1.pdf https://www.dsmz.de/microorganisms/medium/pdf/DSMZ_Medium1.pdf https://doi.org/10.1080/00021369.1989.10869443 https://doi.org/10.5458/jag1972.34.45 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 11–15 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-15 ANALYSIS OF CRUDE OIL IN TERMS OF FOULING AND CORROSION RICHÁRD KATONA1, ROLAND LÓCSKAI1,2, GERGŐ BÁTOR1,2, ANTAL KRÓJER1, AND TIBOR KOVÁCS*1 1Institute of Radiochemistry and Radioecology, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY 2Social Organisation for Radioecological Cleanliness, József Attila u. 7/A, Veszprém, 8200, HUNGARY The present-day oil refining market has the potential to process opportunity crudes. The degree of corrosion and fouling issues has increased in terms of processing opportunity crudes. In the case of changing crude slates, predicting the impacts of crude oil on production is necessary to ensure safe and profitable processing. Crude oil can be characterized by the quantitative and qualitative analyses of saturates, aromatics, resins, asphaltenes and contaminants as well as the determination of its physical and chemical properties. The complementary nature of these analyses is necessary to identify the root causes of fouling and corrosion problems and quantify the impact of them. In this study, analytical methods were developed to investigate the effect of crude oils on production in terms of fouling and corrosion, which can be used to characterize crude oil from different sources. The compatibility, emulsification tendency and fouling potential of crude oil were analyzed by the developed analytical methods. The measurement method to assess the compatibility of crude oils is based on the determination of the flocculation tendency of asphaltenes in oils. A Porla GLX Step Analyzer was used for the compatibility tests. Emulsification and fouling tendencies were measured by our in-house designed laboratory equipment. These measurements were conducted under typical operating conditions. The developed analytical methods were validated by crude oils from five different sources. Keywords: crude oil, fouling, corrosion, emulsion, desalter, compatibility, asphaltenes 1. Introduction Refineries have had to handle problems associated with fouling and corrosion for a prolonged period of time. The degree of corrosion and fouling issues has increased since the processing of opportunity crudes and co-processing of renewable feeds were introduced [1, 2]. It is well known that crude oils contain water, chloride salts and organochlorine compounds [3]. When crude oils are preheated in the heat exchangers of crude units, the hydrolysis of chloride salts and organochlorine com- pounds occurs. In these chemical reactions, gaseous hy- drochloric acid is released. The gaseous hydrochloric acid present in the steam phase can easily dissolve in condensed water, resulting in very low pH values. The concentrated caustic solution can cause general corro- sion of equipment in crude units composed of carbon steel [4, 5]. Desalting is the first line of defense against corrosion induced by the hydrochloric acid [6]. Earlier studies showed that the efficiency of desalter apparatus is influenced by the stability of the emulsion generated. The stability of the emulsion increases and the efficiency of the desalter decreases when the asphaltene content of the crude oil is high [7]. Besides the desalting process, *Correspondence: kt@almos.uni-pannon.hu neutralizers and corrosion inhibitors are also used against corrosion in crude units [8]. Fouling deposits can be categorized into inorganic and organic types. Organic fouling results from as- phaltenes and high molecular weight hydrocarbons which become insoluble in the system [9]. The asphaltene molecules can become unstable due to the effect of blend- ing and heating incompatible crude oils. The precipitated molecules can be deposited on the surface of heat ex- changers [10]. Corrosion products of iron (FeS, Fe2O3) and inorganic contaminants in crude oils can cause inor- ganic fouling. In the crude oil preheat train, the viscosity of the oil is lowered because of heating and the deposi- tion of inorganic contaminants increases. Regardless of the cause, high fouling rates can lead to excessive equip- ment cleaning requirements and costs. Chemicals–such as dispersants and inhibitors–are used to effectively re- duce fouling in several critical areas throughout a refinery [11]. In the case of changing crude slates, predicting the impacts of crude oil on production is necessary for safe and profitable processing [12]. In this study, analytical methods were developed to investigate the effect of crude oils on production with regard to fouling and corrosion, which can be used for characterizing crude oil from dif- https://doi.org/10.33927/hjic-2019-15 mailto: kt@almos.uni-pannon.hu 12 KATONA, LÓCSKAI, BÁTOR, KRÓJER, AND KOVÁCS Figure 1: Schematic drawing of laboratory desalter model equipment. ferent sources. The key causes of crude oil processing is- sues are incompatibility, the emulsification tendency and the fouling potential [3]. 2. Experimental Five crude oils from different sources were used in this study to validate the developed procedure. API gravity, Total acid number (TAN), water content, asphaltene con- tent, inorganic contaminants and the distillation curve of crude oils were measured in accordance with ASTM D4052, ASTM D664, ASTM D4807, MSZ EN 459, MSZ EN ISO 9029 and ASTM D7096. In addition, the total chlorine content of one desalted crude oil was measured by neutron activation analysis. 2.1 Emulsification tendency In the first section of the determination of the emulsifi- cation tendency, 950 ml of crude oil and 50 ml of water were agitated using an ULTRA-TURRAX stirrer at 75 °C and 13, 500 rpm for 5 minutes. Then 80 ml of emul- sion was poured out into a tube that was made according to the standard test method ASTM D96. After that the emulsification tendencies of crude oils were analyzed by our in-house-designed laboratory desalter simulator. The schematic drawing of the equipment used is shown in Fig. 1. The desalter simulator apparatus consists of a ther- moblock, 8 transformers and toroidal transformers. The emulsion separation tubes are located in the aluminium thermoblock. The electric field in the tubes is controlled by the installed transformers. The connection between tubes and toroidal transformers is provided by special caps with two isolated electrodes. The temperature of the thermoblock can be adjusted by the control unit of the apparatus. During the measurement, the temperature was 110 °C and the electrical voltage was 2, 500 V. The emulsion was resolved in the separation tubes in accor- dance with the known method of providing an electric field to polarize water droplets. The volume of water that dropped into each tube was recorded after 10, 20, 30, and 60 minutes. The more water that dropped and the clearer the water phase, the more effective the desalting process [13]. 2.2 Fouling tendency The fouling tendency was measured in the high- performance thermoblock. The equipment is an alu- minium thermoblock in which 8 places for special tubes and an electric heater coil are located. The heating can be regulated by a control unit. The special tubes are pressure vessels 100 ml in volume. Into the tubes, 25 g of crude oil samples were measured. The samples were heated at 250 °C for 48 hours then cooled to 80 °C and filtered under a vacuum through a 0.7 µm membrane filter. The filter along with the residue was washed with petrol, dried and weighed to yield the final result. In order to measure the amount of generated deposits on the wall of the tubes, they were washed with petrol and then the petrol wash was filtered under a vacuum through a 0.7 µm membrane filter too. The deposit was calculated from S = 10000 (m2 −m1)100 25 , (1) where the mass of the filter is denoted by m1, the mass of the filter with the deposit by m2 and the amount of deposits in ppm by S [14]. 2.3 Neutron activation analysis Neutron activation analysis determines the total chlorine content in the sample regardless of its oxidation, inor- ganic or organic states. The measurement is capable of handling inhomogeneous samples of greater masses, thus mitigating the uncertainties of sampling. The measuring range spans several times the magnitude of the range of alternative techniques, and has a small measurement error [15]. Samples were prepared in 1 ml vials and sealed in polyethylene bags. The samples were then placed in larger 5 cm plastic vials that were irradiated. To pro- vide redundancy in the case of a leakage, the outer vial was heat sealed using a hot iron. The masses of the sam- ples were approximately 1 g and were irradiated by a 1.1 MW Mark II TRIGA Reactor at the University of Texas in Austin. In order to examine the total chlorine content, samples were irradiated using a reactor power of 500 kW for 5 minutes. Thermal neutrons were used due to favor- able cross-sections for radiative capture within this range. The following reaction was used to determine the total chlorine content: 37Cl(n, γ)38Cl with its half-life of 37 minutes. Samples were allowed to decay for 15 minutes prior to being counted by a gamma-ray spectrometer. The samples themselves were counted for between 600 and 2, 700 seconds, depending on the net count rate under the 1, 642 keV photopeak for chlorine. Hungarian Journal of Industry and Chemistry ANALYSIS OF CRUDE OIL IN TERMS OF FOULING AND CORROSION 13 Table 1: Summary of the properties of the crude oils Crude oil A B C API 29 46 42 TAN (mg KOH/g) 0.19 0.0448 0.0459 Water Content (%) 0.1 0 0 Asphaltene content (%) 2.42 0.017 0.25 Inorganic contamination (ppm) 260 110 82 Total chlorine content after desalting (ppm) 8.08± 0.18 Yield of fractions (%) C1-C10 18 40 32 C11-C14 9 12 12 C15-C25 23 27 26 C26-C50 25 16 19 C51+ 25 5 11 Emulsification tendency Fig. 2 S (ppm) 1286 472 248 SBN (-) 77.6 43.3 64 IN (-) 44.6 lower lower Crude oil D E API 29 31 TAN (mg KOH/g) 0.17 0.0726 Water Content (%) 0 0 Asphaltene content (%) 2.05 1.6 Inorganic contamination (ppm) 149 198 Yield of fractions (%) C1-C10 18 20 C11-C14 8 5 C15-C25 22 25 C26-C50 25 27 C51+ 27 23 Emulsification tendency Fig. 2 S (ppm) 1187 772 SBN (-) 79.5 68.4 IN (-) 30.8 25.1 2.4 Compatibility The compatibility test of crude oils is based on the Oil Compatibility Model. This model enables the insolubil- ity number (IN) and solubility blending number (SBN) of crude oils to be calculated. Two tests were carried out to measure the solubility parameter. The first test involved determining the maximum volume of n-Heptane that can be added to a given volume of oil without the precipita- tion of asphaltenes. Then the minimum percent of toluene in the mixture with n-Heptane to dissolve asphaltenes was determined. The volume percentage of toluene in the toluene–n-heptane mixture was plotted against the vol- ume ratio of oil to the toluene–n-heptane mixture. The insolubility number is where a line drawn through the two points intercepts the y-axis. The solubility blending num- ber can be calculated by SBN = IN [ 1 + ml heptane ml oil ] . (2) The criterion for the compatibility of the oil mixture is that the volume average SBN is greater than the IN of all types of oil in the mixture. It follows that the higher the SBN, the more compatible the crude oil is [16, 17]. A Porla GLX Step Analyzer was used to analyze the compatibility of crude oils. The laboratory analyzer per- forms the dilution in an aromatic solvent and the titration with paraffin hydrocarbon automatically. The analyzer was used in accordance with the factory-preset program. Insoluble asphaltenes were detected by an optical detec- tor according to the intensity of scattered visible light. The wavelength of the light was selected to be sensitive enough to detect even minute concentrations of asphal- tene particles in the measuring solution. 3. Results and Discussion The results of analytical measurements are summarized in Table 1. Samples B and C are light crude oils, while A, D and 47(2) pp. 11–15 (2019) 14 KATONA, LÓCSKAI, BÁTOR, KRÓJER, AND KOVÁCS Figure 2: The emulsification tendencies and the pictures of emulsions after the measurements. The order of the crude oils in the pictures from left to right is the following: Crude oil A, B, C, D and E. E are medium crude oils according to API gravity. The crude oils in order of asphaltene content from lowest to highest are: B, C, E, D and A. It is recognized that the highest emulsification of water occurs when the asphal- tene content is high [6]. In blueFig. 2, it can be seen that the crude oil with the highest content of asphaltenes (A) formed the most stable emulsions. On the other hand, the emulsification tendencies of crude oils B and C showed that the emulsion of light crude oil can be easily resolved. The quantities of water phases following the measure- ments are presented in Fig. 2. The water phases of stable emulsions are oily, while those of emulsions with lower asphaltene contents are clear. The total chlorine content of desalted crude oil was measured by neutron activation analysis. The total chlo- rine content in the desalted crude oil may be the result of the remaining dissolved salts and organic chlorides. Organic chlorides present a more significant problem in crude oil refineries when compared to inorganic chlo- rides, since they are not removed by the desalting process [18]. The fouling tendency can also be correlated with the asphaltene content. The lowest amount of deposit was measured for the crude oil of the lowest content of as- phaltenes (Crude oil B). From Table 1, it can be seen that the higher the asphaltene content of crude oil, the greater the amount of the deposit. However, the asphal- tene content of Crude oil C is higher than that of B, but the amount of the deposit of Crude oil C is lower than Figure 3: Compatibility test results of a mixture of Crudes B-D. Figure 4: Compatibility test results of a mixture of Crudes A-C. Figure 5: Compatibility test results of a mixture of Crudes C-E. that of B. This is explained by the amount of inorganic contaminants which is higher for Crude oil B than C. Furthermore, as the results of the compatibility test show, lower SBNs belong to the light crude oils with low asphaltene contents. In Figs. 3-5, the evaluation of three crude oil blends is shown. INmax denotes the greater IN of the two respective crude oils that were blended. In all three cases, the volume average SBN is always greater than the INmax, so the mixtures of Crude oils B-D, C-A and C-E are compatible with any blending ratios. The INs of Crude oils B and C were not determined because the software of the analyzer just shows the highest IN among the blended crude oils. 4. Conclusion Fouling and corrosion exact significant economic and operational penalties. The risk of prosecution has in- creased with the processing of opportunity crudes and co- processing of renewable feeds [19]. In this paper, some laboratory methods were developed to investigate the ef- fect of crude oils on production in terms of fouling and Hungarian Journal of Industry and Chemistry ANALYSIS OF CRUDE OIL IN TERMS OF FOULING AND CORROSION 15 corrosion. The emulsification tendency as well as the ef- fect of blending and heating on crude oils were mea- sured by the developed laboratory methods. These meth- ods can be used for characterizing crude oil from differ- ent sources. The results of developed laboratory methods were compared with the composition of crude oils. As is confirmed by the results, asphaltenes and high molecu- lar weight hydrocarbons have a significant impact on the production of crude oils in terms of fouling and corrosion. The cost and complexity of the determination of total chlorine content by neutron activation analysis make it unfit for routine analytical purposes, but the precision of the results makes it a perfect choice as a reference method [20]. Our plan is to propose methods to quantify lower amounts of total chlorine content in crude oils using neu- tron activation analysis as a reference method. REFERENCES [1] Ropital, F.: Current and future corrosion challenges for a reliable and sustainable development of the chemical, refinery, and petrochemical industries, Materials and Corrosion, 2009 60(7), 495–500 DOI: 10.1002/maco.200805171 [2] Gergely, A.; Locskai, R.; Szabó, P.; Krójer, A.; Kristóf, T.: Hydrogen sulphide corrosion of car- bon and stainless steel alloys immersed in mix- tures of renewable fuel sources and tested under co- processing conditions, Hung. J. Ind. Chem., 2016 44, 55–70 DOI: 10.1515/hjic-2016-0007 [3] Teran, C.: Crude fingerprinting and predictive analytics, Digital Refining, 2017 https://www. digitalrefining.com/article/1001416 [4] Clarida, D.; Johnston, J.; McConnell, M.; Strong, R.: Corrosion and fouling experi- ences in crude units using low base strength neutralizers. (Corrosion 97, NACE Interna- tional, Houston, 1997, NACE-97499) https: //www.onepetro.org/conference-paper/NACE-97499 [5] Al-Omari, A. S.; Al-Zahrani, A. M.; Lobley, G. L.; Tems, R. D.; Dias, O. C.: Case Stud- ies and Best Practices of Refinery Caustic In- jection Systems, Saudi Aramco Journal of Tech- nology, 2009, 29-36 https://nvlpubs.nist.gov/ nistpubs/Legacy/SP/nbsspecialpublication674.pdf [6] Gutzeit, J.: Controlling Crude Unit Overhead Cor- rosion – Rules of Thumb for Better Crude Desalt- ing. (Corrosion 2007, NACE International, Hous- ton, 2007, NACE-07567) https://www.onepetro. org/conference-paper/NACE-07567 [7] Fingas, M.: Water in Oil Emulsion Formation: A re- view of Physics and Mathematical modeling, Spill Sci. Technol. B., 1995 2(1), 55–59 DOI: 10.1016/1353- 2561(95)94483-Z [8] Lukács, Z.; Kristóf, T.: Application of 2N design of experiment method for the evaluation of the effi- ciency and cross-effects of oilfield chemicals, Hung. J. Ind. Chem., 2018 46(2), 13–17 DOI: 10.1515/hjic- 2018-0012 [9] Polley, G. T.; Wilson, B. D.; Yeap, B. L.; Pugh, S. J.: Evaluation of laboratory crude oil thresh- old fouling data for application to refinery pre-heat trains, Appl. Therm. Eng., 2002 22(7), 777–788 DOI: 10.1016/S1359-4311(02)00023-6 [10] Wiehe, A. I.: The chemistry of petroleum fouling, (Proceedings of the 4th Interna- tional Conference on Refinery Processing, AIChE, New York, 2001), pp. 204-210 http://www.solublesolutions.com/wp-content/ themes/twentyeleven_child/images/FoulChem2001.doc [11] Wilson, R. M.; Perugini, J.J.: Antifoulants: a proven energy-savings investment, (National Petroleum Re- finers Associaton, Technical Paper, United States, 1985) https://www.osti.gov/biblio/6358216 [12] Nagi-Hanspal, I.; Subramaniyam, M.; Shah, P.; Moretti, S.; Noland, J.: Exploiting opportunities with challenging crudes, (The Refining, Gas & and Petrochemicals Processing Website) 2012/Q4 http: //www.eptq.com/view_article.aspx?intAID=1197 [13] Hajivand, P.; Vaziri, A.: Optimization of demulsifier formulation for separation of water from crude oil emulsions, Braz. J. Chem. Eng., 2015 32 (1), 107– 118 DOI: 10.1590/0104-6632.20150321s00002755 [14] ASTM D4807(2015), Standard Test Method for Sediment in Crude Oil by Membrane Filtration, ASTM International, DOI: 10.1520/D4807-05R15 [15] Wiehe, A. I.; Kennedy, J. R.: Application of the Oil Compatibility Model to Refinery Streams, Energy Fuels, 2000 14(1), 56–59 DOI: 10.1021/ef9901342 [16] Rathore, V.; Brahma, R.; Thorat, S. T.: Assessment of crude oil blends, Digital Refining, 2011 https: //www.digitalrefining.com/article/1000381 [17] Pei, P. X.; Hsu, S.; Fleming, W. R.: Test methods for total chlorine in lubricating base oils, Measure- ment and Standards for Recycled Oil, 1982 4, 271– 280 https://nvlpubs.nist.gov/nistpubs/Legacy/SP/ nbsspecialpublication674.pdf [18] Li, X.; Wu, B.; Zhu, J.: Hazards of organic chloride to petroleum processing in Chinese Refineries and Industrial Countermeasures, Progress in Petrochemical Science, 2018 2(3) DOI: 10.31031/PPS.2018.02.000539 [19] Garrett, T.; Rattanakhambay, A.; Robbins, N.; Wunder, M.; Yeung, T.: The challenges of crude blending, Digital Refining, 2016 https://www. digitalrefining.com/article/1001216 [20] Tagami, K.; Uchida, S.; Hirai, S.; Tsukada, H.; Takeda, H.: Determination of chlorine, bromine and iodine in plant samples by inductively cou- pled plasma-mass spectrometry after leaching with tetramethyl ammonium hydroxide under a mild temperature condition, Anal. Chim. Acta,2006 570(1), 88–92 DOI: 10.1016/j.aca.2006.04.011 47(2) pp. 11–15 (2019) https://doi.org/10.1002/maco.200805171 https://doi.org/10.1002/maco.200805171 https://doi.org/10.1515/hjic-2016-0007 https://www.digitalrefining.com/article/1001416 https://www.digitalrefining.com/article/1001416 https://www.onepetro.org/conference-paper/NACE-97499 https://www.onepetro.org/conference-paper/NACE-97499 https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nbsspecialpublication674.pdf https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nbsspecialpublication674.pdf https://www.onepetro.org/conference-paper/NACE-07567 https://www.onepetro.org/conference-paper/NACE-07567 https://doi.org/10.1016/1353-2561(95)94483-Z https://doi.org/10.1016/1353-2561(95)94483-Z https://doi.org/10.1515/hjic-2018-0012 https://doi.org/10.1515/hjic-2018-0012 https://doi.org/10.1016/S1359-4311(02)00023-6 https://doi.org/10.1016/S1359-4311(02)00023-6 http://www.solublesolutions.com/wp-content/themes/twentyeleven_child/images/FoulChem2001.doc http://www.solublesolutions.com/wp-content/themes/twentyeleven_child/images/FoulChem2001.doc https://www.osti.gov/biblio/6358216 http://www.eptq.com/view_article.aspx?intAID=1197 http://www.eptq.com/view_article.aspx?intAID=1197 https://doi.org/10.1590/0104-6632.20150321s00002755 https://doi.org/10.1520/D4807-05R15 https://doi.org/10.1021/ef9901342 https://www.digitalrefining.com/article/1000381 https://www.digitalrefining.com/article/1000381 https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nbsspecialpublication674.pdf https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nbsspecialpublication674.pdf https://doi.org/10.31031/PPS.2018.02.000539 https://doi.org/10.31031/PPS.2018.02.000539 https://www.digitalrefining.com/article/1001216 https://www.digitalrefining.com/article/1001216 https://doi.org/10.1016/j.aca.2006.04.011 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 17–23 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-16 IMMOBILIZATION OF LACCASE IN ALGINATE BEADS BÉLA VARGA*1, MÓNIKA MEICZINGER2, AND VIOLA SOMOGYI1 1Institute of Environmental Engineering, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY 2Research Institute of Biomolecular and Chemical Engineering, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY The elimination or degradation of micropollutants from wastewater is becoming ever more important nowadays. Using ox- idoreductase enzymes to treat different micropollutants seems a promising solution. However, the viability of the process is highly dependent on the availability and stability of the applied enzymes. In order to improve the stability and provide faster reaction rates, enzymes can be immobilized in various carriers. Properties such as simple production, easy reten- tion and biodegradable carrier material are advantageous, e.g. entrapping laccase in alginate beads. This paper shows the results of the preparation and characterization of immobilized laccase entrapped in calcium alginate beads. The tech- nique of adding a mixture of sodium alginate and laccase dropwise into calcium chloride has been applied, improved and standardized to produce laccase-containing beads of uniform size and activity. For the purpose of characterization, a widely used substrate, 2,2’-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) diammonium salt, was used to evaluate the performance of the laccase-containing alginate beads. In addition to the characterization of the laccase-containing algi- nate beads, the enzyme kinetic constants (KM = 26.43 µM, Vmax = 0.23 µM/min) were determined. The reduction in the activity during storage has been described by a decay constant (0.26 d−1) that provides information concerning the de- sign constraints of the process. Results will be used to test the method in terms of the removal of organic micropollutants in continuous systems. Keywords: laccase, immobilization, alginate, ABTS, micropollutants 1. Introduction Awareness concerning the spread of micropollutants in water bodies has been raised in many countries [1] due to their toxic effects on ecological systems [2]. Additionally, these pollutants can infiltrate into reservoirs of drinking water in several ways resulting in public health concerns [3]. Since conventional treatment systems were designed to easily remove biodegradable compounds [4], the effi- ciency of micropollutant degradation is insufficient [5]. The partial elimination of these substances is driven by physical adsorption by the sludge [6] and cometabolism during microbial growth [7]. Therefore, the development of advanced technologies to improve the efficiency of re- moval and mitigate ecological effects is needed [8]. Organic micropollutants could be converted by oxi- doreductase enzymes such as tyrosinase, peroxidase and laccase. Laccase has received attention since it does not require additional cofactors. It promotes the production of an organic radical which combines with different trans- formation products [9]. The toxicity of these compounds is a subject of ongoing research [10], although a reduc- tion in toxicity has been shown in many cases [11]. *Correspondence: vargabela@almos.uni-pannon.hu Immobilization can be advantageous compared to the free form of the enzyme [12]. It can enhance the stability of the enzymes under more extreme conditions, namely higher pH and temperature [13], and makes the separa- tion of the catalyst from the effluent easier. For example, separation by applying an external magnetic force is pos- sible when laccase is immobilized on magnetic particles [14]. By using macro-sized supports, e.g. alginate beads [15] – unlike nanoparticles – the carriers can be retained by an outlet sieve, thus using them in continuous packed- bed reactors is possible [16]. Alginates extracted from brown algae [17] are good candidates for water treatment because of their advanta- geous properties [18], low-cost production and low health risks [19]. Calcium alginate microspheres, successfully used to immobilize polygalacturonase for the purpose of cleaning edible products, have shown improved opera- tional stability as the optimum temperature increases and even modifies the pH dependence of the enzyme [20]. Calcium alginate beads are formed by reversible cross- linking bonds with bivalent calcium ions between chains of α-L-guluronic acid and β-D-mannuronic acid [19]. Their unique properties allow them to be used as an ad- sorbent for dye removal [21]. Transformation of the phe- https://doi.org/10.33927/hjic-2019-16 mailto: vargabela@almos.uni-pannon.hu 18 VARGA, MEICZINGER, AND SOMOGYI nolic compound with laccase entrapped in alginate was successfully tested for bisphenol A [22] and the decolour- ization of textile dyes [23]. The aim of this study was to investigate the immobi- lization of laccase in calcium alginate beads and deter- mine its kinetic behaviour by measuring the conversion of 2,2’-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) diammonium salt (ABTS). Since this substrate is widely used, the results may be easily compared with other lac- case immobilization methods. 2. Materials and methods 2.1 Materials Laccase from Trametes versicolor (catechol oxidase ac- tivity ≥ 0.5 U/mg) and ABTS were purchased from Sigma-Aldrich. Food grade sodium alginate was pur- chased from Dragonspice Naturwaren. Citric acid and sodium phosphate used for making McIlvaine buffer so- lution [24] were supplied by VWR International. Calcium chloride was purchased from REANAL. Reagents were of analytical reagent grade except for sodium alginate. 2.2 Measurement of enzyme activity Laccase activity was determined by measuring the change in concentration of the ABTS radical formed in the catalysed oxidation. Measurements were performed in a Macherey-Nagel Nanocolor UV/VIS spectropho- tometer at 420 nm in a citrate-phosphate buffer solution at pH = 5 and within the range of laboratory temperatures (20± 3 °C). The formation of the product was registered by using the extinction coefficient of the ABTS radical (ε = 36, 000 M−1cm−1 [25]) and calculated using an equation A = EVsample lεVtotal (1) based on the Beer-Lambert law [26], whereA denotes the activity (U/l),E is a measure of the change in absorbance at 420 nm over 1 minute (min−1), Vsample represents the amount of enzyme contained in the sample used (ml), l is the pathlength of the beam of light (cm), ε stands for the extinction coefficient of the ABTS radical (M−1cm−1), and Vtotal expresses the total volume used (ml). In the case of immobilized laccase, beads were placed into the tubes instead of the liquid samples. The specific activity was calculated by measuring the weight and num- ber of beads before the activity in a closed tube was mea- sured. Between measurements, the beads were stored in distilled water in order to minimize the loss of water. 2.3 Immobilization of laccase inside an algi- nate matrix The immobilization protocol was based on entrapment of the enzyme before the formation of insoluble algi- nate beads according to Daâssi et al. [23]. The sodium- alginate crystals were dissolved in 5 ml of 1 mg/ml lac- case stock solution to produce a gel with 2 w/v% alginate content. This was added dropwise from a height of 25 cm into a 2 m/m% solution of calcium chloride where beads formed, then the beads were left for 30 mins to solidify. After this procedure, the beads were extracted from the solution by filtration and washed with distilled water in order to eliminate the unbonded laccase. The beads were stored in distilled water until used for further investiga- tions. The efficiency of the immobilization was determined by measuring the activities of the laccase stock solution (Ainitial), the solution of calcium chloride after immobi- lization (ACaCl2 ) and the washing liquid (Awashing) used. Immobilization efficiency (η) was calculated using the following formula η(%) = 100 Ainitial −ACaCl2 −Awashing Ainitial . (2) 2.4 Characterization of the immobilized lac- case The weight was measured frequently in a closed vessel by using an analytical balance (Kern ABJ 220-4NM). The size of the drops was determined by measuring the di- ameter of the freshly prepared beads after filtration on a square membrane filter. In order to determine the wa- ter content, twenty pieces of freshly prepared beads were put into an oven at a temperature of 105 °C and dried until their weight became constant. 2.5 Effect of temperature on alginate beads In order to gather information regarding the effect of tem- perature on enzyme activity within the optimal temper- ature range of wastewater treatment plants [27], the re- action rate of the ABTS transformation was determined inside a thermostatic cabinet (20− 35 °C). 2.6 Storage stability of the beads The storage stability of the immobilized laccase was de- termined by storing the laccase-containing beads in ultra- pure water at 20 °C without mixing. Beads were removed every 24 hours from the vessel in order to determine the change in their activity. The exponential decay constant was calculated by fitting the following equation At = Ainitiale −kdt (3) to the measured activities, where Ainitial and At indicate the specific activities [U/bead] at the initial and given time, respectively, kd denotes the decay constant [d−1] and t stands for the time of sampling [h]. 2.7 Reuse of the alginate beads The laccase-containing beads were mixed with a 0.05 mg/ml ABTS solution in a phosphate-citrate buffer so- lution (pH = 5) and the change in the ABTS radical con- centration was determined spectrophotometrically. After Hungarian Journal of Industry and Chemistry IMMOBILIZATION OF LACCASE IN ALGINATE BEADS 19 Figure 1: Change of size during drying of the beads (A) beads freshly prepared, (B) beads after drying, grid size: 3 mm measuring the activity of the beads (after 5 mins), they were filtered and washed with 5 ml of distilled water which was considered to be one cycle. During the fil- tration, a picture was taken to document the change in the size and shape of the beads. After washing, the beads were put into a fresh volume of buffer solution that con- tained ABTS and the same cycle repeated six more times. 2.8 Determination of kinetic parameters with ABTS Ten beads were put into a pH = 5 citrate-phosphate buffer solution, then different amounts of ABTS stock solution were pipetted in order to obtain different ini- tial substrate concentrations ([S] = 23, 46, 91, 228, 456, and 911 µM). The change in absorbance was recorded at 420 nm to calculate the concentration of ABTS radi- cals (P). The Michaelis-Menten constant (KM) and max- imum rate of reaction (Vmax) were determined by fitting the Michaelis-Menten equation. d[P] dt = Vmax[S] KM + [S] (4) with non-linear regression. 3. Results and Analysis 3.1 Characterization of the alginate beads Using the technique described, it was possible to produce a uniform bead weight: 0.0152 g (SD = 0.0016 g) of high immobilization efficiency (98.2 %). The diameters of the beads were between 2 and 3 mm. The weight and density of the beads constantly changed during their handling due to the rapid rate of evaporation of water. Therefore, the number of beads was used to calculate the specific activ- ity instead of weight, in a similar way to Daâssi et al. [23]. A specific activity of 0.005 U/bead was achieved, which equates to 0.005 µM production of ABTS radicals in one minute. The water content of the beads was determined to be 98± 1 %. The size of the beads shrank twofold whilst being dried as is shown in Fig. 1. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0 200 400 600 800 1000 A B T S ra di ca l co nc en tr at io n (µ M ) Time (s) Figure 2: Change of ABTS radical concentration 3.2 Determination of kinetic parameters The change in the concentration of ABTS radicals as a function of time is presented in Fig. 2. A linear curve could be fitted with the results during the initial 2 min- utes indicating a first-order reaction. The rate of the re- action can be determined from the gradient of the curves. By plotting the initial rates of reaction against the con- centration of the substrate, a lower initial substrate con- centration shows a better fit (R2 = 0.96 − 0.99) when compared to the higher range (R2 = 0.95). As a result, common Michaelis-Menten kinetics were shown (Fig. 3). By fitting an equation to the data gathered from three independent experiments with nonlinear least squares regression, the kinetic constants can be deter- mined, namely KM = 26.43 µM, Vmax = 0.23 µM/min. 3.3 Effect of temperature on alginate beads In the range of 20 − 35 °C, the rate of reaction, thus the enzyme activity, increased rapidly as can be seen in Fig. 4. As expected, the highest activity was measured at 35 °C since, in the case of the free laccase from Trametes versicolor, the optimal temperature for maximum activity was approximately 40 °C [26]. Figure 3: Fitting of Michaelis-Menten curve to the mea- sure. Initial concentrations used for calculation: 5 µM, 9 µM, 23 µM, 46 µM, 91 µM. 47(2) pp. 17–23 (2019) 20 VARGA, MEICZINGER, AND SOMOGYI 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 20 25 30 35 A ct iv ity ( U /b ea d) Temperature (�C) Figure 4: Effect of temperature on the activity of the im- mobilized enzyme. 3.4 Storage stability of the beads During storage, a slow decreasing trend can be observed in the activity of the beads, namely a reduction of 33 % over 5 days. By fitting equation to the measured values, the decay coefficient can be determined as 0.26 d−1. The change in activity and the decay curve produced by the equation are presented in Fig. 5. 3.5 Reuse of the beads The activity of the beads also decreased when repeatedly used. As an example, the results of an experiment were presented in Fig. 6. While the activity of the beads de- creased, their weight increased gradually. The change in the size of the beads and accumulation of the green reac- tion products can be seen in Fig. 7. This figure also shows the deformation of the beads in between the experiments. 4. Discussion Using the method described in this study, uniform beads can be produced with regard to their size and weight by controlling the size of their droplets and the height from which they are dropped. Since the alginate pow- der was dissolved directly in the stock solution of lac- case and the thick layer solidified instantly when it came 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0 1 2 3 4 5 6 7 A ct iv ity ( U /b ea d) Storage time (d) Measured Calculated Figure 5: Change of the activity during storage of the im- mobilized laccase. Figure 6: Change of activity and weight of the beads dur- ing reuse into contact with the CaCl2 solution, the entrapment was instant and effective. As a result, a high immobilization efficiency was achieved (98.2 %). Otherwise, the weight of the beads changed quickly whilst being handled after their formation as the beads dried out easily. However, the change in their water content did not influence their activities. Therefore, in order to achieve reliable results for specific activities, the number of beads was calculated rather than their weight. Evaluation of the reaction kinetics resulted in a Michaelis-Menten curve in the case of immobilized lac- case, in a similar manner to its free form [28]. This in- dicates that the mechanism of the reaction was not influ- enced by the entrapment. Nonetheless, the activity of the beads is different from that of free laccase due to steric effects and limitations with regard to the diffusion of re- actants. In order to achieve higher degrees of activity, the temperature could be increased. This would have ad- vantageous effects on the rate of reaction and diffusion. Additionally, by increasing the number of beads used, the rate of reaction could be fine-tuned to suit the re- quirements of their application in a similar way to the use of free enzymes. In order to determine the product yield or space-time yield, further investigations are rec- ommended. Due to the nature of ABTS, from visual inspection of radicals in itself, it was evident that the reaction products accumulate inside the beads (see Fig. 7). This could result in the deactivation of the enzyme by product inhibition. On the other hand, this phenomenon could be used to re- move the transformation products after enzymatic con- version by means of adsorption. As Tanaka et al. [29] concluded, the diffusion properties depend on the molec- ular size of the substrate. Since laccase-catalysed reac- tions tend to result in coupled products of higher molecu- lar weight [9], its release from the beads could be limited by diffusion. In the case of low-molecular-weight com- pounds, e.g. ABTS, free diffusion can be assumed in the alginate matrix [30]. As the experiments concerning repeated usage have shown, the beads suffer from major structural changes Hungarian Journal of Industry and Chemistry IMMOBILIZATION OF LACCASE IN ALGINATE BEADS 21 Cycle 1 Cycle 2 Cycle 3 Cycle 4 Cycle 5 Cycle 6 Cycle 7 Figure 7: Photos of the laccase containing alginate beads during the repeated use in buffer solutions that contain high concentrations of monovalent cations. This is due to the fact that cross- linking with calcium ions is reversible. These structural changes may result in the enzyme leaking and the specific activity decreasing following repeated usage. The reduc- tion in enzyme activity during storage could also be a re- sult of leakage. However, experiments performed in dis- tilled water fitted well with the decay equation, therefore, this can be taken into consideration during the design of the process. Moreover, in the case of a complex solution, e.g. treated wastewater or even tap water, the presence of bivalent cations may limit the deformation of the beads. 5. Conclusions Uniformly sized alginate beads were produced in the lab- oratory using a widely available biodegradable polymer. The mechanism of immobilized laccase was similar to that of its free form with regard to the conversion of ABTS, although the effects of diffusion and adsorption should be taken into consideration. Due to the structure of the alginate beads, the use of monovalent ions that contain buffers for measurements results in major struc- tural changes. The beads produced were reused success- fully over 7 cycles, while a rapid reduction in activity and structural changes were observed. However, the beads can be stored in distilled water and the reduction in ac- tivity estimated by an exponential decay equation which facilitates the design of the process. Although further in- vestigations regarding the efficiency of enzymatic trans- formation are needed, immobilization of laccase in al- ginate beads is a promising technique for the enzymatic treatment of micropollutants in continuous systems. Symbols A Activity [U/L] ACaCl2 activity of calcium chloride solution after immobilization [U/L] Ainitial activity of the laccase stock solution used [U/L] At specific activity at a given time [U/bead] Awashing activity of the washing liquid [U/L] E change in absorbance at 420 nm over 1 minute [min−1] kd decay constant [d−1] KM Michaelis-Menten constant [µM] l path length of light beam [cm] P reaction product, ABTS radical concentration [µM] S initial substrate concentration [µM] t time [min.; d] Vmax maximum reaction rate [µM/min] Vsample amount of enzyme-containing sample used [ml] Vtotal total volume used [ml] ε extinction coefficient of ABTS radicals[M−1cm−1] η immobilization efficiency [%] Acknowledgement The financial support of Széchenyi 2020 under the projects GINOP-2.3.2-15-2016-00016 and TÉT_15_IN- 1-2016-0094 is greatly appreciated. REFERENCES [1] Ebele, A.J.; Abou-Elwafa Abdallah, M.; Har- rad, S.: Pharmaceuticals and personal care prod- ucts (PPCPs) in the freshwater aquatic environ- ment. Emerg. Contam., 2017 3(1), 1–16 DOI: 10.1016/j.emcon.2016.12.004 [2] Shao, Y.; Chen, Z.; Hollert, H.; Zhou, S.; Deutschmann, B.; Seiler, T.-B.: Toxicity of 10 organic micropollutants and their mix- ture: Implications for aquatic risk assessment. Sci. Total Environ., 2019 666, 1273–1282 DOI: 10.1016/j.scitotenv.2019.02.047 [3] Ahmad, J.; Naeem, S.; Ahmad, M.; Usman, A.R.A.; Al-Wabel, M.I.: A critical review on or- ganic micropollutants contamination in wastew- ater and removal through carbon nanotubes. J. Environ. Manage., 2019 246, 214–228 DOI: 10.1016/j.jenvman.2019.05.152 [4] Monteiro, S.C.; Boxall, A.B.A.: Occurrence and Fate of Human Pharmaceuticals in the Environ- ment. In: Reviews of Environmental Contamination 47(2) pp. 17–23 (2019) https://doi.org/10.1016/j.emcon.2016.12.004 https://doi.org/10.1016/j.emcon.2016.12.004 https://doi.org/10.1016/j.scitotenv.2019.02.047 https://doi.org/10.1016/j.scitotenv.2019.02.047 https://doi.org/10.1016/j.jenvman.2019.05.152 https://doi.org/10.1016/j.jenvman.2019.05.152 22 VARGA, MEICZINGER, AND SOMOGYI and Toxicology. (Springer New York: New York, NY, 2010), pp. 53–154 DOI: 10.1007/978-1-4419-1157- 5_2 [5] Grandclement, C.; Seyssiecq, I.; Piram, A.; Wong- Wah-Chung, P.; Vanot, G.; Tiliacos, N.; et al.: From the conventional biological wastewater treatment to hybrid processes, the evaluation of organic microp- ollutant removal: A review. Water Res., 2017 111, 297–317 DOI: 10.1016/j.watres.2017.01.005 [6] Berthod, L.; Roberts, G.; Sharpe, A.; Whit- ley, D.C.; Greenwood, R.; Mills, G.A.: Effect of sewage sludge type on the partitioning be- haviour of pharmaceuticals: A meta-analysis. Env- iron. Sci. Water Res. Technol. 2016 2, 154–163 DOI: 10.1039/c5ew00171d [7] Torresi, E.; Tang, K.; Deng, J.; Sund, C.; Smets, B.F.; Christensson, M.; et al.: Removal of mi- cropollutants during biological phosphorus re- moval: Impact of redox conditions in MBBR. Sci. Total Environ., 2019 663, 496–506 DOI: 10.1016/j.scitotenv.2019.01.283 [8] Bui, X.T.; Vo, T.P.T.; Ngo, H.H.; Guo, W.S.; Nguyen, T.T.: Multicriteria assessment of advanced treatment technologies for micropollutants removal at large-scale applications. Sci. Total Environ., 2016 563–564, 1050–1067 https://doi.org/10.1016/j. scitotenv.2016.04.191 [9] Catherine, H.; Penninckx, M.; Frédéric, D.: Prod- uct formation from phenolic compounds removal by laccases: A review. Environ. Technol. Innov., 2016 5, 250–266 DOI: 10.1016/j.eti.2016.04.001 [10] Becker, D.; Varela Della Giustina, S.; Rodriguez- Mozaz, S.; Schoevaart, R.; Barceló, D.; de Cazes, M.; et al.: Removal of antibiotics in wastewater by enzymatic treatment with fungal laccase – Degrada- tion of compounds does not always eliminate toxi- city. Bioresour. Technol., 2016 219, 500–509 DOI: 10.1016/j.biortech.2016.08.004 [11] Varga, B.; Somogyi, V.; Meiczinger, M.; Kováts, N.; Domokos, E.: Enzymatic treatment and subsequent toxicity of organic micropollutants using oxidore- ductases - A review. J. Clean. Prod., 2019 221, 306– 322 DOI: 10.1016/j.jclepro.2019.02.135 [12] Brandi, P.; D’Annibale, A.; Galli, C.; Gentili, P.; Pontes, A.S.N.: In search for practical advantages from the immobilisation of an enzyme: the case of laccase. J. Mol. Catal. B Enzym., 2006 41(1–2), 61– 69 DOI: 10.1016/j.molcatb.2006.04.012 [13] Bilal, M.; Asgher, M.: Dye decolorization and detoxification potential of Ca-alginate beads immo- bilized manganese peroxidase. BMC Biotechnol., 2015 15, 111 DOI: 10.1186/s12896-015-0227-8 [14] Rouhani, S.; Rostami, A.; Salimi, A.: Preparation and characterization of laccases immobilized on magnetic nanoparticles and their application as a re- cyclable nanobiocatalyst for the aerobic oxidation of alcohols in the presence of TEMPO. RSC Adv., 2016 6, 26709–26718 DOI: 10.1039/c6ra00103c [15] Hegedüs, I.; Nagy, E.: Comparison of the Structure and the Stability of Single Enzyme Nanoparticles. Hungarian J. Ind. Chem., 2009 37(2), 123–130. [16] Niladevi, K.N.; Prema, P.: Immobilization of lac- case from Streptomyces psammoticus and its appli- cation in phenol removal using packed bed reactor. World J. Microbiol. Biotechnol., 2008 24(7), 1215– 1222 DOI: 10.1007/s11274-007-9598-x [17] Gao, F.; Liu, X.; Chen, W.; Guo, W.; Chen, L.; Li, D.: Hydroxyl radical pretreatment for low- viscosity sodium alginate production from brown seaweed. Algal Res., 2018 34, 191–197 DOI: 10.1016/j.algal.2018.07.017 [18] Fazekas, B.; Pitás, V.; Thury, P.; Kárpáti, Á.: Ex- amination of the Use of Support Materials of Nat- ural Origin in Wastewater Treatment. Hung. J. Ind. Chem., 2009 37(1), 5–10 [19] Li, Q.; Li, Y.; Ma, X.; Du, Q.; Sui, K.; Wang, D.; et al.: Filtration and adsorption properties of porous calcium alginate membrane for methylene blue re- moval from water. Chem. Eng. J., 2017 316, 623– 630 DOI: 10.1016/j.cej.2017.01.098 [20] Deng, Z.; Wang, F.; Zhou, B.; Li, J.; Li, B.; Liang, H.: Immobilization of pectinases into cal- cium alginate microspheres for fruit juice appli- cation. Food Hydrocoll., 2019 89, 691–699 DOI: 10.1016/j.foodhyd.2018.11.031 [21] Merakchi, A.; Bettayeb, S.; Drouiche, N.; Adour, L.; Lounici, H.: Cross-linking and modification of sodium alginate biopolymer for dye removal in aqueous solution. Polym. Bull., 2018 76(7), 3535– 3554 DOI: 10.1007/s00289-018-2557-x [22] Olajuyigbe, F.M.; Adetuyi, O.Y.; Fatokun, C.O.: Characterization of free and immobilized laccase from Cyberlindnera fabianii and application in degradation of bisphenol A. Int. J. Biol. Macromol., 2019 125, 856–864 DOI: 10.1016/j.ijbiomac.2018.12.106 [23] Daâssi, D.; Rodríguez-Couto, S.; Nasri, M.; Mechichi, T.: Biodegradation of textile dyes by im- mobilized laccase from Coriolopsis gallica into Ca- alginate beads. Int. Biodeterior. Biodegrad., 2014 90, 71–78 DOI: 10.1016/j.ibiod.2014.02.006 [24] MacIlvaine, T.C.: A buffer solution for colorimetric comparison. J. Biol. Chem., 1921 49, 183–186 [25] Childs, R.E.; Bardsley, W.G.: The steady-state kinetics of peroxidase with 2,2’-azino-di-(3- ethyl-benzthiazoline-6-sulphonic acid) as chro- mogen. Biochem. J., 1975 145, 93–103 DOI: 10.1042%2Fbj1450093 [26] Swinehart, D.F.: The Beer-Lambert Law. J. Chem. Educ., 1962 39(1), 333 DOI: 10.1021/ed039p333 [27] van Haandel, A.; van der Lubbe, J.: Handbook of Biological Wastewater Treatment. second ed. (IWA Publishing) [28] Gupta, N.; Lee, F.S.; Farinas, E.T.: Laboratory evolution of laccase for substrate specificity. J. Mol. Catal. B Enzym., 2010 62(3–4), 230–234 DOI: 10.1016/j.molcatb.2009.10.012 Hungarian Journal of Industry and Chemistry https://doi.org/10.1007/978-1-4419-1157-5_2 https://doi.org/10.1007/978-1-4419-1157-5_2 https://doi.org/10.1016/j.watres.2017.01.005 https://doi.org/10.1039/c5ew00171d https://doi.org/10.1039/c5ew00171d https://doi.org/10.1016/j.scitotenv.2019.01.283 https://doi.org/10.1016/j.scitotenv.2019.01.283 https://doi.org/10.1016/j.scitotenv.2016.04.191 https://doi.org/10.1016/j.scitotenv.2016.04.191 https://doi.org/10.1016/j.eti.2016.04.001 https://doi.org/10.1016/j.biortech.2016.08.004 https://doi.org/10.1016/j.biortech.2016.08.004 https://doi.org/10.1016/j.jclepro.2019.02.135 https://doi.org/10.1016/j.molcatb.2006.04.012 https://doi.org/10.1186/s12896-015-0227-8 https://doi.org/10.1039/c6ra00103c https://doi.org/10.1007/s11274-007-9598-x https://doi.org/10.1016/j.algal.2018.07.017 https://doi.org/10.1016/j.algal.2018.07.017 https://doi.org/10.1016/j.cej.2017.01.098 https://doi.org/10.1016/j.foodhyd.2018.11.031 https://doi.org/10.1016/j.foodhyd.2018.11.031 https://doi.org/10.1007/s00289-018-2557-x https://doi.org/10.1016/j.ijbiomac.2018.12.106 https://doi.org/10.1016/j.ibiod.2014.02.006 https://doi.org/10.1042%2Fbj1450093 https://doi.org/10.1042%2Fbj1450093 https://doi.org/10.1021/ed039p333 https://doi.org/10.1016/j.molcatb.2009.10.012 https://doi.org/10.1016/j.molcatb.2009.10.012 IMMOBILIZATION OF LACCASE IN ALGINATE BEADS 23 [29] Tanaka, H.; Matsumura, M.; Veliky, I.A.: Diffu- sion characteristics of substrates in Ca-alginate gel beads. Biotechnol. Bioeng., 1984 26(1), 53–58 DOI: 10.1002/bit.260260111 [30] Puguan, J.M.C.; Yu, X.; Kim, H.: Diffusion char- acteristics of different molecular weight solutes in Ca–alginate gel beads. Colloids Surfaces A Physicochem. Eng. Asp., 2015 469, 158–165 DOI: 10.1016/j.colsurfa.2015.01.027 47(2) pp. 17–23 (2019) https://doi.org/10.1002/bit.260260111 https://doi.org/10.1002/bit.260260111 https://doi.org/10.1016/j.colsurfa.2015.01.027 https://doi.org/10.1016/j.colsurfa.2015.01.027 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 25–30 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-17 MULTISTAGE GRAVITY BENEFICIATION OF RUTILE IN A TAR-FREE SAND RESIDUE RUKAYAT AKANDE*1 AND ABRAHAM ADELEKE1 1Department of Materials Science and Engineering, Obafemi Awolowo University, Ile-Ife, Nigeria This study reports on the concentration of rutile in the sand recovered from tar sand in Ondo State in Nigeria. The tar-free sand residue, approximately 90 % of which passes through a sieve with a pore size of 355 µm, was subjected to sieve analysis as well as sequences of panning gravity pre-concentration and shaking-table concentration at a slurry density of 25 % solids to improve the rutile content. The sand residue recovered in addition to the panned pre-concentrate and shak- ing table concentrates were also subjected to reflected light microscopy as well as transmitted light microscopy, counting using ImageJ software and X-ray fluorescence spectroscopy. The micrographs obtained showed that the samples contain rutile, dark-brown in color, interlocked with the major silica content and the content of rutile estimated by ImageJ software increased in the pre-concentrate from 7.90 % to 19.23 % in the final concentrate. X-ray fluorescence spectroscopy also showed that the rutile content increased in the pre-concentrate from 1.43 % to 31.02 % in the final concentrate. Therefore, the rutile content was successfully increased by the cheap gravity techniques of panning and shaking tables. Keywords: Tar sand, rutile, panning, shaking table, pre-concentrate, concentrate 1. Introduction Tar sands, also known as oil sands, are a combination of clay, sand and water saturated in a dense and extremely viscous form of petroleum technically referred to as bitu- men. Tar sands are impregnated sands that yield mixtures of liquid hydrocarbons, which require further processing other than mechanical blending before becoming finished petroleum products. Tar sand deposits are found in var- ious parts of the world including Canada, Madagascar, Venezuela, Russia, the United States and Nigeria [1]. Tar sand is exceedingly rich in oil as well as other valuable minerals and metals in varying proportions [2]. Tar sand deposits are composed primarily of quartz sand, silt, clay, water and bitumen along with trace amounts of metallic minerals such as rutile, pyrite, zircon and gemstones like tourmaline. Rutile is a major source of the element tita- nium and it is industrially used as a white pigment for paint, a ceramic glaze and in optical equipment [3, 4]. It is also used in sunscreen products due to its ability to re- flect ultraviolet light [5]. Nigeria has a large deposit of natural bituminous tar sand, which is estimated to be the fourth largest in the world after Canada, Russia and Venezuela [6]. It is esti- mated that about 34 to 45 billion barrels of heavy oil is trapped in tar sand deposits in Ondo State in Southwest- ern Nigeria alone, with more reserves in Edo and Ogun States. The reserve of tar sand in Ondo State alone is es- *Correspondence: akanderukayat2@gmail.com timated to be 31 billion metric tons [7]. The tar sand here is located in the eastern part of the Dahomey Basin, a coastal sedimentary basin that extends from the Ghana- Ivory Coast border through Togo and the Republic of Benin to Western Nigeria [8]. Adeleke et al. [9] leached tar sand from Ondo State using sodium hydroxide and sodium carbonate to strip off bitumen from the admixture of sand and also conducted a solubility test on the tar sand using toluene. Sodium hydroxide was found to yield a higher percent recovery of bitumen than sodium carbonate over the same leach contact time. Adebiyi et al.[10] determined the elemental composition of oil sand from Southwestern Nigeria by Total Reflection X-ray Fluorescence and detected twelve elements, namely K, Ca, Ti, V, Cr, Mn, Fe, Ni, Cu, Zn, As and Pb. Ignasiak [11] reported that a Nigerian deposit of bitumen-free sand residue contains only 0.06 wt% TiO2. Gravity beneficiation techniques make use of the dif- ference in the specific gravity of minerals to separate them. Therefore, heavy minerals such as rutile can be separated from the associated lighter minerals like silica. In the gravity method of shaking tables, a flowing film of water effectively separates coarse light particles from small dense particles [12]. Mixed rutile and ilmenite ores were subjected to coarse concentration using a mineral jig [13]. In this research, the lean content of the mineral rutile in a residue of tar sand from Ondo was increased over a sequence of gravity beneficiation by panning and shaking https://doi.org/10.33927/hjic-2019-17 mailto:akanderukayat2@gmail.com 26 AKANDE AND ADELEKE tables. 2. MATERIALS AND METHODS 2.1 The Collection of Samples The basic raw materials used in this work were lumps of bituminous tar sand obtained from Ondo State, Nige- ria. Lumps of the semi-solid black tar sand of about 50 kg in weight were collected from the villages of Agbabu and Loda in the Local Government Area of Irele in Ondo State. The Local Government Area of Irele is situated to the east of the Local Government Area of Okitipupa which lies at 40° 3” in longitude east of the Greenwich meridian and at 50° 45” and 80° 15” in latitude north of the equator. The samples were collected from two open pits at a depth of about 3 feet. After field sampling, the lumps of tar sand extracted were stored in large plastic bags. 2.2 Sample preparation The black lumps of tar sand were crushed and ground thoroughly with a steel rod. The ground tar sand was then homogenized in a mortar by pounding it with a pestle. The resultant mass of tar sand was composed of a finer mixture of particles and bitumen. 2.3 Stripping of bitumen by leaching whilst being stirred The tar sand was leached on a magnetic hotplate stirrer using a 32 factorial design, that is, at three different tem- peratures and by applying variables concerning the con- centration process. The 250 ml beaker reactor contain- ing the slurry of tar sand and 25 ml of 0.5 M sodium hydroxide at a pulp density of approximately 40 g/L was initially homogenized for 5 minutes. The beaker was placed on a Stuart magnetic hotplate stirrer, model num- ber R000101019, which was set to 50 °C and 90 rpm to yield the T1C1 test combination. The magnetic stir- rer was switched on and the pulp of tar sand allowed to react with the leachant, sodium hydroxide, for 30 min- utes. After 30 minutes, the concentrate of bitumen was skimmed off and the leachant filtered to remove the re- maining bitumen. The procedure was repeated for com- binations T1C2, T1C3, T2C1, T2C2, T2C3, T3C1, T3C2, and T3C3. 2.4 Stripping of bitumen by leaching in the ab- sence of stirring A slurry of tar sand with a pulp density of approximately 40 g/litre was prepared by mixing 25 ml of 1.5 M sodium hydroxide and 1 g of tar sand in a 500 ml beaker. The beaker containing the slurry of tar sand was then placed on a HP-11electric hot plate. The slurry was stirred con- tinuously for 30 minutes with a glass rod to agitate it to facilitate the recovery of bitumen from the tar sand. Sub- sequently, the bitumen was decanted off and the residue washed with water to remove the bitumen. The residue was washed again with toluene to dissolve traces of bitu- men in the tailings. The residue was then dried at 70 °C in an oven and the procedure repeated using 20, 3, 4, 5, 6, 7, 8, 9, and 10 g of tar sand which resulted in pulp den- sities of approximately 80, 120, 160, 200, 240, 280, 320, 360, and 400 g/L, respectively. The same procedure was repeated using 30 g of tar sand and the slurry prepared by mixing 300 ml of 1.5 M sodium hydroxide and 30 g of tar sand in a 500 ml beaker. The procedure was repeated several times until 7 kg of tar-free sand was recovered. A blank test was conducted using 25 ml of distilled water and 1 g of tar sand in a 500 ml beaker reactor. The beaker containing the slurry of tar sand was then placed on the HP-11 electric hot plate. The slurry was stirred continuously with a glass rod to agitate it to facilitate the recovery of bitumen from the tar sand. 2.5 Panning of the tar-free sand About 50 g of the sample was thoroughly mixed and ag- itated with 50 ml of water in a 500 ml plastic bowl. The panning process caused the lighter mineral particles to float while the denser particles sank. The floating light particles were decanted off while the denser residue was dried and weighed. The procedure was repeated but with 100, 150, 200, 250, 300, and 350 ml of water. Following preliminary panning, the remaining sam- ple that weighed 36.25 g was subjected to further pan- ning using the same 350 ml of water. The procedure was repeated ten more times by weighing the wet residue to- gether with a content of 4.02 g after each panning. The same procedure was employed for the shaking tabling of approximately 138 g of the sample. 966 ml of water was used for panning the 138 g of tar-free sand ten times and the final product weighed and dried before proceeding to the shaking table. 2.6 Determination of the specific gravity of tar-free sand A density bottle was washed, dried and its weight, m1, measured. Subsequently, 5 g of tar-free sand residue was placed inside the dried bottle and its weight, m2, deter- mined. Then, 50 ml of distilled water was poured into a measuring cylinder and weighed to give m4. Lastly, dis- tilled water was added to the tar-free sand in a covered bottle and shaken very well for thorough mixing, more- over, additional distilled water was added to increase the volume to 50 ml and then weighed to give m3. To deter- mine the specific gravity, thsese values were inserted into equation [12] G = m2 −m1 (m4 −m1)− (m3 −m2) (1) where m1 is the mass of empty bottle, m2 is the mass of bottle + dry soil, m3 is the mass of bottle + dry soil + Hungarian Journal of Industry and Chemistry MULTISTAGE GRAVITY BENEFICIATION OF RUTILE 27 Figure 1: Weight percentage of tar sand recovered follow- ing stripping on an electric hot plate. water, and m4 is the mass of bottle filled with water only. After deriving the specific gravity, the ratio of solids to water required to make the slurry was also determined as the feed should be composed of about 25 % solids by weight. 2.7 Shaking tabling A slurry composed of 25 % solids by weight was pro- duced by mixing 108 g of the panned sample in 324 ml of water which was introduced via the feed box of a shak- ing table, model number ED808148566Q. The slurry was distributed evenly over the table, wash water was dis- persed along the balance length of the feed launder. The table was then vibrated longitudinally using a slow for- ward stroke and a rapid return strike which caused the mineral particles to ‘crawl’ along the deck parallel to the direction of motion. As a result, the minerals were sub- jected to two forces, one due to the motion of the table and the other at right angles to it due to the flowing film of water. Consequently, the particles moved diagonally across the deck from the feed end and fanned out on the table – the smaller, denser particles rode on top towards the concentrate launder at the far end. The larger lighter par- ticles were washed into the tailings launder, which runs along the length of the table. An adjustable splitter at the concentrate launder was used to separate the product into three fractions, namely a high-grade concentrate and two medium-grade fractions. The concentrate continued to be shaken nine more times and the products weighed wet after every shaking. The final products were dried sepa- rately and then weighed [12]. 2.8 Particle size analysis using a sieving method The eight sieves with pore sizes of between 850 and 63 µm were chosen based on the square root of two rule. The weights of the empty sieves were recorded and stacked on top of each other with the coarsest sieve with a pore size of 850 µm on the top and the finest of 63 µm at the bottom. A fitting receiver was placed below the bottom sieve to collect any undersized particles and a lid placed on top of the coarsest sieve to prevent the sample from escaping. Fig. 1 shows the sieving machine and the setup of the sieves. A 300 g sample of the recovered tar-free sand was poured onto the uppermost, coarsest sieve, and the nest of sieves placed in an Endecotts test sieve shaker, model number 9205. The nest of sieves was vibrated in a vertical plane. The period of shaking was set to 30 mins using an automatic timer. During the shaking, the under- sized material dropped through successive sieves until it was retained on a sieve with a pore size slightly smaller than the diameter of the particles. In this way, the sam- ple was separated into size fractions. After the required time, the nest of sieves was dismantled and the amount of material retained on each sieve weighed [12]. 2.9 Reflected light microscopy The sample was deposited on a slide which was placed on the stage of an ACCU-SCOPE microscope, serial number 0524011. The diaphragm was then used to vary the inten- sity and size of the cone of light that was projected up- wards through the slide. The microscope was connected to a laptop to capture the plane section of the image which was captured using image capture software on a laptop as soon as a clear image was obtained. Reflected light mi- croscopy was used to view the raw sample, panned sam- ple, samples obtained from the particle size analysis and those from the shaking table. ImageJ was used to analyze the images obtained from reflected light microscopy. 2.10 Thin section microscopy The sample was impregnated since it was a loose sample by mixing equal amounts of araldite (resin and hardener) and sample together which was placed on a glass slide. The glass slide was placed on a hot plate for 5 mins and left to cool overnight. The cooled glass slide was ground by a lapping machine, transferred to a lapping plate for further grinding and viewed intermittently several times under the microscope till it became thin enough for light to penetrate through. The glass slide was dried on the hot plate, covered in araldite and viewed under a Brunel pet- rographic microscope, model no. 1287599. 2.11 X-ray Fluorescence spectroscopy A 3 g sample was pulverized into a fine homogenous mass and pelletized. An X-ray fluorescence spectrometer, model no. EDX3600B, was used to analyze the sample. The desired method was selected on the measuring instru- ment and the sample carefully placed onto the instrument according to the set-up of its benchtop measurement po- sition. The compartment containing the sample was cov- ered to prevent X-ray radiation from being scattered. The measurement conditions were set, the details of the sam- ple entered and the complete spectrum recorded which 47(2) pp. 25–30 (2019) 28 AKANDE AND ADELEKE Table 1: Screened distribution of sand recovered from tar sand. Pore-size range of sieve (µm) Weight of sieve fraction (g) Weight of sieve fraction (%) Nominal aperture size (µm) Cumulative Undersize Distribution (%) Cumulative Oversize Distribution (%) +850 1.02 0.37 850 99.63 0.37 850 to +500 0.92 0.33 500 99.26 0.74 500 to +355 5.05 1.84 355 97.45 2.55 355 to +212 20.44 7.45 212 89.97 10.03 212 to +150 178.7 65.14 150 24.86 75.14 150 to +125 50.51 18.41 125 6.45 93.55 125 to +90 16.62 6.06 90 0.39 99.61 90 to +63 0.96 0.35 63 0.04 99.96 < 63 0.11 0.04 lasted 60 to 120 seconds per sample, moreover, all de- tectable elements were measured simultaneously. Raw quantitative spectra and quantified results were acquired and saved. 3. RESULTS AND DISCUSSION It can be seen from Fig. 1 that the weight percent of recovered tar-free sand reduces as the mass of the tar sand increases which indicates that the percent recovery decreases as the pulp density increases. Therefore, pulp density is one of the factors that affects leaching, hence the laboratory-scale leaching of about 40 g/L was con- ducted. The decrease in the leaching rate as the pulp den- sity increases is due to increased particle crowding and a lower concentration of leachant available per unit vol- ume of the slurry [14]. It was observed that the smallest bitumen recovery of 55 % obtained from stirred leach- ing on the magnetic stirrer hotplate was greater than the highest recovery of 45 % for the unstirred leaching on the ordinary hotplate. This confirms that slurry stirring is an important process variable in the hydrometallurgical approach to leaching [14]. The sieve analysis showed that the majority of the sample was found within the 212 to 150 µm pore-size range, namely 65.14 % of the sample. The results ob- tained suggest that the residue of the tar sand is fairly coarse as the fraction smaller than 63 µm is insignificant. Coarse-sized ores are preferred for gravity concentration as they are more efficiently treated than finer ones [12]. The results are presented in Table 1. The lighter minerals from the 350 ml volume of wa- ter were further reduced following nine sequential pan- nings. The residue obtained after ten sequential pannings was found to contain many of the denser minerals given that the majority of the lighter minerals had been pre- viously removed as shown in Table 2. This table shows that the mass of the denser sand residue obtained after panning generally decreased as the volume of water used for panning increased. The results obtained suggest that panning efficiency increases as the volume of water in- creases. This may be due to the fact that by increasing Table 2: The recovery of denser materials as a function of the volume of water used during panning. Volume (ml) Residue (g) Mass recovered (%) 50 47.43 94.86 100 45.42 90.84 150 45.38 90.76 200 44.01 88.02 250 40.18 80.36 300 38.63 77.26 350 36.25 72.5 the water content, a more dilute slurry subject to less par- ticle crowding, a better degree of free flow and settling of mineral particles is produced, thus enhancing the separa- tion of lighter particles from denser ones. Fig. 2 indicates that generally the mass of denser residues decreased as the sequence of panning stages pro- gressed. The results suggest that as the sequence of pan- ning advances, the minerals composed of denser parti- cles are progressively concentrated into the concentrates of denser residues [12]. It was observed from Fig. 3 that the more the shak- ing table charge was subjected to a sequence of shaking tabling, the greater the proportion of denser rutile that was separated from the lighter silica mineral. The weight of the dry concentrates and tailings obtained after the se- Figure 2: Residue obtained after panning the sand residue with various volumes of water. Hungarian Journal of Industry and Chemistry MULTISTAGE GRAVITY BENEFICIATION OF RUTILE 29 Table 3: Summary of the numbers of total and rutile counts as well as the estimated percentage of the raw and panned samples in addition to the concentrates and tailings from the shaking table obtained from the XRF spectroscopy of each sample. Type of sample Total count Rutile count ImageJ estimated % of rutile The % of rutile by XRF spectroscopy Raw sand residue 1265 192 7.9 1.43 Panned sand residue 913 192 10.04 17.5 Shaking table concentrate of sand residue 1061 204 19.23 31.02 Shaking table tailings of sand residue 2041 99 4.85 5.01 quence of ten shaking tablings of the charge were 10.55 and 48.29 g, respectively. This figure shows that the mass of the concentrates of denser particles recovered gener- ally decreased while that of the tailings increased with some minor exceptions. The results obtained confirm the concentration of rutile in the concentrate increased as the sequence advanced. The sequential stages of shaking progressively improved the grade of the concentrates in terms of denser minerals. The plane-polarized and cross-polarized images ob- tained from reflected light microscopy and thin section microscopy of the raw and panned samples, in addition to the concentrate and tailings from the shaking table were analyzed. According to Francis [15], rutile was translu- cent and dark red-brown in color, while silica was trans- parent and whitish in color under a thin section micro- scope. Further observations showed that pure crystalline rutile was not abundant in the sample as most of the ru- tile under cross-polarized light was observed to be inter- locked with silica. The plates further indicated that the rutile content of the concentrate on the shaking table ex- ceeded that of the panned concentrate. A summary of the number of total and rutile counts as well as the percent distribution estimated for the raw and panned samples in addition to the concentrates and tailings of the shaking table is shown in Table 3. ImageJ counting estimated the rutile content in the recovered sand residue as well as in the concentrates of the shaking table and following pan- Figure 3: Recovered masses of concentrates and tailings from ten sequential shaking tables. ning as 7.9, 19.23 and 10.04 %, respectively. Fig. 4 shows the results from the X-ray fluorescence spectroscopy of the recovered and panned sand residues in addition to the concentrates from the shaking table. The percentage of silica in the recovered sand residue was reduced from 49.9 % to 47.6 % in the panned pre- concentrate and further to 41.4 % in the concentrate of the shaking table. Moreover, the titanium content of the recovered sand residue was enhanced from 0.86 % to 10.6 % in the panned pre-concentrate, while that of the concentrate from the shaking table in the panned pre- concentrate was further improved to 18.3 %. Therefore, the results obtained showed that the rutile content of the recovered sand residue was enhanced from 1.43 % to 17.50 % in the panned concentrate and finally to 31.02 % in the concentrate from the shaking table. The decrease in the silica content and increase in the titanium con- tent strongly suggest that the panning and shaking table were efficient methods to improve the rutile content of the residue of tar sand. The results obtained from ImageJ counting also compare favorably with that of XRF spec- troscopy except for the sample as received where both appeared to be significantly different. The flow diagram of the beneficiation process is shown in Fig. 5. Figure 4: Percent recoveries of the elements Ti, Si, Fe, Sn and Al in the concentrates. 47(2) pp. 25–30 (2019) 30 AKANDE AND ADELEKE Figure 5: Flow diagram of the Bench-scale Beneficiation of Tar Sand from Ondo. 4. Conclusions The rutile content of the tar-free residue obtained from tar sand extracted from Ondo was successfully improved by gravity beneficiation of the sand residue. The re- sults obtained show that a sequence of panning pre- concentrations followed by a sequence of shaking tabling improved the rutile content from 1.43 % to 31.02 %. REFERENCES [1] Falebita, O. A.; Koul, S.: Sustainable development of oil sands and host communities: preliminary sys- tem dynamics assessment. (Managing Intellectual Capital and Innovation for Sustainable and Inclu- sive Society: Proceedings of the MakeLearn and TIIM Joint International Conference 2, ToKnow- Press, 2015), pp. 2095–2109 [2] Allen, C.: Mineralogy of Tar Sand in the Upper Part of the Green River Formation in the Eastern Utah Basin, Utah. U.S. Geological Survey Open-file Re- port 76-381, pp. 27. [3] Erdogan, L. Dielectric Properties of Oil sands at 2.45 GHz Determined with Rectangu- lar Cavity Resonator. M.Sc. Thesis, Uni- versité de Montréal, Montréal, QC, Canada. https://publications.polymtl.ca/732/1/2011_leventErdogan.pdf Retrieved on 12 November 2019 [4] Lavat, A. E.; Gayo, G. X.: New Environmental Friendly Yellow Ceramic Pigments of the Type (Fe111MV)-Ti02. J. Chem. Chem. Eng. 2014 8: 1026–1035 DOI: 10.17265/1934-7375/2014.11.003 [5] Reddy, K. M.; Manorama, S. V.; Reddy, A. R.: Bandgap studies on anatase titanium dioxide nanoparticles. Mat. Chem. Phys. 2003 78(1): 239– 245 DOI: 10.1016/s0254-0584(02)00343-7 [6] Adegoke, O. S.; Ako, B. D.; Enu, E. I.: Geotech- nical Investigations of the Ondo State Bitumi- nous Sands, Geology and Reserves Estimate. Un- published manuscript. Department of Geology, Obafemi Awolowo University, Ile-Ife, Nigeria. [7] Emmanuel, E.; Ajibade, O. M.: Elemental Compo- sition and Geochemistry of Heavy Oil in Parts of Eastern Dahomey Basin, Southwestern Nigeria. J. Env. Earth Sci. 2014 4(12): 18–25 [8] Enu, E. I.: Textual Characteristics of the Nigeria Tar Sands. Sedimentary Geology, 44: 65–81. DOI: 10.1016/0037-0738(85)90032-6 [9] Adeleke, A. A.; Ibitoye, S. A.; Oladokun, C. B.; Oluwasegun, K. M.: An Evaluation of the Hot Aqueous Caustic Leaching of Nigerian Ondo Tar Sand. Petroleum and Coal, 2011 53(4): 320–323 [10] Adebiyi, F. M.; Asubiojo, O. I.; Ajayi, T. R.; Obia- junwa, E. I.: Trace element and physico-chemical characteristics of the sand and water fractions of Nigerian bituminous sands. Chemistry and Ecology, 2005 21(5): 369–380 DOI: 10.1080/02757540500291725 [11] Ignasiak, A. K.: Microscopic Structure of Athabasca Oil Sand. Can. J. Chem. Eng., 0982 60(4): 538–545 DOI: 10.1002/cjce.5450600416 [12] Wills, B. A., Napier-Munn, T.: Mineral processing technology (Elsevier, 2015) ISBN: 978-0-08-097053-0 [13] 911 Metallurgist. Metallurgists & Mineral Process- ing Engineers. https://www.911metallurgist.com Retrieved on 25 August 2018 [14] Ghosh, A.; Ray, H. S.: Principles of extractive metallurgy. (New Age International, 1991) ISBN: 9788122403220 [15] Francis, C.: Characterization and Upgrad- ing of a Rutile Concentrate Produced from the Oil Sands Tailings. M.Sc. Thesis, Univer- sity of Alberta, Edmonton, AB, Canada, pp. 95-96. https://era.library.ualberta.ca/items/79044641-a67d- 410f-b072-4bf5b4cadb0d/view/f051e5bd-3406-4262-abc8- dd479e66f11a/MQ81377.pdf Retrieved on 12 November 2019 Hungarian Journal of Industry and Chemistry https://publications.polymtl.ca/732/1/2011_leventErdogan.pdf https://doi.org/10.17265/1934-7375/2014.11.003 https://doi.org/10.1016{/}s0254-0584(02)00343-7 https://doi.org/10.1016/0037-0738(85)90032-6 https://doi.org/10.1016/0037-0738(85)90032-6 https://doi.org/10.1080/02757540500291725 https://doi.org/10.1002/cjce.5450600416 https://www.911metallurgist.com https://era.library.ualberta.ca/items/79044641-a67d-410f-b072-4bf5b4cadb0d/view/f051e5bd-3406-4262-abc8-dd479e66f11a/MQ81377.pdf https://era.library.ualberta.ca/items/79044641-a67d-410f-b072-4bf5b4cadb0d/view/f051e5bd-3406-4262-abc8-dd479e66f11a/MQ81377.pdf https://era.library.ualberta.ca/items/79044641-a67d-410f-b072-4bf5b4cadb0d/view/f051e5bd-3406-4262-abc8-dd479e66f11a/MQ81377.pdf HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 31–42 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-18 GPU-ACCELERATED SIMULATION OF A ROTARY VALVE BY THE DIS- CRETE ELEMENT METHOD BALÁZS FÜVESI1 AND ZSOLT ULBERT*1 1Department of Process Engineering, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY The rotary valve is the most frequently used piece of equipment that is suitable for the controlled feeding or discharging of products in powdered or granular form. It is usually connected to silos, hoppers, pneumatic conveying systems, bag filters or cyclones. In this paper, a simulation study is presented on the discharge of solid particles from a silo through a rotary valve. The discrete element method (DEM), which accounts for collisions between particles and particle-wall collisions, was used to model and simulate the motion of individual particles. The diameter of the simulated silo was 0.2 m and a total of 245, 000 particles were calculated. In the simulations, the effect of the geometric and operational parameters of the rotary valve on the mass outflow rate was investigated. The diameter of the rotary valve varied between 0.06 and 0.12 m and the rotational speed of the rotor was changed between 0.5 and 5 s−1. The simulations showed that the mass outflow rate of the particles from the rotary valve changes periodically due to its rotary cell structure. Within the lower range of rotational speeds of the rotor, the mass outflow rate of particles changes linearly in correlation with the rotational speed. The identification of this linear section is important in terms of control as this would facilitate the implementation of control devices by applying well-established linear control algorithms. Adjacent to the linear section, the dependence of the average mass outflow rate on the rotational speed was found to be nonlinear. Within the upper range of examined rotational speeds for each diameter of the rotary valve, the mass outflow rate reaches a maximum then decreases. The simulations were performed using GPU hardware. The application of parallel programming was an essential aspect of the simulations and significantly decreased the calculation time of simulations. In the treatment of particle-wall contacts, a novel flat triangular-based geometric representation technique was used which allows the particle-wall contacts to be calculated more effectively and their treatment implemented more easily into the parallel programming code. Using the calculated particle positions, the particles were visualized to view the effect of the interactions between the particles and rotor blades on particle motion. The simulation results showed that the discrete element method is capable of determining the detailed flow patterns of particles through the rotary valve at various rotational speeds. Keywords: silo, rotary valve, simulation, discrete element method, GPU 1. Introduction The handling of solids is part of many industrial tech- nologies in the chemical, pharmaceutical and food indus- tries amongst others. From simple storage (e.g. storage silo) and transport operations (e.g. pneumatic transport) to gas-solid or liquid-solid two-phase flows (e.g. gas- solid fluidization, gas-solid catalytic reactions), a number of technological operations involve solid and particulate materials. A common feature of gas-solid two-phase flow is that the interaction between the solid and gas phases can facilitate effective phase mixing in addition to inten- sive mass and heat transfer for a variety of operational purposes (e.g. gas-solid chemical reactions). The move- ment of the solid phase may result in different flow mech- anisms and flow patterns, depending on the frequency of collisions between the particles and effects as a result of the gas stream. The experimental investigation of the unit *Correspondence: ulbert@fmt.uni-pannon.hu operations of particulate solids is time-consuming and costly. Measurements in this type of systems are often in- accurate due to the limitations of the measuring devices and difficulty of taking measurements in a solid flow without changing its flow properties. Although a number of advanced measurement techniques for measuring the flow properties of solids, e.g. phase doppler particle an- alyzer, high-speed video processing, magnetic resonance tomography, positron emission particle tracking, etc., ex- ist which provide useful information on the movement of particles, many factors that significantly influence the flow structure such as the forces between the colliding particles are not yet available. However, experiments can be performed more easily using computer simulation if a sufficiently detailed model is available to investigate the system. The modeling and simulation of the flow of particulate solids are becom- ing increasingly popular and intensive research has been undertaken in this field. The flow of particulate solids is https://doi.org/10.33927/hjic-2019-18 mailto:ulbert@fmt.uni-pannon.hu 32 FÜVESI AND ULBERT always coupled with the flow of the interacting gas phase and in this sense the entire particulate flow can be con- sidered to be a gas-solid two-phase flow. However, the construction of the model depends on the degree of in- teraction between the solid and gas phases. The effect of the gas phase on the solid particles may vary from one unit operation to another. For example, in terms of the gravity discharge of solids from a silo, the effect of the gas phase is negligible on the motion of the particles and the calculation of any interaction between the solid par- ticles and gas phase and the calculation of the gas phase flow are unnecessary. On the other hand, when pneumatic transportation is used to transport solids, the interaction between the two phases is significant, therefore, the mod- elling and calculation of flow characteristics of both the solid particles and gas phase as well as their interactions is necessary. Two main methods for the treatment of solid particles with regard to gas-solid two-phase flows are proposed by the authors. The Euler method considers both the solid and fluid phases as a continuous phase and these meth- ods are referred to as two-fluid models (TFM) in the lit- erature. In these models, both phases are calculated by volume-averaged flow equations that are supplemented by terms which describe the interaction between the two phases. The disadvantage of these methods is that the solid particles are treated as a continuous phase and no information is provided about the motion of individual particles. In another method, the gas phase is treated as a continuous phase and described by volume-averaged flow equations, however, the solid particles are treated as dis- crete systems with their own mass, independent motion and velocity. The flow equations for the gas phase also in- clude terms that describe the interaction between the solid particles and the gas phase. To model and calculate the forces that result from particle-particle and particle-wall collisions, the discrete element method (DEM) based on the soft- or hard-collision model is used. These calcula- tion methods are referred to as the Euler-Lagrange, Euler- discrete element method (EULER-DEM) or the Compu- tational Fluid Dynamics Discrete Element Method (CFD- DEM). If the effect of the gas phase on the particles is negligible, the motion of the particles can be calculated by applying the discrete element method alone without taking into account any interactions with the gas phase. The discrete element methods have the advantage of pro- viding detailed information about the flow characteris- tics of solid particles. However, since all the particles are handled independently, they require significant comput- ing capabilities in the case of a large number of particles. In the application of “two-fluid” theory, Anderson and Jackson [1] first introduced interaction terms in the volume-averaged Navier-Stokes equation that is used to describe the particle-fluid and particle-particle interac- tions. In the 1990s, a significant amount of research using TFMs in the modeling and simulation of two-phase flow in fluidized beds was conducted by Bouillard et al. [2], Ding and Gidaspow [3], Kuipers et al. [4–6] and Nieuw- land et al. [7]. McKeen et al. [8] simulated the movement of catalyst particles in a fluidized-bed catalytic cracking reactor. Wang et al. [9] developed a TFM including chem- ical reaction source terms that are capable of describ- ing homogeneous and heterogeneous reactions. Pougatch et al. [10] simulated particle attrition in a fluidized bed introducing an attrition term into the model. However, given the complexity of two-fluid models, they are un- able to describe the discrete particle dynamics and flow properties of individual particles. The rapid development of computers and parallel computing over recent decades has led to the populariza- tion of discrete element methods in terms of the simula- tions of solid particulate flows and gas-solid two-phase flows. In the application of discrete element methods, the movement of each individual particle in the particulate flow is calculated. The motion of particles is calculated in accordance with Newton’s three laws of motion by taking into account the forces acting on the particles as a result of collisions with other particles or from their interactions with the gas phase. Over the last two decades, numerous papers have been published regarding the application of the discrete el- ement method. Tsuji et al. [11] and Kawaguchi et al. [12,13] presented simulations of gas-solid fluidization by applying a soft collision model developed by Cundall and Strack [14]. In these simulations, the mixing of particles inside the fluidized bed was examined. Hoomans et al. [15, 16] used a hard-sphere collision model in their sim- ulations and studied the effect of gas velocity on the flow of particles in a fluidized bed. Kaneko et al. [17] also used the CFD-DEM simulation method in the simulation of a polymer reactor. Using the temperature and velocity field of the gas phase, they calculated the kinetics of the poly- merization reaction and the reaction heat for each indi- vidual particle. Ulbert et al. [18] presented a CFD-DEM simulation study on the simulation of a reactive gas-solid two-phase flow in a circulating Mediator Recirculation Integrating Technology (MERIT) fluidized bed combus- tor. The model proposed by the authors contains the bal- ance equations for the components of the reaction and the shrinking core model was used to describe the gas- solid chemical reaction for each particle. Tsuji et al. [19] presented a simulation study that calculated the motion of 4.5 million particles in a fluidized bed using parallel computation. Yang et al. [20] examined the effect of the surface energy of particles on the transitions between dif- ferent flow types in a fluidized bed. Fan et al. [21] pre- sented a simulation study on the simulation of cohesive particles. Accordingly, the cohesive surface interaction between the particles was taken into account when the interacting forces were calculated. Hou et al. [22] simu- lated a fluid bed equipped with heat transfer tubes using the CFD-DEM simulation method. The authors examined the heat transfer between the tubes and the fluidized bed. Hou et al. [23] studied the effect of the particle-gas in- teraction, particle-particle collisions, gravity and friction on the transfer of energy between particles in particulate Hungarian Journal of Industry and Chemistry GPU-ACCELERATED SIMULATION OF A ROTARY VALVE 33 Figure 1: Cross-sectional view of a silo equipped with a rotary valve. solid flows and gas-solid two-phase flows. As can be seen from the aforementioned examples, in addition to the calculation of the motion of discrete particles, the discrete element method is suitable to treat the discrete nature of the solid particle phase by taking into account further characteristics of discrete particles, e.g. temperature, concentration, etc. The subject of this paper is the simulation of the op- erational characteristics of a rotary valve that is often connected to a silo that stores solids (Fig. 1). During the discharge of solids from silos, the effect of the gas phase on the movement of particles is negligible, i.e. it is unnecessary to consider the interaction between solid particles and the gas phase. Many examples are docu- mented in the relevant literature related to the simulation of discharging solids from silos using the discrete ele- ment method. Masson and Martinez [24] simulated silo feeding and discharge using the discrete element method. In their simulation studies, the effect of the mechanical parameters of the particles on particle flow was studied. It has been shown that the friction and elasticity of the parti- cles have a significant effect on the properties of particle flows. Yang and Hsiau [25] simulated silo feeding and discharge. In their research, different types of inserted el- ements were used and tested to improve the quality of the outflow of solid material from a silo. Goda and Ebert [26] presented a simulation study that investigated the effect of particles on a wall using the discrete element method. Gonzalez-Montellano et al. [27] examined a laboratory- sized silo experimentally and by simulation to identify the coefficient of friction of the particles. Zeng et al. [28] examined the fluctuation in the velocity of the particles being discharged from a silo. They showed that the ve- locity of particles fluctuates above a critical value of the coefficient of friction. Few reports in the literature are related to experimen- tal and simulation studies on rotary valves used to dis- charge or feed solids. A rotary valve connected to a silo was experimentally studied by Al-Din and Gunn [29]. The authors investigated the dependence of particulate mass flows delivered by the rotary valve on the rotational speed of the rotor and a linear relationship was found be- tween the rotational speed and mass flow only within a certain range of rotational speed. Kirkwood et al. [30] used a rotary valve to discharge solids from a silo and examined the mixing and segregation of solids inside the silo by analyzing the residence times. It was found that at low rotational speeds the particles segregated inside the silo and dead spaces developed in which the parti- cles moved down significantly more slowly. However, at higher rotational speeds the entire particle bed of the silo moved evenly downwards, thereby ensuring that the solid layers that were fed into the silo did not mix. The application of rotary valves occurs in all indus- trial areas where solids are stored in silos. The controlled discharge of the content of silos is important from both technological and control aspects. A deeper understand- ing of operational characteristics can significantly con- tribute to their geometric and operational design parame- ters. The simulation method used in our simulation study is the discrete element method, which provides a detailed insight into the flow of the solid phase inside a rotary valve. The purpose of the simulations is to determine how the geometric and operational design parameters of the rotary valve influence the mass outflow rate of solid ma- terial through the rotary valve. The following chapters of this paper describe the dis- crete element method used in the modeling and simula- tion of particle flow, the triangular geometric representa- tion technique used in the analysis of particle-wall colli- sions, the main features of the GPU hardware and parallel programming used in calculations, and the simulation re- sults of the investigation of the effect of geometric and operational parameters of rotary valves on the mass out- flow rate of solid particles. 2. Discrete element method The discrete element method has become a popular method in the modeling and simulation of particulate flows. It is capable of calculating the motion of individ- ual particles. The interactions between the particles and the boundary of the system are the collisions between particles and the particle-wall collisions which result in contact forces. The total contact force that acts on a sin- gle particle is the sum of the contact forces with regard to the forces that proceed from all the neighboring con- tacts. The individual particles exhibit two types of mo- tion: translational and rotational. Translational motion is 47(2) pp. 31–42 (2019) 34 FÜVESI AND ULBERT caused by the contact forces and gravity. Rotational mo- tion is generated only by the contact forces. In the case of when the effect of the gas phase is negligible, the three- dimensional translational and rotational motion of parti- cles is expressed by Newton’s three laws of motion as follows: miai = j=qi∑ j=1 (fcn,ij + fct,ij) +mig, (1) Iω̇i = j=qi∑ j=1 (Tij +Mij) , (2) where qi denotes the number of particles simultaneously in contact with particle i. mi, ai, and ω̇i are the mass, translational acceleration, and angular acceleration of particle i, respectively. fcn,ij and fct,ij stand for the nor- mal and tangential components of the contact force be- tween particles i and j, Tij is the torque at the point of contact between particles i and j due to the tangential contact force, Mij denotes the rolling resistance acting on particle i with regard to particle j, I represents the moment of inertia of particle i, g stands for the gravita- tional force vector, and ri is a vector from the center of particle i to the point of contact. The torque Tij is cal- culated using the tangential contact force by the vector cross product Tij = ri × fct,ij . (3) The rolling resistance Mij is expressed by the following equation (Wensrich and Katterfeld [31]): Mij = −µr rirj ri + rj |fcn,ij | (ωiri − ωjrj) , (4) where µris the rolling resistance coefficient and ri and rj denote the radii of particles i and j, respectively. When the particle collides with a wall, Eq. 4 is simplified as Mij = −µrri |fcn,ij | (ωiri) . (5) After the contact forces and torque have been calculated, the translational and angular accelerations can be ex- pressed from Eqs. 1 and 2 and the translational and an- gular velocities as well as position of particles can be cal- culated by numerical integration. The calculation of contact forces can be divided into two groups, namely the calculation of normal and tangen- tial contact forces. In our simulation, the contact forces were calculated in accordance with the soft-sphere con- tact model developed by Cundall and Strack [14]. In this model, the contact forces are calculated on the bases of simple mechanical models such as spring, dashpot and friction elements (Fig. 2 ). These mechanical elements influence particle motion through the parameters of the spring constant k, damping coefficient η and coefficient of friction µ. The sum of the normal forces is yielded by the normal forces generated by the spring and dashpot element: fcn,ij = (−knδn,ij − ηnvr,ij · nij)nij . (6) Figure 2: Soft-sphere contact model: (a) normal direction, (b) tangential direction. Similarly, in the tangential direction, the sum of tangen- tial forces acting on particle i is fct,ij = −ktδt,ij − ηtvrs,ij , (7) where vector nij denotes a unit vector from the center of particle i to the center of particle j, vr,ij represents the relative translational velocity of the colliding particles, δn,ijstands for the normal particle displacement, δt,ij is the tangential particle displacement vector and vrs,ij de- notes the relative slip velocity of the particles. 2.1 Calculation of the normal and tangential displacement of particles During elastic collisions between particles, at first par- ticles are deformed and compressed, then expand and start to move in the opposite direction. This mechanism is modeled in the soft-sphere contact model by assuming the shape of particles is unchangeable. When two parti- cles collide, they overlap with each other and the distance between their centers becomes smaller than the sum of their radii. Based on this overlapping, a normal displace- ment δn,ij , as shown in Fig. 3, can be defined as δn,ij = (di + dj)/2− |pi − pj | , (8) where di and pi in addition to dj and pj are the diameters and position vectors of particles i and j, respectively. The calculation of contact forces in the tangential di- rection is also based on a displacement defined in the tangential direction. As the particles first make contact, the points of impact A and B are the same points of the space (Fig. 3). The location of two points is continuously changing due to the translational and rotational motion of particles and a distance (C) between the two points evolves which is equal to the absolute value of the tan- gential displacement vector δt,ij . The direction of δt,ij is perpendicular to the unit vector nij and extends from the actual point of contact. As the collision proceeds, the absolute value and direction of the displacement vector Hungarian Journal of Industry and Chemistry GPU-ACCELERATED SIMULATION OF A ROTARY VALVE 35 Figure 3: Illustration of normal and tangential particle dis- placement. change according to the instantaneous relative slip veloc- ity at the point of contact. Moreover, the directions of the unit vector nij and tangential displacement vector δt,ij also change in accordance with the translational move- ment of particles during the collision (Fig. 4). For this reason, the value of the tangential displacement vector is corrected by two steps. At first, the original value of the tangential displacement vector is transformed into the new position of particles defined by the unit vector nij . In the second step, the transformed tangential displacement vector is corrected by the time integration of the relative slip velocity. 2.2 Calculation of the slip velocity The tangential component of the relative translational ve- locity vrt,ij can be determined as vrt,ij = vr,ij − (vr,ij · nij)nij , (9) where vr,ij = vp,i − vp,j is the relative translational velocity. The circumferential velocity vc at the point of contact is calculated by the angular velocity ω of the par- ticles. For particle i, it is given as vc,i = (|ri|ωi)× nij . (10) By using the circumferential velocities, the relative slip velocity can be determined from vrs,ij = (vc,i − vc,j) + vrt,ij . (11) 2.3 Calculation of the slide of particles When two particles come into contact, the normal force fcn,ij presses them together. The sliding motion of sur- faces introduces a tangential force of friction that acts on the surface in the direction which opposes the motion. In our simulations, the Coulomb’s law of friction is ap- plied to calculate the slide of particles. The magnitude of the tangential contact force vector is verified against the value of µ |fcn,ji|, where µ represents the coefficient of friction between particles i and j. If the absolute value of Figure 4: Change in particle orientation during the colli- sion. fct,ij exceeds that of µ |fcn,ji|, the particles slide and the value of the tangential contact force is changed in accor- dance with fct,ij = µ |fcn,ij | tij , (12) where tij = fct,ij/|fct,ij | defines a unit vector. When the particles slide, the tangential contact force fct,ij is con- sidered to be exerted only by the spring. Therefore, a new value of tangential displacement vector δt,ij , which is identical to the new value of the tangential contact force, must be calculated: δt,ij = µ |fcn,ij | k δij |δij | . (13) 2.4 Calculation of the damping coefficient The damping coefficient η of a dashpot can be related to the coefficient of restitution e according to η = −2 √ mk ln e√ π2 + (ln e) 2 (14) developed by analytically solving the equation that de- scribes a damped oscillatory system containing a mass, spring and dashpot (Tsuji et al. [11]). 2.5 Treatment of particle-wall collisions Particle-wall contacts can be deduced from particle- particle contacts. On the external side of the wall the par- ticle comes into contact with, a mirror particle is placed with an angular velocity of zero (Fig. 5). The unit vector nij , from the center of particle i to the center of mirror particle j, is always perpendicular to the wall plane. In the case of a motionless wall, the translational velocity of a mirror particle is zero. Consequently, the relative velocity between the contacting and mirror particles is the same as the velocity of the contacting particle. When the particle comes into contact with a moving wall, the translational velocity of the mirror particle is identical to the velocity of the point of impact. Under these conditions, the con- tact forces for particle i can be calculated in a similar way to in the case of particle-particle contacts. 47(2) pp. 31–42 (2019) 36 FÜVESI AND ULBERT Figure 5: The placement of a mirror particle during a particle-wall collision. 3. Conditions of simulation Simulation of the three-dimensional motion of particles was conducted in a silo equipped with a rotary valve as shown in Fig. 1. The dimensions of the silo and rotary valve are given in Table 1. The number of spherical par- ticles of uniform size used in the simulation study was 245, 000. The particle properties used in the simulation are summarized in Table 2. The initial particle bed com- pletely filled the silo above the rotary valve as shown in Fig. 1. In the simulation study, the effect of the geometric de- sign and operational parameters of the rotary valve on the mass outflow rate was examined. In total 48 simulation runs were accomplished. The diameter of the rotary valve varied between 0.06 and 0.12 m. The rotational speed of the rotor was adjusted between 0.5 and 5 s−1. The pa- rameters of the simulation runs are summarized in Table 3. The specifications of the computer and GPU used in the simulation study are shown in Table 4. 4. Geometric representation of the silo and rotary valve During the discharge of the silo through the rotary valve, particles come into contact with each other and the walls. The walls of the silo are stationary, on the other hand, in the case of the rotary valve, the walls of the body are stationary and the walls of the rotating rotor move. The detection and treatment of particle-wall collisions require the calculation of the distance between the parti- cle and wall which is based on the normal vector of the Table 1: Geometric and operational parameters. Height of silo, m 0.49 Diameter of silo, m 0.2 Half-angle of cone 45◦ Height of rotary valve, m 0.14 Diameter of rotary valve, m 0.06, 0.072, 0.084, 0.096, 0.108, 0.12 Rotational speed, s−1 0.5, 1, 1.5, 2, 2.5, 3, 4, 5 Table 2: Particle properties used in the simulations. Diameter, m 0.004 Density, kg/m3 1190 Spring constant, normal direction, N/m 800 Spring constant, tangential direction, N/m 800 Coefficient of restitution, e 0.6 Friction coefficient, µ 0.43 Rolling resistance coefficient, µr 10−4 Table 3: Parameters of simulation runs. Number of particles 245, 000 Simulation time, s 4 Timestep, s 5× 10−5 Interval for data saving, s 0.002 Number of standing wall elements 500− 600 Number of moving wall elements 24 wall surface. In the case of flat walls, the normal vec- tor can be determined easily. However, in the case of curved surfaces, the calculation of the normal vector of the surface and contact point is much more computation- ally complex. The number of calculations can also be in- creased significantly when collisions with moving walls are taken into consideration. Therefore, in our simulation technique, the equipment walls are approached in accordance with a set of flat tri- angles which allows wall collisions to be detected and contact forces to be easily calculated. This has two advan- tages over the conventional treatment of particle-wall col- lisions using the mathematical definition of curved sur- faces. On the one hand, any curvature wall can be easily defined by a set of flat triangles using methods developed in the field of computer graphics. The surfaces and ge- ometries created using computer graphics software can be saved and imported into the simulation program along with the other simulation parameters. As a result, the ge- ometry of the equipment can be altered in the program without changing the programming code. In the case of the application of moving walls, the description of the movement of the triangles is also produced by computer graphics software and imported into the simulation pro- gram. Another advantage is that the GPU can work more Table 4: Details of the computer hardware. CPU Intel Core i5-3570 3.4GHz Memory size 16GB Motherboard Gigabyte GA-Z77-D3H GPU Nvidia GeForce GTX 1080 Ti Number of CUDA cores 3, 584 GPU memory size 11GB GPU gigaFLOPS (Single precision) 10, 609 Hungarian Journal of Industry and Chemistry GPU-ACCELERATED SIMULATION OF A ROTARY VALVE 37 Figure 6: The flat triangular approximation of the silo ge- ometry. efficiently if the threads need to conduct similar calcu- lations simultaneously. Thus, it is advantageous for the GPU to process the treatment of many triangles of the same structure as opposed to just a few but different struc- tures of equations for the surface. A further advantage of the flat triangular approximation is that its accuracy can be increased where required by using more triangles, e.g. the approximation of the silo geometry by a set of flat triangles is shown in Fig. 6. In the treatment of particle-wall collisions, it is neces- sary to determine the position of mirror particles. When using the flat triangular approximation, firstly it is neces- sary to check whether the given particle can collide with the given triangle. If the particle collides, the point of the triangle closest to the particle is determined using the method developed by Mark W. Jones [32]. Then the posi- tion of the particle is mirrored on the opposite side of the triangle in the direction given by the normal vector. 5. Application of the GPU To apply the discrete element method effectively, the computer hardware used for the calculations has to be capable of parallel computing the motion of individual particles. In our simulation study, GPU cores and paral- lel programming were used to calculate the particles in parallel. A GPU core was assigned to each particle. A core can only write data in the memory associated with the particle but can read all the data from the memory as- sociated with the other particles to analyze collisions be- tween particles. As a consequence, all the data that results from the calculation of a collision between two particles is stored in two memory places. This solution had to be applied to avoid problems related to the writing of dupli- cates. Additionally, because of the way data is read, as- sociated problems may arise when new data is calculated from the data in both the previous and current iterations. To overcome these problems, two data units were used, one always contained the data of the previous iteration and the other stored that of the current iteration. The number of particles simulated may exceed the number of cores in the GPU, therefore, particles are cal- culated in consecutive batches. After calculating a batch of particles assigned to the GPU cores, a new group of particles is added to the GPU cores and they are calcu- lated. This process is repeated until each particle has been calculated, then the next iteration is executed. It is also important to note that during the parallel ex- ecution all the calculation data is stored in the memory of the GPU and that the CPU only gives instructions to the GPU about which calculations to perform. Therefore, at the beginning of the simulation, all the initial data con- cerning particles and walls are copied into the memory of the GPU and the CPU is used to instruct the GPU to perform the calculations. Since the copying of data be- tween the memories of the GPU and computer (host) is slow, this would increase the run time of the program, so it is important that the copying of data between the two memories is limited. For this reason, two computational cycles were defined. During the internal cycle, only the GPU is working, it executes a defined number of itera- tions. At the end of the internal cycle, control is passed on to the external cycle and the data is copied to the host memory and saved as a file. It is also important to note that executing the code on the GPU is effective if all the GPU cores conduct the same calculation. Attention should be paid to the use of branch instructions as they can significantly impair per- formance. 6. Simulation results and discussion In the aforementioned simulations, the startup process of silo discharge and the effects of the diameter of the ro- tary valve and rotational speed of the rotor on the mass outflow rate of particles leaving the rotary valve were in- vestigated. 6.1 The startup process of silo discharge The startup process of silo discharge through the rotary valve is shown in a cross-sectional view in Figs. 7 and 8 at different rotational speeds. Fig. 7 shows the startup process at a rotational speed of 0.5 s−1. The positions of the particles were visualized six times per revolution. At this rotational speed, the particles have sufficient time to fill the entire volume of the cell. Particles entering the cell are forced to pass by the blades of the rotor. When a cell reaches the outlet, all the particles fall out of the cell. Fig. 8 shows the startup process of discharge at a rotational speed of 5 s−1. In this case, the blades of the rotor rotate so fast that the particles do not have sufficient time to fill the entire volume of the cell. Approximately only one third of the volume of the cell is filled with particles. The 47(2) pp. 31–42 (2019) 38 FÜVESI AND ULBERT Figure 7: The startup process of silo discharge at a rotational speed of 0.5 s−1. blades push the particles so strongly that they collide with the wall of the outlet channel and accumulate. As a result, not all the particles are capable of falling out of the cell and some remain in the rotary valve. 6.2 Analysis of the mass outflow rate through the rotary valve During silo discharge through a rotary valve, the mass outflow rate exhibits a pulsating, periodic change over time as a consequence of the rotating cellular structure of the rotary valve. Following a finite period of time after the rotary valve starts to operate, a periodic change in the outflow of particles develops of constant amplitude and frequency. Depending on the geometric and operational parameters, the amplitude and frequency of this periodic change may vary. To compare the effect of the various geometric and operational parameters on the mass outflow rate, it is nec- essary to calculate the time-averaged pulsating mass out- flow rate. Therefore, the simulation time was sufficient to obtain at least ten cycles of change of constant amplitude and frequency. The time-averaged mass outflow rate was calculated from the simulation data saved after certain in- tervals of time. First, the number of particles discharged from the silo during the time interval between data saves was deter- mined. Next, the resulting data was multiplied by the volume and density of the particles and divided by the time interval between data saves to obtain the periodic time variation of the mass outflow rate. Finally, the time- averaged mass outflow rate was obtained by averaging the values throughout the simulation. For example, Figs. 9 and 10 show the dynamic change in the mass outflow rate of particles at rotational speeds of 0.5 and 5 s−1 when the diameter of the rotary valve was 0.12 m. When the rotational speed was 5 s−1, only the initial part of the entire simulation observed a high- frequency periodic change in the mass outflow rate more clearly. By comparing the two figures, it can be seen that Figure 8: The startup process of silo discharge at a rotational speed of 5 s−1. Hungarian Journal of Industry and Chemistry GPU-ACCELERATED SIMULATION OF A ROTARY VALVE 39 Figure 9: The dynamic change in the mass outflow rate at a rotational speed of 0.5 s−1. as the rotational speed increases, a periodic change of higher frequency and smaller amplitude was observed. Fig. 11 shows the dependence of the time-averaged mass outflow rate on the rotational speed of the rotor when the diameter of the rotary valve was varied. The curves represent different diameters of the rotary valve. It can be observed that as the diameter increases, the av- erage mass outflow rate rises. In the first section of the curves, the mass outflow rate increases linearly as the ro- tational speed of the rotor rises. This section is significant in terms of its control. Linear industrial controllers are well established and suitable for the application of control theory for linear systems. The linear control theory ap- plies to systems whose output is proportional to their in- put. These types of systems are governed by simple linear differential equations, while nonlinear systems are often governed by nonlinear differential equations and mathe- matical techniques developed to handle them are much more complex and far less general. Due to the linear de- pendence of these sections, the mass outflow rate can be easily controlled by common industrial control systems Figure 10: The dynamic change in the mass outflow rate at a rotational speed of 5 s−1. Figure 11: The dependence of the rotational speed on the time-averaged mass outflow rate using different rotary valve diameters. based on linear control algorithms. In Fig. 11, it can also be observed that the length of linear sections depends on the diameter of the rotary valve. The linear section in the case of the smallest diameter of 0.06 m occurs between rotational speeds of 0.5 − 3 s−1, while in the case of the largest diameter of 0.12 m it occurs between rotational speeds of 0.5− 1.5 s−1. Since a desired mass outflow rate can be more easily controlled in a linear section with a smaller gradient, it is advisable to keep the diameter as small as possible. Following the initial linear section, the mass outflow rate reaches a maximum value, the larger the diameter, the greater the maximum value of the mass outflow rate. If a mass outflow rate in excess of the maximum value is required, this can only be achieved by increasing the di- ameter and not by raising the rotational speed. As the di- ameter of the rotary valve decreases, the rotational speed at which the maximum is observed increases, as shown in Fig. 11. After the maximum is achieved, a decreasing sec- tion follows where the time-averaged mass outflow rate begins to decrease. Since the location of the maximum shifts towards higher rotational speeds as the diameter of the rotary valve increases, the decreasing section also shifts. Between the rotational speeds of 4 and 5 s−1, the mass outflow rates for each diameter converge. The optimal value of the diameter of the rotary valve can be determined based on the length of the linear sec- tion. A given mass outflow rate can only be achieved with a sufficiently large diameter. However, the use of a larger diameter is not recommended as the gradient of the lin- ear range increases as the diameter increases and control of the system becomes less precise. By considering the range beyond the maximum, it is clearly unnecessary to operate a rotary valve of a given diameter at higher rota- 47(2) pp. 31–42 (2019) 40 FÜVESI AND ULBERT Figure 12: The flow patterns of particulate flows at different rotational speeds of the rotor. tional speeds to increase the level of performance, since the maximum mass outflow rate cannot be exceeded, in fact it may even decrease. 6.3 Flow patterns of the particles at different rotational speeds of the rotor As was shown in Figs. 9 and 10, the mass outflow rate exhibits a constant periodic change after the startup pro- cess of discharge. In order to understand the phenomena that occur inside the rotary valve with a greater degree of accuracy, the positions of the particles were visualized to obtain the instantaneous particle flow patterns visible at various rotational speeds of the rotor (Fig. 12). In all cases, a cross-sectional view of the rotary valve is presented when the rotor blades were in the same po- sition and the operation was already within the range of the constant periodic change. Fig. 12 shows the gradual decrease in the charge of cells as the rotational speed in- creases. At a rotational speed of 0.5 s−1, the cells still have sufficient time to fully recharge. Between rotational speeds of 1.0 to 1.5 s−1, an empty volume fraction re- mains in the cells after being charged. At these rotational speeds, this empty volume fraction is still small and the accelerated feed and discharge periods can compensate for the decreasing mass outflow rate. At a rotational speed of 2 s−1, the time-averaged mass outflow rate reaches its maximum. When the rotational speed was 2.5 s−1, the re- sulting empty volume fraction of the cells becomes even bigger but the particles rest on the lower blade. At a ro- tational speed of 3 s−1 a significant change can be ob- served in the operation. At this rotational speed, the par- ticles seem to float between the blades of the cell, i.e. the speed of the blades is equal to the rate at which the par- ticles drop. At a rotational speed of 4 s−1, the particles are already pushed forward by the top blades. By further increasing the rotational speed (5 s−1), the time the cells have to be recharged is so small that the mass outflow rate begins to decrease even more. 7. Conclusions Generally speaking, it can be said that by using the dis- crete element method in particulate flow simulations, it is possible to observe particle flows in detail and to bet- ter understand related phenomena. The simulation study presented in this paper focused on the simulation of the discharge of solid particles from a silo controlled by a rotary valve using the discrete element method. In the simulations, the mass outflow rate of particles was ana- lyzed along with different geometric and operational pa- rameters of the rotary valve. By applying the discrete el- ement method, it was possible to dynamically simulate the movement of individual solid particles through the ro- tary valve and examine the mass outflow rate of particles. By visualizing the individual particles at their positions, the effect of the interaction between the particles and the blades of the rotary valve on the flow patterns of particles at different rotational speeds was observed. The simulations showed that the mass outflow rate of the particles leaving the silo changes periodically due to the rotating cell structure of the rotary valve. The depen- dence of the time-averaged mass outflow rate on the ro- tational speed of the rotary valve was found to be nonlin- ear within the upper range of examined rotational speeds. However, for each diameter of the rotary valve, a range of rotational speeds in which the mass outflow rate changes Hungarian Journal of Industry and Chemistry GPU-ACCELERATED SIMULATION OF A ROTARY VALVE 41 linearly with the rotational speed was observed. Within these ranges of linear speed, the rotary valve can be used more advantageously to perform different control tasks. Following the linear range, for each diameter of the ro- tary valve, the mass outflow rate reaches a maximum then decreases. The simulations were performed using GPU hard- ware. The application of parallel programming is an es- sential part of the simulations and can significantly de- crease the calculation time of simulations from a few days to just a few hours. By considering the treatment of particle-wall contacts, a novel flat triangular-based geo- metric representation technique was proposed which also made it possible to calculate the particle-wall contacts in parallel by GPU. Symbols p position vector of a particle, m v velocity vector of a particle, m/s a acceleration vector of a particle, m/s2 g gravitational acceleration vector, m/s2 f force vector, N T torque vector, Nm I moment of inertia of a particle, kgm2 d diameter of a particle, m r radius of a particle, m e coefficient of restitution, dimensionless k spring constant, N/m m mass of a particle, kg n unit vector, dimensionless t unit vector, dimensionless Greek letters δ normal displacement, m δ tangential displacement vector, m η damping coefficient, kg/s µ coefficient of friction, dimensionless µr coefficient of rolling resistance, dimensionless ρ density of a particle, kg/m3 ω angular velocity vector, s−1 Subscripts c circumferential cn normal contact ct tangential contact n normal direction, normal component r relative rt relative tangential rs relative slip t tangential direction, tangential component Acknowledgement This research was supported by the State of Hun- gary within the framework of the project 20385- 3/2018/FEKUSTRAT. REFERENCES [1] Anderson, T.B., Jackson, R.: A fluid mechanical description of fluidized beds. I&EC Fundamentals 1967, 6(4), 527–539 DOI: 10.1021/i160024a007 [2] Bouillard, J.X., Lyczkowski, R.W., Gidaspow, D.: Porosity distribution in a fluidized bed with an im- mersed obstacle. AIChE J. 1989, 35(6), 908–922 DOI: 10.1002/aic.690350604 [3] Ding, J., Gidaspow, D.: A bubbling fluidization model using kinetic theory of granular flow. AIChE J. 1990, 36(4), 523–538 DOI: 10.1002/aic.690360404 [4] Kuipers, J.A.M., Prins, W., van Swaaij, W.P.M.: Theoretical and experimental bubble formation at a single orifice in a two-dimensional gas-fluidized bed. Chem. Eng. Sci. 1991, 46(11), 2881–2894 DOI: 10.1016/0009-2509(91)85157-S [5] Kuipers, J.A.M., van Duin, K.J., van Beckum, F.P.H., van Swaaij, W.P.M.: A numerical model of gas-fluidized beds. Chem. Eng. Sci. 1992, 47(8), 1913–1924 DOI: 10.1016/0009-2509(92)80309-Z [6] Kuipers, J.A.M., van Duin, K.J., van Beckum, F.P.H., van Swaaij, W.P.M.: Computer simulation of the hydrodynamics of a two-dimensional gas- fluidized bed. Comput. Chem. Eng. 1993, 17(8), 839–858 DOI: 10.1016/0098-1354(93)80067-W [7] Nieuwland, J.J., Veenendaal, M.L., Kuipers, J.A.M., van Swaaij, W.P.M.: Bubble formation at a single orifice gas-fluidized beds. Chem. Eng. Sci. 1996, 51(17), 4087–4102 DOI: 10.1016/0009- 2509(96)00260-6 [8] McKeen, T., Pugsley, T.: Simulation and experi- mental validation of a freely bubbling bed of FCC catalyst. Powder Technology 2003, 129(1–3), 139– 152 DOI: 10.1016/S0032-5910(02)00294-2 [9] Wang, X., Jin, B., Zhong, W.: Three-dimensional simulation of fluidized bed coal gasification. Chem. Eng. Proc.: Proc. Int. 2009, 48(2), 695–705 DOI: 10.1016/j.cep.2008.08.006 [10] Pougatch, K., Salcudean, M., McMillan, J.: Simu- lation of particle attrition by supersonic gas jets in fluidized beds. Chem. Eng. Sci. 2010, 65(16), 4829– 4843 DOI: 10.1016/j.ces.2010.05.031 [11] Tsuji, Y., Kawaguchi, T., Tanaka, T.: Discrete particle simulation of two-dimensional fluidized bed. Powder Technology 1993, 77, 79–87 DOI: 10.1016/0032-5910(93)85010-7 [12] Kawaguchi, T., Tanaka, T., Tsuji, Y.: Numerical simulation of two-dimensional fluidized beds using the discrete element method (comparison between the two- and three-dimensional models). Powder Technology 1998, 96, 129–138 DOI: 10.1016/S0032- 5910(97)03366-4 [13] Kawaguchi, T., Sakamoto, M., Tanaka, T., Tsuji, Y.: Quasi three-dimensional numerical simulation of spouted beds in cylinder. Powder Technology 2000, 109, 3–12 DOI: 10.1016/S0032-5910(99)00222-3 47(2) pp. 31–42 (2019) https://doi.org/10.1021/i160024a007 https://doi.org/10.1002/aic.690350604 https://doi.org/10.1002/aic.690360404 https://doi.org/10.1016/0009-2509(91)85157-S https://doi.org/10.1016/0009-2509(91)85157-S https://doi.org/10.1016/0009-2509(92)80309-Z https://doi.org/10.1016/0098-1354(93)80067-W https://doi.org/10.1016/0009-2509(96)00260-6 https://doi.org/10.1016/0009-2509(96)00260-6 https://doi.org/10.1016/S0032-5910(02)00294-2 https://doi.org/10.1016/j.cep.2008.08.006 https://doi.org/10.1016/j.cep.2008.08.006 https://doi.org/10.1016/j.ces.2010.05.031 https://doi.org/10.1016/0032-5910(93)85010-7 https://doi.org/10.1016/0032-5910(93)85010-7 https://doi.org/10.1016/S0032-5910(97)03366-4 https://doi.org/10.1016/S0032-5910(97)03366-4 https://doi.org/10.1016/S0032-5910(99)00222-3 42 FÜVESI AND ULBERT [14] Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. Géotechnique 1979, 29(1), 47–65 DOI: 10.1680/geot.1979.29.1.47 [15] Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., van Swaaij, W.P.M.: Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidized bed: A hard sphere approach. Chem. Eng. Sci. 1996, 51(1), 99–118 DOI: 10.1016/0009- 2509(95)00271-5 [16] Hoomans, B.P.B., Kuipers, J.A.M., van Swaaij, W.P.M.: Granular dynamic simulation of seg- regation phenomena in bubbling gas-fluidized beds. Powder Technology 2000, 109, 41–48 DOI: 10.1016/S0032-5910(99)00225-9 [17] Kaneko, Y., Shiojima, T., Horio, M.: DEM simula- tion of fluidized beds for gas-phase olefin polymer- ization. Chem. Eng. Sci. 1999, 54, 5809–5821 DOI: 10.1016/S0009-2509(99)00153-0 [18] Ulbert, Zs., Hatano, H., Ohya, H., Endoh, S.: Com- puter simulation of gas-solid chemical reaction and particle motion in MERIT (Mediator Recirculation Integration Technology) circulating fluidized bed combustion system. Proc. The 6th SCEJ Sympo- sium on Fluidization, The Society of Chemical En- gineers, Japan, pp. 189–196, 2000. [19] Tsuji, Y., Yabumoto, K., Tanaka, T.: Sponta- neous structures in three-dimensional bubbling gas- fluidized bed by parallel DEM–CFD coupling sim- ulation. Powder Technology 2008 184(2), 132-140, DOI: 10.1016/j.powtec.2007.11.042 [20] Yang, F., Thornton, C., Seville, J.: Effect of surface energy on the transition from fixed to bubbling gas- fluidized beds. Chem. Eng. Sci. 2013, 90, 119–129 DOI: 10.1016/j.ces.2012.12.034 [21] Fan, H., Mei, D., Tian, F., Cui, X., Zhang, M.: DEM simulation of different particle ejection mechanisms in a fluidized bed with and without cohesive inter- particle forces. Powder Technology 2016, 288, 228– 240 DOI: 10.1016/j.powtec.2015.11.023 [22] Hou, Q.F., Zhou, Z.Y., Yu, A.B.: Gas–solid flow and heat transfer in fluidized beds with tubes: Ef- fects of material properties and tube array set- tings. Powder Technology 2016, 296, 59–71 DOI: 10.1016/j.powtec.2015.03.028 [23] Hou, Q.F., Kuang, S.B., Yu, A.B.: A DEM-based approach for analyzing energy transitions in granu- lar and particle-fluid flows. Chem. Eng. Sci. 2017, 161, 67–79 DOI: 10.1016/j.ces.2016.12.017 [24] Masson, S., Martinez, J.: Effect of particle mechani- cal properties on silo flow and stresses from distinct element simulations. Powder Technology 2000, 109, 164–178 DOI: 10.1016/S0032-5910(99)00234-X [25] Yang, S.-C., Hsiau, S.-S.: The simulation and ex- perimental study of granular materials discharged from a silo with the placement of inserts. Powder Technology 2001, 120, 244–255 DOI: 10.1016/S0032- 5910(01)00277-7 [26] Goda, T.J., Ebert, F.: Three-dimensional dis- crete element simulations in hoppers and si- los. Powder Technology 2005, 158, 58–68 DOI: 10.1016/j.powtec.2005.04.019 [27] Gonzalez-Montellano, C., Ramírez, Á., Gallego, E., Ayuga, F.: Validation and experimental calibration of 3D discrete element models for simulation of the discharge flow in silos. Chem. Eng. Sci. 2011, 66, 5116–5126 DOI: 10.1016/j.ces.2011.07.009 [28] Zeng, Y., Jia, F., Zhang, Y., Meng, X., Han, Y., Wang, H.: DEM study to determine the relationship between particle velocity fluctuations and contact force disappearance. Powder Technology 2017, 313, 112–121 DOI: 10.1016/j.powtec.2017.03.022 [29] Al-Din, N., Gunn, D.J.: Metering of Solids by a Rotary Valve Feeder. Powder Technology 1983, 36, 25–31 DOI: 10.1016/0032-5910(83)80005-9 [30] Kirkwood, C., Roberts, A.W., Craig, D.A.: The antisegregation and blending characteristics of a massflow hopper and rotary valve. ICBMH ’98 proceedings: 6th Int. Conf. on Bulk Materials, Storage, Handling and Transportation, Institute of Engineers, Australia, 1998 3, pp. 469–476 ISBN:1858256967 [31] Wensrich, C.M., Katterfeld, A.: Rolling friction as a technique for modelling particle shape in DEM. Powder Technology 2012, 217, 409–417 DOI: 10.1016/j.powtec.2011.10.057 [32] Jones, M.W.: 3D Distance from a Point to a Triangle. Technical report CSR-5-95, Depart- ment of Computer Science, Swansea University, 1995, pp. 1-5 http://cs.swan.ac.uk/~csmark/ publications/1995_3D_distance.html Hungarian Journal of Industry and Chemistry https://doi.org/10.1680/geot.1979.29.1.47 https://doi.org/10.1016/0009-2509(95)00271-5 https://doi.org/10.1016/0009-2509(95)00271-5 https://doi.org/10.1016/S0032-5910(99)00225-9 https://doi.org/10.1016/S0032-5910(99)00225-9 https://doi.org/10.1016/S0009-2509(99)00153-0 https://doi.org/10.1016/S0009-2509(99)00153-0 https://doi.org/10.1016/j.powtec.2007.11.042 https://doi.org/10.1016/j.ces.2012.12.034 https://doi.org/10.1016/j.powtec.2015.11.023 https://doi.org/10.1016/j.powtec.2015.03.028 https://doi.org/10.1016/j.powtec.2015.03.028 https://doi.org/10.1016/j.ces.2016.12.017 https://doi.org/10.1016/S0032-5910(99)00234-X https://doi.org/10.1016/S0032-5910(01)00277-7 https://doi.org/10.1016/S0032-5910(01)00277-7 https://doi.org/10.1016/j.powtec.2005.04.019 https://doi.org/10.1016/j.powtec.2005.04.019 https://doi.org/10.1016/j.ces.2011.07.009 https://doi.org/10.1016/j.powtec.2017.03.022 https://doi.org/10.1016/0032-5910(83)80005-9 https://doi.org/10.1016/j.powtec.2011.10.057 https://doi.org/10.1016/j.powtec.2011.10.057 http://cs.swan.ac.uk/~csmark/publications/1995_3D_distance.html http://cs.swan.ac.uk/~csmark/publications/1995_3D_distance.html HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 43–51 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-19 RESIDENCE TIME DISTRIBUTION-BASED ANALYSIS OF AN INDUSTRIAL-SCALE BIOGAS FERMENTER ANDRÁS TANKOVICS1, DÁVID TAKÁCS1, JUDIT SZENDEFY2, BÉLA CSUKÁS1, AND MÓNIKA VARGA*1 1Institute of Methodology, University of Kaposvár, Guba Sándor u. 40, Kaposvár, 7400, HUNGARY 2Biogas Plant, Magyar Cukor Zrt., Pécsi u. 10, Kaposvár, 7400, HUNGARY Residence Time Distribution (RTD) measurement-based analysis of mixing conditions on an industrial-scale (13, 000 m3) anaerobic digester of pressed sugar-beet slices at Kaposvár Sugar Factory of Magyar Cukor Zrt. was studied. The lithium salt tracing technique was applied, while the quantity of the lithium chloride tracer and the sampling of the effluent were designed by a preliminarily studied simulation model of mixing. The lithium concentration at the outlet was analysed by Inductively coupled plasma–optical emission spectroscopy (ICP-OES). Taking into account the geometrical arrangement, the biogas flow produced and the cyclically changing recycle flow, various mixing models were generated with different compartmentalization and flow structures by applying the method of Programmable Process Structures. The simulation-based approximate identification of the mixing model was accomplished by a heuristic approach that took into consideration multiple structures with changing mixing flows. A model with an advantageously small number of compartments and parameters was sought which satisfies the measured RTD. The results suggest the intensive mixing of upper levels with a poorly mixed lower level, which contributes to the long tail in the RTD. The actual set-up supports a good horizontal distribution of the sugar-beet slices and the microbial biomass, while the limited degree of vertical mixing helps to avoid the elutriation of the useful microbiome. The suggested mixing model will be combined with the formerly elaborated model involving 9 bacterial groups. Keywords: anaerobic digestion, mixing of digester, Residence Time Distribution, lithium tracing, Programmable Process Structures 1. Introduction The mixing of anaerobic digesters is a critical issue be- cause it should support the uniform distribution of raw materials to be digested and the bacterial biomass along the cross-section of the unit. Moreover, the excessive sed- imentation of the solid (bacterial) phase at the bottom of the digester should be avoided. However, an unnec- essarily high degree of mixing may elutriate the bacterial biomass that decreases the effectiveness of transforma- tions and may cause surplus environmental load. The computational modelling of anaerobic digesters was developed from the modelling of wastewater treat- ment and degradation [1, 2] and from the ADM models designed by IWA [3, 4]. A comprehensive review from 2013 is available [5]. Knowledge about mixing in large biogas digesters is still in its infancy, so the objective of this study is to broaden this by determining the minimum retention time of substrates fed into anaerobic digesters and estimate the distribution time of substrates before being extracted from the investigated digester. *Correspondence: varga.monika@ke.hu Over the last century, compartmentalized models were successfully used given the lack of computational fluid dynamics (CFD) models available, e.g. for study- ing an interesting scale-up problem, where the partially mixed pilot performed better than the perfectly mixed laboratory unit [6]. Moreover, at that time, a review com- paring compartmentalization with CFD was published [7] which concluded that CFD is the most scale–independent method for mixing and scale-up studies, however, it is hampered by limitations such as high computational de- mands and the inadequacy of submodel consideration, e.g. biokinetics. In a recently published detailed review [8], the authors concluded that “Compartmental models allow multi-scale modelling with low computational time compared to a full coupled model (e.g. reactive numerical simula- tions). Thanks to these main characteristics, compart- mental models are able to model complex full-size in- dustrial systems. An effective compartment model could handle multiple, multiphysics phenomenological models (detailed kinetic reaction scheme, complex heat and mass transfer model, population balance, etc.) that could not be included in CFD analysis. Observed deviation between https://doi.org/10.33927/hjic-2019-19 mailto:varga.monika@ke.hu 44 TANKOVICS, TAKÁCS, SZENDEFY, CSUKÁS, AND VARGA fully detailed CFD model and compartment model show very small results deviation despite of very significant re- duction in calculation time (3 orders of magnitude).” An interesting new example of a hypothesis-driven compartmental model is given in [9]. Two full-scale digesters from a biogas plant (2, 000 m3 and 1, 500 m3) equipped with different mixing sys- tems and filled with different substrates were investigated by Kamarad et al. [10]. To characterize the substrate distribution, solutions of lithium hydroxide monohydrate were used in tracer tests at concentrations of 45− 50 mg Li+/kg TS in digesters. The tracer concentration in the effluent of the digester was measured. Although the data calculated by CFD methods were in very good agreement with the results, a full comparison was not made Kaposvár Sugar Factory of Magyar Cukor Zrt. de- veloped an internationally straightforward anaerobic fer- mentation technology to generate onsite used and surplus energy with a reduced amount of waste emitted by pro- ducing fuel gas of high methane content from the pressed sugar beet slices. In a former PhD thesis [11], a detailed simulation model was developed and validated for an ap- proximately perfectly mixed pilot unit that took into con- sideration 9 pseudogroups of bacteria by applying the available earlier version of Direct Computer Mapping- based modelling methodology. The scaling up of the model to an industrial scale [12] requires more detailed knowledge about the hydro- dynamic conditions of the appropriately compartmental- ized volume. The objective of this work was to measure and analyse the Residence Time Distribution of a large in- dustrial unit by a lithium tracer technique. The final goal of this analysis was to embed the formerly developed de- tailed digestion model into the compartmentalized mix- ing model to enhance the investigated fermenter. 2. Experimental work The industrial digester (see Fig. 1) was a cylindrical unit with a diameter of 25 m and height of 28 m. The effluent was removed from a volume of approximately 13, 000 m3 at a level of 20 m. Spatially uniform feeding and appropriate mixing was ensured by the recycle flow from a level 20 m in height from the base of the unit and fed into the annularly placed 6 subsequent bottom segments with a prescribed cyclic change of the flow. The pressed sugar beet slices were also fed into this recycle flow via a screw feeder. Mixing was enhanced significantly by the increasing upward flow of the generated biogas. Moreover, three small mechan- ical mixers ensured the raw materials and bacteria were uniformly distributed in the lower third of the unit. The slowly accumulating inorganic residue could be removed from the middle of the bottom of the digester by means of a slowly rotating agitator. Figure 1: Schematic diagram of the industrial-scale fer- menter of 13, 000 m3 in volume 2.1 Tracing and Measurement Technique For the tracing of the flow, 40 kg of lithium chloride p.a. was dissolved in 150 litres of tap water. Considering the contamination of sugar beet slices with many ions, lithium was chosen given the sensitivity of its measure- ments. A blank sample was tested and a calibration pre- pared by adding known amounts of lithium chloride to the blank solution of realistic composition. This amount was added to the sludge flow of fresh sugar beet slices at a rate of 43 m3/h on average. The duration of tracing was 150 seconds, while the main recycle flow was switched off during this time. Before and after the tracer inlet, the recycle flow was maintained at 200 m3/h. Differing from the geometrically arranged order of annular segments 1-2-3-4-5-6, the actually applied or- der of feeding was 1-2-(3)-6-4-5, while during the exper- iments the inlet to segment #3 was closed (because of a malfunctioning valve). Accordingly, the cyclic feeding sequence was changed to 1-2-6-4-5 over 360 seconds. The tracer was fed into segment 4, while the samples were taken from the effluent and discharged from seg- ment 5 at a height of 20 m. The draft samples of 1-2 litres in volume were filtered before being analysed. The concentration of Li+ was measured by the ICP- OES spectrometer at Kaposvár University. 2.2 Modelling and Simulation Methodology For the modelling and simulation of the flow structure of the fermenter, the methodology of Programmable Process Structures (PPS) was applied [9, 10]. In PPS (see Fig. 2), the locally programmable struc- ture of process models could be generated from two gen- eral meta-prototypes and from the standardized descrip- tion of the actual process network automatically, resulting in the dynamic structure of unified state and transition Hungarian Journal of Industry and Chemistry RESIDENCE TIME DISTRIBUTION-BASED ANALYSIS 45 Figure 2: Schematic diagram of the generation of the process model by the method of Programmable Process Structures elements. The prototype elements, which describe the functionalities of the model, could also be derived from the general meta-prototypes. The freely editable, actual prototype programs contain symbolic input, parameter and output variables as well as a locally executable pro- gram code. In various applications, many state and tran- sition elements can be modelled with the same or similar, reusable local programs (referred to as actual prototypes). The state and transition elements of the actual model can be parameterized and initialized concerning their case- specific prototypes. Its execution, namely connection- based communication between the state and transition elements of the programmed structure, is solved by the general-purpose kernel program of the method. During the simulation, the actual elements are ini- tialized by initial conditions and parameters, moreover, the output values are recalculated stepwise by taking into consideration input and parameter data according to the associated local program prototype. The identified input and output connections of the extensive/intensive prop- erties and signals make the combined execution of the balance-based and signal-based functionalities possible. 3. Results and Analysis 3.1 Experiments The quantity of the lithium chloride tracer (40 kg) and the sampling of the effluent during the measurement of the RTD were designed by some preliminarily studied, ap- proximate mixing models. The concentration of lithium chloride in the effluent was measured every 0.5−2 hours during the first period and daily or even less frequently over the following 23 days (when production slowed down by ending campaign). The data are summarized in Table 1 and in Fig. 3. 3.2 Possible Flow Structures and Parameters By taking into consideration the geometrical layout, the produced biogas flow and cyclically changing recycle flow, various mixing models with different compartmen- talization and flow structures were generated by applying the method of Programmable Process Structures. Table 1: Measured concentrations of lithium chloride in the effluent Sample ID Date Time, hours Li+, mg/l 0 Blind 0 <0.100 1 2018.12.18 0.5 <0.100 2 2018.12.18 1 0.226 3 2018.12.18 2 0.347 4 2018.12.18 4 0.311 5 2018.12.18 6 0.237 7 2018.12.18 8 0.811 8 2018.12.19 26 0.627 9 2018.12.19 38 0.538 10 2018.12.20 48 0.625 11 2018.12.22 96 0.536 12 2018.12.24 144 0.401 12B 2018.12.26 192 0.516 13 2018.12.28 240 0.423 14 2019.01.02 360 0.325 15 2019.01.04 408 0.361 16 2019.01.07 480 0.384 17 2019.01.10 552 0.323 The most detailed compartmentalization (D1) is shown in Fig. 4. Here the sludge zone of the process unit was divided into three vertical layers which were divided further into an inner cylinder, surrounded by an annulus that was divided into six parts. The mixing flows between the 21 (approximately per- fectly mixed) compartments are as follows: - Circflowmix: cyclically changing recycle flow, su- perposed by a bidirectional mixing flow between the an- nular slices above each other; - Gasflowmix: bidirectional mixing flow, induced by the upward biogas flow between the compartments above 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 5000 10000 15000 20000 25000 30000 35000 40000 C o n ce n tr a ti o n , m g /L Time, hour Measured tracer concentration in effluent Figure 3: Measured tracer concentration of Li+ in the ef- fluent of the industrial-scale fermenter (40 kg of lithium chloride dissolved in 150 litres of tap water was added to a sludge of sugar beet slices that entered segment #4 at a flow rate of 43 m3/h for 150 seconds). 47(2) pp. 43–51 (2019) 46 TANKOVICS, TAKÁCS, SZENDEFY, CSUKÁS, AND VARGA Figure 4: The most detailed schematic diagram of the fer- menter (Scheme D1). each other (for practical considerations, circflowmix and gasflowmix can be integrated into a single flowmix); - Permix: peripheral bidirectional mixing flow be- tween the subsequent annular slices of the same layer , horizontally; - Radmix: radial bidirectional mixing flow between the annular slices and the inner cylinder of the same layer, horizontally. The flow rate of the biogas increased as the height of the unit increased. Accordingly, gasflowmix was higher between layers 2 and 3 than between 1 and 2. It is worth mentioning that the outflow rate of the ef- fluent was considerably less than the inflow rate of the sludge composed of pressed sugar beet slices because a significant proportion of the raw material was converted into biogas that escaped through the upper gas dome. In Scheme D2 (Fig. 5 which was finally proven to be the best solution), the inner cylinder and the related ver- tical and horizontal mixing flows were removed because they slowed down the mixing in Scheme D1 by an unfea- sible degree. As a further simplification, the upper layer was not decomposed at all. It should be noted that an interaction between the changes in the structures and their parame- ters exists. As an example, Scheme D3 (Fig. 6) behaved in a similar way to D2 with the highest horizontal mixing flows in the upper layer. Planned future research will fo- cus on using less compartments which might also be ad- vantageous in terms of embedding the detailed dynamic model of anaerobic digestion into the flow structure. Ac- cordingly, Scheme D4 (Fig. 7) was also tested, where only the lower layer was decomposed. As a marginal solution, the case of Scheme D5 (Fig. 8) was also studied with only vertical decomposition of the process unit. Figure 5: Schematic diagram of the best performing fer- menter (Scheme D2). 3.3 Generation of the Simulation Models The process network of the various schemes can be de- fined by a general declarative program. A characteristic part of the definitions of Scheme D2 is as follows: states([c11],[sludge]). states([c12,[sludge]). ... states([c36,[sludge]). states([env,[sludge,feed]). transitions([c11],[flowmix,permix]). transitions([c12],[flowmix,permix]). ... transitions([c26],[flowmix,permix]). transitions([c31],[permix]). transitions([c32],[permix]). Figure 6: Simplified schematic diagram of the fermenter without decomposition of the upper layer (Scheme D3). Hungarian Journal of Industry and Chemistry RESIDENCE TIME DISTRIBUTION-BASED ANALYSIS 47 Figure 7: A more simplified schematic diagram of the fer- menter (Scheme D4). Figure 8: Trial in the absence of horizontal compartmen- talization (Scheme D5). ... transitions([c35],[permix, feeding\_recirc]). transitions([c36],[permix]). dcode(sludge,[li]). dcode(feed,[li]). where states and transitions are declared by the state and transition elements in the given compartments, respec- tively. Meanwhile, the lowest level components of the state elements are described by the dcode() predicates. The transition-based representation of the flow struc- ture is defined by the facts of predicate: trans(TransitionName,Compartment, InputComponents,InputSigns, OutputComponents,OutputSigns). For example, trans(flowmix,[c13], [n([c13],[sludge]), n([c23],[sludge])], [n([c13],[sludge]), n([c23],[sludge])],[],[]). generates bidirectional connections for mixing flows be- tween compartments [c13] and [c23] of the first and second layers, respectively. Similarly, trans(permix,[c32], [n([c32],[sludge]), n([c33],[sludge])], [n([c32],[sludge]), n([c33],[sludge])],[],[]). generates peripheral flows between compartments [c32] and [c33] in the third layer. The initial concentrations and parameters of the vari- ous state and transition elements were recorded in an MS Excel file from where they were transformed into a tex- tual form of declarative predicates. According to Fig. 2, the automatic generation of the Programmable Process Structures was conducted: • from the two general prototypes [9, 10], • from the textual description of the actual process network, and • from the textual declaration of initial values and pa- rameters. The automatically generated structural model of Scheme D2 is illustrated in Fig. 9. 3.4 Local Programming of the Generated Process Structure The functionalities of the flow-structure models can be represented by locally executable programs embedded in the prototype elements. The prototype elements can be prepared from the copies of the meta-prototypes, e.g. the local program for the calculation of flowmix is presented below: Figure 9: Programmable Process Structure of Scheme D2. 47(2) pp. 43–51 (2019) 48 TANKOVICS, TAKÁCS, SZENDEFY, CSUKÁS, AND VARGA { program ( ’ p e r m u t a t i o n ( InpConcs , [ d ( [ Coord1 , s ludge , l i ] , [ C1 ] , g_m3 ) , d ( [ Coord2 , s ludge , l i ] , [ C2 ] , g_m3 ) ] ) , p e r m u t a t i o n ( P a r a m e t e r s , [ d ( vgas , [ Vgas ] , m3_h ) , d ( v c i r c , [ V c i r c ] , m3_h ) ] ) , g ( dt , DT) , p l u s f l o w ( Coord1 , Vci rc , Vplus ) , DM1 i s ( Vgas+Vplus ) / 6 0 * ( C2−C1 )*DT, DM2 i s ( Vgas+Vplus ) / 6 0 * ( C1−C2 )*DT, OutComps = [ d ( [ Coord1 , s ludge , l i ] , [DM1] , g ) , d ( [ Coord2 , s ludge , l i ] , [DM2] , g ) ] , Ou tS igns = [ ] , Re po r t = [ ] . p l u s f l o w ( [ Coord1 ] , Vci rc , Vplus ) :− sub_atom ( Coord1 , 2 , 1 , 0 , Column ) , g ( s e l e c t e d , Column ) , ! . p l u s f l o w ( _ , _ , 0 ) . ’ ) } This program illustrates how the input, parame- ter and output data are represented by the unified d(Identifier,List_of_Values,Dimension) triplets. The input data for calculations originated from the re- spective compartments Coord1 and Coord2. The pa- rameters define the Vgas mixing flow associated with the gas flow and the Vplus mixing flow generated by the recycle upflow between the given vertical segments. The actually active segment was determined by the inte- ger value of Column in the g(selected,Column) global predicate (which was actualized by another local program, responsible for the cyclic switching of the recy- cle flow). The auxiliary clause plusflow() considers mixing flow Vcirc only if the recirculation is actually associ- ated with the given compartments, otherwise the surplus mixing is equal to zero. The OutComps output of the program forwarded the changes in the amount of tracing component DM1 and DM2 in the compartments Coord1 and Coord2, respec- tively. 3.5 Evaluation of the RTD Measurements The mass flow of the produced biogas was calculated from the measured total volumetric flow rate of the three parallel operating digesters, divided according to the in- dividually measured loads of sugar beet slices of the par- allel lines whilst taking into consideration the measured composition of the biogas. The slightly changing output 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1000 2000 3000 4000 5000 6000 C o n ce n tr a ti o n , p p m Time, min Characteristic part of Shemes D1 and D5 Measured values D1 scheme D5 scheme Figure 10: Initial period of Scheme D1 that was insuffi- ciently mixed and Scheme D5 that was overmixed. flow of the effluent was calculated from the changing load and biogas mass flows. Regarding the identification and validation of the mul- tiple structure, it must be emphasized that this task was underdetermined. Considering the multiple interactions between the structures and parameters, instead of a rigor- ous optimization procedure, a heuristic trial and error ap- proach was applied, controlled by the main features first and then the values of the normalized root mean square error (NRMSE) in the refinement. The evaluation was ef- fectively aided by monitoring the change in concentra- tions in each compartment. First, D1 (being insufficiently mixed) and D5 (being overmixed) were excluded as can be seen in Fig. 10. Af- terwards, D2-4 were studied and based on the calculated NRMSE values, D3 and D4 were stepwise rejected. Fi- nally, the parameters of D2 were refined. The simulated and measured data are illustrated in Figs. Figs. 11–14 in addition to the calculated NRMSE values. Fig. 10 (and many trials using other parameters) showed that even Schemes D3-5 were unable to model the evident initial peaks of the measurements. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 5000 10000 15000 20000 25000 30000 35000 C o n ce n tr a ti o n , p p m Time, min Simulated effluent concentration for scheme D3 Simulated Measured Figure 11: Trial using horizontal compartmentalization in two lower layers only (Scheme D3, NRMSE = 12.10 %). Hungarian Journal of Industry and Chemistry RESIDENCE TIME DISTRIBUTION-BASED ANALYSIS 49 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 5000 10000 15000 20000 25000 30000 35000 C o n ce n tr a ti o n , p p m Time, min Simulated effluent concentration for scheme D2 Simulated Measured Figure 12: Trial using horizontal compartmentalization in three layers (Scheme D2, NRMSE = 12.08 %). The horizontal decomposition in all layers (without an inner tube) is able to express the main features using many parameters (an example is shown in Fig. 11). De- composition Scheme D2 (compartment volume is 685.3 m3) with improved parameterization describes the mea- sured data better (see Fig. 12). The change in the concentration in the various com- partments (except for the traced bottom and output top ones) is shown in Fig. 13. This shows that because of the changing input bottom compartment of the recycle flow, the initial peak in concentration does not appear in each compartment. The long tail is partly a consequence of the great dif- ference between the input feed and output effluent flows as well as of the limited degree of mixing between the layers and annular zones (especially in the bottom layer). It is also worth noticing the constant difference in concen- tration in the compartments. Unfortunately, the continu- ation of measurements was limited by the closing end of the campaign. Moreover, were they to continue, the sensi- tivity of the concentration measurements would decrease far below 0.3 ppm. The computation time, depending on the number of 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 5000 10000 15000 20000 25000 30000 35000 C o n ce n tr a ti o n , p p m Time, min Effluent concentration for scheme D2 (modified parameters) Calculated Measured Figure 13: Scheme D2 with a feasible parameter set (NRMSE = 21 %). 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 5000 10000 15000 20000 25000 30000 35000 C o n ce n tr a ti o n , p p m Time, min Calculated concentrations in the compartments of scheme D2 sludge [c22] li sludge [c23] li sludge [c24] li sludge [c25] li sludge [c26] li sludge [c21] li sludge [c32] li sludge [c33] li sludge [c35] li sludge [c34] li sludge [c11] li sludge [c12] li sludge [c13] li sludge [c15] li sludge [c16] li sludge [c36] li sludge [c31] li Figure 14: Change in concentration in the various com- partments of the experiment shown in Fig. 12. compartments and the flow structure, was between 2 and 90 minutes but was typically less than 10 minutes. A large difference in the volume of the compartments leads to an increase in computational efforts, while an equidistant volume distribution accelerates computation. 4. Discussion The results prove that the fermenter is thoroughly mixed but by no means perfectly mixed. Accordingly, detailed modelling needs to take into consideration hydrodynam- ics, e.g. by using mixing flows between the compartments of the volume. It should be emphasized that an interaction between the fermentation model and mixing exists because mixing is considerably enhanced by the produced biogas. In due course, this coupling has to be taken into consideration in the final implementation of the model. It is a specific feature of the Residence Time Distribu- tion of anaerobic digesters that the output liquid effluent flow is much less than the input load of the fresh sugar beet sludge (because of the production of biogas). This contributes to the slowly decreasing long tail of the RTD. By considering the approximately identified parame- ters, the mixing conditions can be characterized by the following features: • In the bottom layer, no bidirectional, horizontal mix- ing flow occurs between the annular segments. The horizontal mixing flow between the annular seg- ments can also be neglected in the middle layer. However, in the intensively bubbling upper layer, significant mixing occurs characterized by a flow rate of 200 m3/h between the adjacent segments. • Regardless of the very limited horizontal mixing in the bottom and middle layers, the appropriate distri- bution of the fluid and solid phases is achieved by the frequent change in location of the recycle flow (and feed) between the bottom segments. 47(2) pp. 43–51 (2019) 50 TANKOVICS, TAKÁCS, SZENDEFY, CSUKÁS, AND VARGA • The estimated mixing flow rate of the vertical sludge, generated by the upward biogas flow, was 150 m3/h and 300 m3/h between the bottom & mid- dle and middle & upper layers, respectively. • The estimated vertical mixing flow rate in the chang- ing active sites, generated by the recycle flow, was 200 m3/h between the vertically connected segments (i.e. the mixing ratio was 0.5). 5. Conclusions The Kaposvár Sugar Factory of Magyar Cukor Zrt. has developed an internationally straightforward anaerobic fermentation technology to generate onsite used and sur- plus energy with decreased emissions of waste by pro- ducing fuel gas of high methane content from pressed sugar beet slices. The development of an improved com- puter model as well as the analysis-based development of the technology required more detailed knowledge of the hydrodynamic conditions found within the process unit. The objective of this work was to measure and analyse the RTD of the appropriately compartmentalized indus- trial unit of 13, 000 m3 in volume. Spatially uniform feeding and appropriate mixing was ensured by the recycle flow, fed into the annularly placed six subsequent bottom segments with cyclic changes in location. The pressed sugar beet slices were also fed into this recycle flow via a screw feeder. Mixing was enhanced significantly by the upward increasing upstream flow of generated biogas and three mechanical agitators. By con- sidering the geometrical arrangement, produced biogas flow and cyclically changing recycle flow, various mixing models were generated with different compartmentaliza- tion and flow structures via the method of Programmable Process Structures. A lithium salt-based tracer technique was applied, while the quantity of the lithium chloride tracer and the sampling of the effluent during the measurement of the Residence Time Distribution were designed by some preliminarily studied simulation models of mixing. The lithium chloride concentration at the outlet was initially measured every 2 hours and daily or less frequently dur- ing longer periods. The lithium chloride concentration was analyzed by ICP-OES. For the simulation-based approximate identification of the mixing model, a heuristic approach was used that took into consideration multiple structures with chang- ing mixing flows. A model of an advantageously smaller number of compartments and parameters was sought which satisfies the measured Residence Time Distribu- tion. The results suggest the intensive mixing of upper lev- els with a limitedly mixed lower part that contributes to the long tail of the RTD. The applied arrangement sup- ports the horizontal distribution of the sugar beet slices to be digested and the multiple groups of bacteria. More- over, limited vertical mixing helps to avoid the elutriation of the bacteria and undigested organic materials into the environment. The suggested mixing model will be combined with the formerly elaborated model of the anaerobic fermen- tation process that consists of 9 groups of bacteria in our future work. Abbreviations ADM Anaerobic Digestion Model CFD Computational Fluid Dynamics ICP-OES Inductively Coupled Plasma Optical Emission Spectroscopy IWA International Water Association NRMSE Normalized Root Mean Square Error PPS Programmable Process Structures RTD Residence Time Distribution TS Total Solid content Acknowledgement The authors are grateful for the financial support from the project EFOP-3.6.1-16-2016-00007. REFERENCES [1] Batstone, D.J.; Keller, J.; Newell, R.B.; New- land, M.: Modelling anaerobic degradation of com- plex wastewater. I: model development. Biore- sour. Technol., 2000, 75, 67–74, DOI: 10.1016/S0960- 8524(00)00018-3 [2] Batstone, D.J.; Keller, J.; Newell, R.B.; Newland, M.: Modelling anaerobic degradation of complex wastewater. II: parameter estimation and validation using slaughterhouse effluent. Bioresour. Technol., 2000, 75, 75–85, DOI: 10.1016/S0960-8524(00)00019-5 [3] Blumensaat, F.; Keller, J.: Modelling of two-stage anaerobic digestion using the IWA Anaerobic Di- gestion Model No. 1 (ADM1). Water Res., 2005, 39, 171–183, DOI: 10.1016/j.watres.2004.07.024 [4] Fezzani, B.; Cheikh, R.B.: Modelling of the mesophilic anaerobic co-digestion of olive mill wastewater with olive mill solid waste us- ing anaerobic digestion model No. 1 (ADM1). Bioresour. Technol., 2008, 99, 6565–6577, DOI: 10.1016/j.biortech.2007.11.035 [5] Lauwers, J.; Appels, L.; Thompson, I.P.; Degrčve, J.; Impe, J.F.V.; Dewil, R.: Mathematical modelling of anaerobic digestion of biomass and waste: Power and limitations. Prog. Energ. Combust., 2013, 39, 383–402, DOI: 10.1016/j.pecs.2013.03.003 [6] Csukás, B.; Balogh, S.: Combining Genetic Pro- gramming with Generic Simulation Models in Evo- lutionary Synthesis. Computers in Industry, 1998, 36, 181–197, DOI: 10.1016/S0166-3615(98)00071-2 Hungarian Journal of Industry and Chemistry https://doi.org/10.1016/S0960-8524(00)00018-3 https://doi.org/10.1016/S0960-8524(00)00018-3 https://doi.org/10.1016/S0960-8524(00)00019-5 https://doi.org/10.1016/j.watres.2004.07.024 https://doi.org/10.1016/j.biortech.2007.11.035 https://doi.org/10.1016/j.biortech.2007.11.035 https://doi.org/10.1016/j.pecs.2013.03.003 https://doi.org/10.1016/S0166-3615(98)00071-2 RESIDENCE TIME DISTRIBUTION-BASED ANALYSIS 51 [7] Guillard, F.; Tragardh, C.: Modeling of the Per- formance of Industrial Bioreactors with a Dy- namic Microenvironmental Approach: A Critical Review. Chem. Eng. Technol., 1999, 22, 187– 195, DOI: 10.1002/(SICI)1521-4125(199903)22:3<187::AID- CEAT187>3.0.CO;2-9 [8] Jourdan, N. et al.: Compartmental Modelling in chemical engineering: A critical review. Chem. Eng. Sci., 2019, 210, 115–196, DOI: doi.org/10.1016/j.ces.2019.115196 [9] Norregaard, A. et al.: Hypothesis-driven compart- ment model for stirred bioreactors utilizing com- putational fluid dynamics and multiple pH sen- sors, Chem. Eng. J., 2019, 356, 161–169, DOI: 10.1016/j.cej.2018.08.191 [10] Kamarad, L.; Pohn, S.; Bochmann, G.; Harasek, M.: Determination of mixing quality in biogas plant di- gesters using tracer tests and Computational Fluid Dynamics. Acta Univ. Agric. Silvic. Mendelianae Brun., 2013, 61(5), 1269–1278, DOI: 10.11118/ac- taun201361051269 [11] Varga, M.: Economic optimization of sustainable complex processes (model based optimization un- der uncertain cost parameters for industrial scale anaerobic fermentation of sugar beet slice). PhD thesis, Kaposvár University 2009, (in Hungarian) [12] Csukás, B.; Varga, M.: Préselt répaszelet anaerob fermentációján alapuló biogáz előállítás számítógépi modellezésének előzetes vizsgálata. 2008. Témajelentés. Folyamatinformatika Kutató – Fejlesztő Bt., Kaposvár (in Hungarian) 47(2) pp. 43–51 (2019) https://doi.org/10.1002/(SICI)1521-4125(199903)22:3<187::AID-CEAT187>3.0.CO;2-9 https://doi.org/10.1002/(SICI)1521-4125(199903)22:3<187::AID-CEAT187>3.0.CO;2-9 https://doi.org/doi.org/10.1016/j.ces.2019.115196 https://doi.org/doi.org/10.1016/j.ces.2019.115196 https://doi.org/10.1016/j.cej.2018.08.191 https://doi.org/10.1016/j.cej.2018.08.191 https://doi.org/10.11118/actaun201361051269 https://doi.org/10.11118/actaun201361051269 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 53–61 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-20 PROGRAMMABLE PROCESS STRUCTURE BASED ANALYSIS OF HY- DROGEN SUPPLY CHAINS MÓNIKA VARGA *1, CATHERINE AZZARO-PANTEL2, JOSÉ MANUEL FLORES-PEREZ2, AND BÉLA CSUKÁS1 1Institute of Methodology, University of Kaposvár, Guba Sándor u. 40, Kaposvár, 7400, HUNGARY 2Laboratoire de Génie Chimique, Université de Toulouse, 4 Allée Emile Monso, Toulouse, 31030, FRANCE The planning and operation of Hydrogen Supply Chains (HSC) often require easily extensible, generic dynamic simulation tools. In this paper, the non-conventional modeling and simulation methodology of Programmable Process Structures is applied for the description of these process systems. Programmable Process Structures of HSC models are generated from the two general functional meta-prototypes of the method and from the description of the studied HSC network. The actual program prototypes of production, transformation, transportation, utilization and intermediate storage are copied from the meta-prototypes and filled with the locally executable declarative program code for the various classes of elements. The actual state and transition elements are parameterized and initialized according to their case-specific prototypes. The execution of the programmed HSC structures is solved by the general purpose kernel program. The application of methodology, developed in other fields, is illustrated by a fictitious, simplified HSC example. Analysis of this example model illustrates that the coordination of the hectically changing energy production of renewable resources with its seasonally and tendentiously changing demands is in need of dynamic simulation-based planning. Keywords: Hydrogen Supply Chains, dynamic simulation, model generation, simulation-based plan- ning, Programmable Process Structures 1. Introduction According to the recently published special section of the AIChE Journal about “Hydrogen Deployment” [1], the importance of hydrogen energy is expected to grow in the future because stakeholders have just started to discover the potential advantages of it worldwide [2]. Hydro- gen Supply Chains (HSCs) contain geographically dis- tributed, heterogeneous elements, while many of them are time-specifically constrained by meteorological or sea- sonal characteristics. Having analyzed some comprehen- sive case studies [3–5], it is obvious that the zero emis- sions of hydrogen burning and the relatively easy storage of hydrogen is perspectivic. At present, a deadlock ex- ists because a lack of infrastructure inhibits production, while a lack of production impedes the development of infrastructure. Dynamic modeling and simulation-based design tools may play an important role in accelerating developments. Primary renewable resources of hydrogen production are solar radiation, wind power and hydroelectric power produced by environmental systems. The utilization of solar energy requires natural or cultivated land, which is a special finite resource that is also generated for use in other industries, e.g. food-producing agriculture. *Correspondence: varga.monika@ke.hu The direct utilization of solar energy is achieved by various photovoltaic processes with a limited degree of efficiency and capacity. Intensive research and develop- ments are in progress that make use of many new re- sources involving special organisms (e.g. microalgae [6]) or artificial biomimetic systems (e.g. artificial leaves [7]). The indirect utilization of solar energy is realized by the construction of power plants, which also compete for land suitable for food production. Wind farms and some smaller hydroelectric power stations produce electrical energy periodically (often ran- domly), so on-site conversion to electrical energy is nec- essary for their effective utilization. However, electrical grids do not tolerate such erratic changes and no good solutions for temporary electricity buffering and storage exist. Primary non-renewable resources are nuclear raw materials, natural gas and coal. Nuclear raw materials need to be transported to specifically located nuclear power stations. Natural gas is available from widespread pipelines. On the contrary, coal-based processing requires stations to be constructed in the vicinity of mining sites. Large hydroelectric and nuclear power plants produce easily transportable electrical energy that can contribute to buffered hydrogen production by on-site electrolysis. https://doi.org/10.33927/hjic-2019-20 mailto:varga.monika@ke.hu 54 VARGA, AZZARO-PANTEL, FLORES-PEREZ, AND CSUKÁS The most frequently applied process, Steam Methane Reforming (SMR), competes with other applications of natural gas but can produce hydrogen-containing gas any- where in the vicinity of pipelines. On the contrary, coal gasification has to be conducted near to the resource, even though it produces more expensive hydrogen-containing gas from a resource, which is less convenient for tra- ditional applications and petrochemical production, than natural gas. In addition, everything depends on short-term prices, however, thorough preparation for the future re- quires longer-term decisions. Intermediate energy-containing products and byprod- ucts exist in HSC networks, which originate from many distributed districts. Biomass energy plants as well as some agricultural byproducts and waste appear season- ally, however, municipal biological wastes are constantly emitted. Biomass is processed by gasification, while biological wastes produce biogas by anaerobic fermentation. The enrichment of the hydrogen content in the fermentation process is supported by some experimental results [8, 9]. Energy-consuming, direct hydrogen-producing re- sources are electrical energy resulting from various chains, and hydrogen-containing gases that originate from SMR and gasification processes. The two corre- sponding kinds of hydrogen production processes are electrolysis and the different methods of gas separa- tion for hydrogen enrichment, e.g. membrane separation, pressure swing adsorption, etc. The produced raw hydrogen has to be conditioned in production sites by liquefaction or compression. The re- sultant liquidized or compressed hydrogen product needs a production or storage site for liquid or gaseous hydro- gen, respectively. Hydrogen can be transported from production to con- sumption sites by tanker trucks, tube trailers, railway tankers, railway tube cars or pipelines. Pure liquid or gaseous hydrogen needs to be stored at the site of con- sumption. The utilization of hydrogen is connected to these stor- age sites. In addition to fueling cars and other vehicles, especially aircraft, hydrogen can be used to generate ther- mal energy as well as in various industrial processes, e.g. ammonia production, synthesis of organic products, metal processing, etc. In the long-term development of energy supply, hy- drogen can serve as a clean and safe energy carrier that contributes to the adoption of a low-carbon economy. Considering the slow but definite development towards a hydrogen economy, the optimal design, planning and operation of a Hydrogen Supply Chain have come to the fore of research over the last decade. Supply chain network design is a dynamic field of re- search, not only in the context of HSC but in general. A comprehensive set of various fields of applications was reviewed by Govindan [10]. Regarding the generally applied modeling methods for the planning and design of HSC, a review from 2012 lists the main approaches [11] such as mathemat- ical optimization as well as Geographic Information Sys- tem (GIS) based models and frameworks. In conclusion, mainly mathematical programming-based optimization methods are used to design and plan HSCs worldwide. In 2017, another paper [12] provided an overview once more about the actual situation on the basis of almost 100 references. Various mathematical optimiza- tion methods such as mixed-integer linear program- ming (MILP), mixed-integer nonlinear programming (MINLP), mathematical programming problem (MPP), mathematical optimization (MOP) as well as GIS-based approaches and transition models illustrate progress over the last five years. Regardless of some alternatively ap- plied system dynamics (SD) and agent-based modeling (ABM) solutions, mathematical programming-based op- timization methods are still the most widely used tools in planning and design. This preferred approach is based on Mixed-Integer Linear Programming (MILP) in which binary variables denote the selection of technologies and establishment of transportation links, while continuous ones denote mass and energy flows, capacities of the supply chain nodes, as well as the cost and environmental performance metrics. These decision variables are optimized subject to mass balance constraints, capacity limitations and calculations of objective functions [13, 14]. Life Cycle Analysis (LCA) is an appropriate method- ology to evaluate the environmental performance and po- tential impacts associated with a product system across its entire life cycle. Among various LCA studies that evaluated the environmental benefits of using hydrogen over conventional energy sources, the so-called well-to- wheels (WTW) studies - an LCA approach to evaluate the environmental advantages of alternative fuel vehicles (AFVs) over conventional engine vehicles across the en- tire automotive fuel pathway - are studied including hy- drogen supply chains [15]. Multiple spatial and temporal scales imply a big chal- lenge. Samsatli et al. [16] highlights the need for the com- bined consideration of various temporal scales. A simple example of this is that the representation of energy stor- age requires a short (e.g. hourly) timescale, so taking into consideration years-long planning may result in computa- tionally intensive models. Dynamic simulation-based ap- proaches (including the actually studied Programmable Process Structures) may contribute towards solutions to these problems. Another big challenge is the question of managing un- certainty throughout the design optimization procedure. Accounting for multiple criteria, e.g. costs, environ- mental impacts, safety, etc., is an obvious expectation of green supply chain design, while the incorporation of these issues generates another interesting task to be solved [17]. Many aspects of these challenges in design, deploy- ment and operation in the field of HSC have been sum- marized in a recently published book [18]. Hungarian Journal of Industry and Chemistry ANALYSIS OF HYDROGEN SUPPLY CHAINS 55 Dynamic simulation Automatic generation Meteo data, changing demands Evaluated results of HSC Two general meta-prototypes Local program of specific classes Description of HSC structure Initial values and parameters Actual HSC to be studied Figure 1: Scheme of HSC model generation by the method of Programmable Process Structures. 2. Applied modeling method: Pro- grammable Process Structures By considering the aforementioned overview of chal- lenges that computational model-supported HSC plan- ning and operation face, a simple trial with our non- conventional approach was conducted. Accordingly, the methodology of Programmable Process Structures (PPS) [19, 20] was applied for the simulation-based analysis of a fictitious but typical Hydrogen Supply Chain. In PPS (see Fig. 1), the locally programmable struc- ture of process models can be generated from the GraphML-based graphical implementation of one state and one transition defining meta-prototypes, as well as from the unified textual description of the studied HSC network automatically resulting in a dynamic structure of the standardized state and transition elements. The pro- grammable prototype elements that describe production, transformation, transportation, utilization and intermedi- ate storage in the model can also be derived from the two general meta-prototypes. The freely editable, actual pro- gram of the prototypes contain symbolic input, parameter and output variables, as well as a locally executable GNU Prolog [21] code. In the various cases, the actual state and transition elements can be modeled with the same or a similar reusable local program, while the number of pro- gram prototypes is considerably less than the number of actual elements. The state and transition elements of the Pro- grammable HSC Structures can be parameterized and ini- tialized concerning their case-specific prototypes. During the simulation, the actual elements are executed by the prototypes associated with them. This execution results in cyclically stepwise, connection-based communication of the state and transition elements of the programmed structure with the general-purpose kernel that works like an operational system. When the simulation is run, the actual elements start with the initial conditions and parameters, while the out- put values are recalculated stepwise by taking into con- sideration the knowledge of input and parameter data in accordance with the associated local program prototypes. Distinguished input and output connectors as well as connections exist for the conservation measures (exten- sive/intensive properties) and signals, respectively. This solution supports the combined execution of the balance- based and signal-based functionalities, i.e. the hybrid models of hybrid processes. 3. Modelling and analysis of Hydrogen Supply Chains by Programmable Pro- cess Structures 3.1 Description of the example HSC process network In the following sections, the following notations of Pro- log syntax shall be used: Anything* = list of Anything; [H|T] = list of any elements or functors with a head H and a list of tail T; [] = empty list; full stop = end of logical sentence; if (or :-) implication in the sense of Horn clauses; comma in the program body of clauses = and. The general declarations of multi-compartment (multi- level, multi-scale) state and transition elements are as fol- lows: states(Coordlist,Statelist). transitions(Coordlist,Transitionlist). Coordlist = Coord* Statelist = State* Transitionlist = Transition* The transition-based declaration of the structures is de- scribed by the individual transitions as follows: trans(Transition,InpConcs,OutComps,InpSigns, OutSigns). InpConcs; InpSigns; OutComps; OutSigns= Specified_State* Specified_State = n(Coordlist,State) where the n() functors of InpConcs, InpSigns, OutComps and OutSigns contain the coordinate and name of the state connected with the given input mea- sure, input signal, output measure or output signal, re- spectively. This solution allows the same state names in the various compartments, levels or scales to be used be- cause they are identified by the spatial coordinate and name together. Our simplified example process network contains 38 state and 39 transition elements altogether in four regions, i.e. coordinates. The interpretation and abbreviation of the elements is explained in Table 1. The declaration of the process network is illustrated by a characteristic part of the definitions as follows: states([0], [raw_renew,raw_gas,raw_nucl,h2heating,co2,o2]). states([1], [h2cons,demand,h2carfuel,eleng,hydrogen,h2prod ]), etc ... transitions([0],[]). transitions([1], [t_h2refuel1,t_photovolt1,t_nucl1, t_electrolysis1,t_compression1, t_tankertruck11,t_h2heat1]), 47(2) pp. 53–61 (2019) 56 VARGA, AZZARO-PANTEL, FLORES-PEREZ, AND CSUKÁS Table 1: List of involved state and transition elements. Type of element Name of element Short symbol Renewable resources raw_renew Cultivated land raw_land Natural gas raw_gas Raw materials for nuclear en- ergy raw_nucl Energy-producing biomass biomass Biological waste biowaste Electrical energy eleng State Hydrogen-containing gas h2contgas element Raw hydrogen produced hydrogen Compressed H2 produced h2prod Compressed H2 available for consumers as a fuel h2cons Fuel demand of cars demand H2 utilized as a motor fuel h2carfuel Heat energy h2heating Carbon dioxide co2 Oxygen o2 Cultivation of land t_cultivation Photovoltaic process t_photovolt Wind power station t_wind Hydroelectric power station t_hydro Nuclear power station t_nucl Steam methane reforming t_smr Transition Biomass gasification t_gasify element Anaerobic fermentation t_ferment Electrolysis t_electrolysis Gas separation t_gassep Compression t_compression Transportation (tanker truck) t\tankertruck Refueling cars (vehicles) t_h2refuel Heat production t_h2heat etc ... trans(t_h2refuel1,[n([1],h2cons)],[n([1],h2cons ), n([1],h2carfuel)],[n([1],demand)],[]). trans(t_h2refuel2,[n([2],h2cons)],[n([2],h2cons ),n([2],h2carfuel)],[n([2],demand)],[]). trans(t_h2refuel3,[n([3],h2cons)],[n([3],h2cons ),n([3],h2carfuel)],[n([3],demand)],[]). trans(t_h2refuel4,[n([4],h2cons)],[n([4],h2cons ),n([4],h2carfuel)],[n([4],demand)],[]). trans(t_tankertruck11,[n([1],h2prod)],[n([1], h2prod),n([1],h2cons)],[n([1],h2cons)],[]). trans(t_tankertruck22,[n([2],h2prod)],[n([2], h2prod),n([2],h2cons)],[n([2],h2cons)],[]). trans(t_tankertruck33,[n([3],h2prod)],[n([3], h2prod),n([3],h2cons)],[n([3],h2cons)],[]). trans(t_tankertruck31,[n([3],h2prod)],[n([3], h2prod),n([1],h2cons)],[n([1],h2cons)],[]). trans(t_tankertruck32,[n([3],h2prod)],[n([3], h2prod),n([2],h2cons)],[n([2],h2cons)],[]). trans(t_tankertruck44,[n([4],h2prod)],[n([4], h2prod),n([4],h2cons)],[n([4],h2cons)],[]). trans(t_tankertruck41,[n([4],h2prod)],[n([4], h2prod),n([1],h2cons)],[n([1],h2cons)],[]). trans(t_tankertruck42,[n([4],h2prod)],[n([4], h2prod),n([2],h2cons)],[n([2],h2cons)],[]). trans(t_compression1,[n([1],hydrogen)],[n([1], hydrogen),n([1],h2prod)],[],[]). trans(t_compression2,[n([2],hydrogen)],[n([2], hydrogen),n([2],h2prod)],[],[]). trans(t_compression3,[n([3],hydrogen)],[n([3], hydrogen),n([3],h2prod)],[],[]). trans(t_compression4,[n([4],hydrogen)],[n([4], hydrogen),n([4],h2prod)],[],[]). etc ... 3.2 Meta-prototypes-based generation and analysis of the example model Generation of the Programmable Process Structure orig- inates from the GraphML description of general meta- prototypes that belong to the methodology (see their il- lustration in the upper part of Fig. 2) and from the afore- mentioned description of the actual process network. The initial values and parameters of state elements as well as the parameters of transition elements can be added automatically from an Excel based text file. The gener- ated general-purpose algorithm results in a GraphML file that contains the two meta-prototypes and structure of the whole HSC model in addition to the initial values and pa- rameters. The example structure is illustrated in Fig. 2. The GraphML model can be further edited by various graph editors (yEd Graph Editor to be precise). The structure also contains the two empty meta-prototypes. The generated Programmable Process Structure de- scribes all of the state (storage of resources, intermedi- ates and products) and transition elements that represent the transformation and transportation in the various com- partments (spatial scales) individually, because they may have different parameters and connections. By default, the transportations belong to the sending compartment. This representation does not result in the multiplication of the program code because many elements are executed by the same meta-prototypes that contain the local pro- grams (see the next Section). The connections from the state to transition elements represent the causally right transfer of quantitative and qualitative amounts and signals necessary for the calcu- lation of the transitions. The connections from the tran- sition to state elements represent the increases and de- creases in the corresponding additive measures as well as the calculated signals. Accordingly, the model is pre- pared for the common but distinguished representation of the extensive/intensive properties as well as of the signals and rules. Consequently, the Programmable Process Structures make both the elements and structure visible during the simulation process. This supports the easy (optionally multiscale) temporal consideration of external data, e.g. related to meteorology, the climate and supply & demand. Moreover, the flux routes and influence routes can be seen and evaluated explicitly. Furthermore, this (apparently too abundant) knowl- edge representation makes the unified execution of quite Hungarian Journal of Industry and Chemistry ANALYSIS OF HYDROGEN SUPPLY CHAINS 57 Figure 2: Programmable Process Structure of the investigated HSC. different hybrid dynamical models possible by applying the same kernel program. Figure 3: Some meta-prototypes-based case-specific pro- totypes of HSC networks. 3.3 Meta-prototypes-based functional proto- types of the example model First, the programmable prototypes have to be pre- pared by making copies from the application-independent meta-prototypes. Afterwards, the declaration of the actual local program follows. The implementation of state and transition prototype elements of the HSC model needs creative collaboration between the modeller and field ex- perts. The preparation of the case-specific local programs is supported by the meta-knowledge embedded in the meta-prototypes, in accordance with the generated Pro- grammable Process Structure. The basic meta-prototypes and some actual programmed prototypes are illustrated in Fig. 3. The example prototypes in Fig. 3 are as follows: Transition prototypes: prot_compr: compression of hydrogen; prot_electrol: electrolysis; prot_nucl: nuclear reactor-based electrical energy supply; prot_photovolt: photovoltaic system-based electrical energy supply; prot_refueling: hydrogen utilization by the refueling of cars (planes); prot_truck: transportation of hydrogen by truck ; prot_wind: wind power system-based electrical energy supply; State prototypes: prot_demand: scheduled hydrogen demand; 47(2) pp. 53–61 (2019) 58 VARGA, AZZARO-PANTEL, FLORES-PEREZ, AND CSUKÁS prot_storage: hydrogen storage; prot_condstorage: hydrogen storage, associated with condition signaling capabilities. As an example, see the simplified local program of prot_photovolt that calculates the produced elec- trical energy as a function of the actual meteorological conditions. This textual code is to be embedded into the respective prototype element as follows: Program\prot\photovolt = ’permutation(Parameters, [d(area,[Area],m2),d(coeff,[Coeff],nd)\textbar \], time(\,T,\), abs\time(T,Y,M,D), meteo(Y,M,D,Radiation,\,\,\,\,\), g(dt,DT), Energy is Coeff*Area*Radiation*DT, Report = [d(energy,[Energy],kWh)], OutComps = [d(eleng,[Energy],kWh)], OutSigns = [],!.’ 3.4 Generation and simulation of the example HSC system Programming of the process structure is followed by its automatic interpretation into the executable dynamic databases of the given model. Now the general model generator of the kernel reads the GraphML file and gen- erates: • the facts that describe the actual state & transition elements and those that determine state → transition and transition → state connections in the “user” file; • the clauses that declare the program prototypes in the “expert” file. The resultant user and expert files contain a detailed de- scription of the process model, including the locally exe- cutable program prototypes as well as the actual elements and connections that contain the initial data and parame- ters. Having generated the model, the executable code has to be supplemented by the (optionally scale-specific) global data, e.g. g(dt,DT ), that also declares the main time step in the previous example code. Dynamic simulation is executed by the general ker- nel program that comprises the cyclically repeated pro- cessing of state elements, state → transition connections, transition elements and transition → state connections. For example, the execution of all transitions associ- ated with the prototype of prot_photovolt are exe- cuted by the the unification of them, e.g. a(y,t_photovolt2,[2],prot_photovolt,[c( antecedent,dl,[]),c(consequence,dl,[]),c( cond,dl,[]),c(param,dl,[d(area,[1000],m2),d (coeff,[0.20000000000000001],nd)])],[i(inp, dl,[]),i(conc,dl,[])],[o(comp,dl,[]),o(out, dl,[]),o(report,dl,[])],[],[],[]). is unified with the prototype v(y,Spatial,prot_photovolt,[c(antecedent,dl, Antecedent),c(consequence,dl,Consequence),c (cond,dl,Conditions),c(param,dl,Parameters) ],[i(inp,dl,InpSigns),i(conc,dl,InpConcs) ],[o(comp,dl,OutComps),o(out,dl,OutSigns),o (report,dl,Report)],Temporal,Possibilities, Evaluations) :- Program_prot_photovolt. where Program_prot_photovolt refers to the ex- ample code above. The variables in prototypes are described by al- phanumeric symbols that start with a capital letter. The variables and constants in the functors i(), c() and o() are described by lists of the triplets d(Identifier, Variable_list, Instructions). 3.5 Analysis of the simulated results Some characteristic parameters in these illustrative exam- ples were the following: coefficient of photovoltaic model: 0.2 photovoltaic areas: 1000, 2000, 3000, 800 m2 coefficient of wind-power units: 20 kW h3/km3/ unit number of units: 100, 200 power consumption of electrolysis: 52.5 kWh/kg power consumption of compression: 3.0 kWh/kg Some of the simulated results are illustrated with a spe- cific part of a simplified HSC model. In this case, a sea- sonally and randomly changing continuous increase in the demand of cars for hydrogen fueling over the next decade was assumed, as can be seen in Fig. 4. In terms of production, the available capacity of pho- tovoltaic power stations and wind farms in the studied region was investigated. These renewable resources pro- duce electrical energy depending on the actual meteoro- logical conditions (solar radiation and wind speed). How- ever, the available climate scenarios did not contain data concerning radiation and wind speed, so meteorological data from the past ten years was used to imitate realis- tic changes in radiation and wind speed. The calculated electrical energy production of the investigated system is illustrated in Fig. 5. 0 2000 4000 6000 8000 10000 12000 14000 16000 k g/ d a y Time, day Daily hydrogen fuel consumption Figure 4: Hypothetical data for the seasonally changing and increasing demand for hydrogen. Hungarian Journal of Industry and Chemistry ANALYSIS OF HYDROGEN SUPPLY CHAINS 59 0 500000 1000000 1500000 2000000 2500000 3000000 3500000 0 5000 10000 15000 20000 25000 30000 35000 40000 P h o to vo lt a ic , k W h /d Time, day Photovoltaic and wind based electric energy PV production, kWh/d Wind W in d , k W h /d Figure 5: Electrical energy production of the involved wind power and photovoltaic systems (to be converted into hydrogen by electrolysis). In the studied area, wind power can produce more electrical energy than photovoltaic systems, however, it fluctuates very erratically. The contribution of photo- voltaic systems is smoother with characteristic seasonal changes. Increasing demands for hydrogen production can be met by more renewable energy production or additional hydrogen production from the available electrical energy (e.g. originating from nuclear power stations) and ad- ditional hydrogen production (e.g. from Steam-Methane Reforming, biomass gasification, etc.). In the case of lower fuel demands for cars (over a short initial period), the solar energy-based renewable resources can produce the necessary energy for hydro- gen production. Furthermore, some extraordinary surplus amounts can be produced over randomly appearing short periods by wind farms, as can be seen in Figs. 5 and 6. The results underline that dynamic simulation models need to be applied to solve planning problems in the case of rapidly changing supply and demand. 4. Discussion As has been summarized in the Introduction, mainstream HSC planning and operation applies mathematical pro- gramming, where a simplified model is usually embed- ded in a superstructure-based optionally multi-objective optimization program. The suggested non-conventional approach might have some potential advantages. Having analyzed the charac- teristics of Hydrogen Supply Chains and their typical ele- ments, Programmable Process Structures may offer some useful options: • the methodology tolerates case-specific multiscale extensions (e.g. when a proportion of the elements is in need of a more detailed model and/or refined spatial or temporal resolution). • by taking into consideration the intensive devel- opment of new hydrogen-producing methods, Pro- -20000 -10000 0 10000 20000 30000 40000 50000 60000 k g H 2/ d a y Time, day Uncovered daily H2 demand Figure 6: The change in hydrogen demand produced by the available direct renewable resources. grammable Process Structures are capable of em- bedding any detailed prototypes of nonlinear dy- namic models into existing structures. • the extensibility of Programmable Process Struc- tures favours stepwise structural and functional changes in supply and demand without reformula- tion of the model. • optional multiscale compartmentalization helps to represent multiple districts and sites via the optional connection between spatial coordinates of the ele- ments with an external GIS tool. • individual elements and connections can be associ- ated with a time-driven execution at given points in time and/or over prescribed time intervals. • material and energy balances can be combined with signal and rule-based, event-driven processes. • a Programmable Process Structure has built-in capa- bilities for robust communication with external (e.g. meteorological) databases and collaborating (e.g. metaheuristic optimization) programs. • the method supports the handling of the extraordi- nary uncertainty of coefficients of cost and some other parameters over a longer time horizon. The embedded elementary evaluations and elements of evaluation may be prepared for a multi-objective case involving a set of time-invariant natural objec- tive functions. 5. Conclusion Regardless of the well-established mathematical programming-based methods of the planning and op- eration of Hydrogen Supply Chains, problem-solving in terms of these systems might also require an easily modifiable, generic dynamic simulation-based approach. As was also shown by the studied illustrative exam- ple, the representation of energy storage corresponding to renewable resources requires a short (e.g. hourly, daily) 47(2) pp. 53–61 (2019) 60 VARGA, AZZARO-PANTEL, FLORES-PEREZ, AND CSUKÁS timescale as opposed to a years-long planning time hori- zon. The reason for this is that these large-scale, long-term supply chains are not only affected by uncertainly chang- ing demands, but by erratically changing supply, partly originating from natural resources, determined by errati- cally changing extreme weather conditions and a slowly changing climate. By considering these challenges with regard to the planning and operation of a Hydrogen Supply Chain, the non-conventional modeling and simulation methodology of Programmable Process Structures was implemented as developed in other multidisciplinary fields. The suggested Programmable Process Structures may contribute to the solution to these problems. Pro- grammable Process Structures of HSC models can be generated from the description of a (optionally geograph- ically determined and multiscale) process network and two general functional meta-prototypes, automatically. The case-specific, functional program prototypes can be derived from the two meta-prototypes. Simultaneously, the actual elements of the process network can be param- eterized and initialized according to the case-specific pro- totypes. The simulation of the model can be solved by a general-purpose kernel program. The method also allows different time steps for the various classes of elements to be used. The application of Programmable Process Structures was illustrated by a fictitious example of a simple Hydro- gen Supply Chain. Planned future work, in collaboration with HSC ex- perts, will focus on the application of a repeated for- ward/backward simulation-based iterative methodology to reach a reasonable compromise between changing de- mands and supplies. Furthermore, Programmable Process Structures will be combined with metaheuristic optimiz- ers for the multi-objective optimization of HSC planning. Notations GIS Geographic Information System GraphML A comprehensive and easy-to-use file format for graphs HSC Hydrogen Supply Chain MILP Mixed Integer Linear Programming MINLP Mixed-Integer Nonlinear Programming MPP Multi-period problems MOP Multi-objective problems PPS Programmable Process Structure Acknowledgement This research was supported by the TÉT_16-1-2016- 0116 French-Hungarian bilateral project. REFERENCES [1] Hydrogen Deployment. Special Section of Chem- ical Engineering Progress, Journal of the Ameri- can Institute of Chemical Engineers, August 2019, aiche.org/cep [2] Calloway, B.; McWhorter, S.; James, W.: Advance- ments in hydrogen deployment. In: Special Section of CEP Magazine, August 2019, https://www. aiche.org/resources/publications/cep/2019/ august/advancements-hydrogen-deployment [3] De-Leon Almaraz, S.: Multi-objective optimisation of a Hydrogen Supply Chain. Ph.D. Thesis, 2014, Université de Toulouse. [4] De-Leon Almaraz, S.; Azzaro-Pantel, C.: Chapter 4 – Design and Optimization of Hydrogen Sup- ply Chains for a Sustainable Future. In: Hydro- gen Economy, Eds. Scipioni, A., Manzardo, A., Ren, J., (Academic Press, 2017), pp. 85–120 DOI: 10.1016/B978-0-12-811132-1.00004-3 [5] Varga, M.; Csukás, B.: Meta-modeling of Hydrogen Supply Chains: a Programmable Structure based representation. In: (Azzaro-Pantel, C. ed.) Hydro- gen Supply Chain: Design, Deployment and Oper- ation, (Elsevier, 2018) DOI: 10.1016/B978-0-12-811197- 0.00014-2 [6] Show, K.Y.; Yan, Y.; Zong, C.; Guo, N.; Chang, J.S.; Lee, D.J.: State of the art and challenges of bio- hydrogen from microalgae, Bioresource Technol., 2019, 289, 121747 DOI: 10.1016/j.biortech.2019.121747. [7] Sivaram, V.: Taming the Sun: Innovations to Har- ness Solar Energy and Power the Planet, (MIT Press, 2018), pp. 392 ISBN-13: 978-0262037686 [8] Wang, J.; Yin, Y.: Principle and application of different pretreatment methods for enriching hydrogen-producing bacteria from mixed cultures, Int. J. Hydrogen Energ., 2017, 42(8), 4804–4823 DOI: 10.1016/j.ijhydene.2017.01.135. [9] Bakonyi, P.; Kumar, G.; Koók, L.; Tóth, G.; Rózsen- berszki, T.; Bélafi-Bakó, K.; Nemestóthy, N.: Mi- crobial electrohydrogenesis linked to dark fermen- tation as integrated application for enhanced biohy- drogen production: A review on process characteris- tics, experiences and lessons, Bioresource Technol., 2018, 251, 381–389 DOI: 10.1016/j.biortech.2017.12.064. [10] Govindan, K.; Fattahi, M.; Keyvanshokooh, E.: Supply chain network design under uncertainty: A comprehensive review and future research direc- tions, Eur. J. Oper. Res., 2017, 263, 108–141 DOI: 10.1016/j.ejor.2017.04.009 [11] Dagdougui, H.: Models, methods and approaches for the planning and design of the future hydrogen supply chain, Int. J. Hydrogen Energy, 2012, 37, 5318–5327 DOI: 10.1016/j.ijhydene.2011.08.041 [12] Maryam, S.: Review of modelling approaches used in the HSC context for the UK, Int. J. Hy- drogen Energy, 2017, 42(39), 24927–24938 DOI: 10.1016/j.ijhydene.2017.04.303 [13] Grossmann, I.E.; Guillen-Gosalbez, G.: Scope for the application of mathematical programming tech- niques in the synthesis and planning of sustainable Hungarian Journal of Industry and Chemistry aiche.org/cep https://www.aiche.org/resources/publications/cep/2019/august/advancements-hydrogen-deployment https://www.aiche.org/resources/publications/cep/2019/august/advancements-hydrogen-deployment https://www.aiche.org/resources/publications/cep/2019/august/advancements-hydrogen-deployment https://doi.org/10.1016/B978-0-12-811132-1.00004-3 https://doi.org/10.1016/B978-0-12-811132-1.00004-3 https://doi.org/10.1016/B978-0-12-811197-0.00014-2 https://doi.org/10.1016/B978-0-12-811197-0.00014-2 https://doi.org/10.1016/j.biortech.2019.121747. https://doi.org/10.1016/j.ijhydene.2017.01.135. https://doi.org/10.1016/j.biortech.2017.12.064. https://doi.org/10.1016/j.ejor.2017.04.009 https://doi.org/10.1016/j.ejor.2017.04.009 https://doi.org/10.1016/j.ijhydene.2011.08.041 https://doi.org/10.1016/j.ijhydene.2017.04.303 https://doi.org/10.1016/j.ijhydene.2017.04.303 ANALYSIS OF HYDROGEN SUPPLY CHAINS 61 processes, Comput. Chem. Eng., 2010, 34(9), 1365– 1376 DOI: 10.1016/j.compchemeng.2009.11.012 [14] Guillen-Gosalbez, G.; Mele, F.D.; Grossmann, I.E.: A bi-criterion optimization approach for the de- sign and planning of hydrogen supply chains for vehicle use; AIChE J., 2010, 56(3), 650–667 DOI: 10.1002/aic.12024 [15] Kudoh, Y; Ozawa, A: Life Cycle Assessment of Hy- drogen Supply Chain: A Case Study for Japanese Automotive Use. In: Azzaro-Pantel, C.: Hydrogen Supply Chains - Design, Deployment and Opera- tion (Academic Press, Elsevier, 2018), Chapter 15, pp. 499 DOI: 10.1016/B978-0-12-811197-0.09988-7 [16] Samsatli, S.; Staffell, I.; Samsatli, N.J.: Opti- mal design and operation of integrated wind- hydrogen-electricity networks for decarbonising the domestic transport sector in Great Britain, Int. J. Hydrogen Energy, 2016, 41, 447–475 DOI: 10.1016/j.ijhydene.2015.10.032 [17] Nurjanni, K.P.; Carvalho, M.S.; Costa, L.: Green supply chain design: A mathematical modeling approach based on a multi-objective optimization model, Int. J. Prod. Econ., 2017, 183, 421–432 DOI: 10.1016/j.ijpe.2016.08.028 [18] Azzaro-Pantel, C.: Hydrogen Supply Chains - De- sign, Deployment and Operation (Academic Press, Elsevier, 2018), pp. 588 DOI: 10.1016/B978-0-12-811197- 0.09988-7 [19] Varga, M.; Csukás, B.: Programmable Process Structures, Generated from a Network and from Functional Meta-Prototypes. In: 2017 AIChE An- nual Meeting Proceedings, Computing and Sys- tems Technology Division, Software Tools and Implementations for Process Systems Engineer- ing. Minneapolis (MN), USA, 2017. Paper 448e https://www.aiche.org [20] Varga, M.; Prokop, A.; Csukás, B.: Biosystem models, generated from a complex rule/reac- tion/influence network and from two functional- ity prototypes, BioSystems, 2017, 152, 24–43 DOI: 10.1016/j.biosystems.2016.12.005 [21] http://www.gprolog.org/ 47(2) pp. 53–61 (2019) https://doi.org/10.1016/j.compchemeng.2009.11.012 https://doi.org/10.1002/aic.12024 https://doi.org/10.1002/aic.12024 https://doi.org/10.1016/B978-0-12-811197-0.09988-7 https://doi.org/10.1016/j.ijhydene.2015.10.032 https://doi.org/10.1016/j.ijhydene.2015.10.032 https://doi.org/10.1016/j.ijpe.2016.08.028 https://doi.org/10.1016/j.ijpe.2016.08.028 https://doi.org/10.1016/B978-0-12-811197-0.09988-7 https://doi.org/10.1016/B978-0-12-811197-0.09988-7 https://www.aiche.org/conferences/aiche-annual-meeting/2017/proceeding/paper/448e-programmable-process-structures-generated-network-and-functional-meta-prototypes https://doi.org/10.1016/j.biosystems.2016.12.005 https://doi.org/10.1016/j.biosystems.2016.12.005 http://www.gprolog.org/ HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 63–70 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-21 INVESTIGATION OF MIXING IN TANKS OF A SPECIAL GEOMETRY BÁLINT LEVENTE TARCSAY*1, ATTILA EGEDY1, JANKA BOBEK2, AND DÓRA RIPPEL-PETHŐ2 1Department of Process Engineering, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY 2Department of Chem. Eng. Sci.s, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY Mixing is one of the most crucial processes in the chemical industry. Homogeneity is a requirement for all feedstocks and industrial products. The degree of mixing depends on the hydrodynamic properties of the fluid in the units. The Residence Time Distribution (RTD) was investigated in a tank of a special geometry. Mixing was investigated using various geome- tries of the tank by applying the Heaviside function in step-response experiments. After obtaining experimental results, the RTD function was calculated. The flow structure in the tank was approximated by fitting black-box transfer function models onto the RTD function of the system. Two general model structures were defined and their fitness compared. By evaluating the fitted models, a relationship was established between the flow structure in the tank and its geometry. Keywords: tanks with of special geometry, mixing, flow pattern, transfer function model 1. Introduction In industrial processes, the thorough homogeneity of ma- terials is a necessity to ensure the quality of the fin- ished products and the safety of the production process. Therefore, various industrial-scale techniques for the ho- mogenisation of materials of various densities and phases have been developed over the years. In the case of fluids, homogenization can be achieved through the application of mixers in a device. The optimal type and geometry of these stirred tanks can be approximated by specific mod- els which take the properties of the materials to be mixed and the desired degree of homogeneity into considera- tion [1]. In the case of the homogenisation of large vol- umes of liquids (petrochemical and radiochemical indus- tries), external recirculation can be more advantageous since it is considerably more inexpensive in terms of both installation and operation. The operation of stirred tanks has been thoroughly investigated by using mathematical and experimental methods, however, research on mixing achieved by external recirculation is lacking [2]. Opti- mization of the mixing process to achieve the most ho- mogenous product possible in an external recirculation tank can only be achieved after the flow structure of the unit is known. The aim of this research was to analyze the flow structure in a tank of special geometry, moreover, to determine the connection between the flow structure and geometry of the tank. This research forms the basis for future investigations where the flow structure and homogenization in a tank will be tested through an external recirculation mixing *Correspondence: tarcsaybalint95@gmail.com process. Mixing in stirred tanks can traditionally be ana- lyzed by step-response experiments. Step-response experiments are a widely accepted and simple way of determining the flow structure of process units [3]. During these experiments, the input of the sys- tem is a tracer fed according to the Heaviside function. The time evolution of the tracer concentration in the out- put of the tank holds information regarding the flow and mixing properties of the system. The Heaviside function is the integral of the Dirac delta function δ(t): S(t) = ∫ ∞ −∞ δ(t)dt = { 0 if t < t∗ 1 if t ≥ t∗ (1) used as an input in pulse-response experiments [4]. The Heaviside function, S(t), as a function of time takes on a value of 0 if the investigated time t is less than the time of the tracer feed t∗ and a value of 1 during and after being fed. After experimental data was collected on the response function of the system, the residence time distribution function (RTD) could be calculated. The RTD function approximates the probability that the residence time of the fed tracer agent is less than a specified residence time [4]. The RTD function, F , was calculated from F (t) = Rt −Rt−1 Rsteady −R0 , (2) whereRsteady andR0 are the values of the response func- tion in the steady state and at the beginning of the exper- iment, respectively. By analyzing the RTD curve in the studied unit, mix- ing and flow can be characterised and compared to the https://doi.org/10.33927/hjic-2019-21 mailto:tarcsaybalint95@gmail.com 64 TARCSAY, BOBEK, EGEDY, AND RIPPEL-PETHŐ Figure 1: Response functions of various ideal flow models [3]. flow structure of ideal hydrodynamic models [5, 6]. In this study, the ideal hydrodynamic models, which were used to describe the response, include the plug flow reac- tor (PFR) and tanks-in-series (TIS) model. The PFR model is useful for modelling systems with a continuous flow of cylindrical geometry. It works under the assumption that the flow through a PFR is perfectly mixed in the radial direction but not in the axial direc- tion. It is also assumed that the residence time of each fluid unit within the system is equal to the mean residence time, τ . Because of these conditions, the step response function of a system with PFR tendencies is identical to the step function which was used as an input but with a time delay equal to the mean residence time in the unit. An altered version of the PFR model is the laminar flow reactor (LFR) model. The LFR model assumes a laminar flow structure within the system. Therefore, the residence time of the fluid phase within the system has a distribu- tion. This causes the shape of the response function to be- come slightly altered compared to the step function. The TIS model is an extension of the continuous stirred tank reactor model (CSTR) which defines a system that is per- fectly mixed and, therefore, homogeneous. Upon entering the unit, the tracer is dispersed throughout the whole tank, therefore, it immediately appears in the output without a time delay. The TIS model examines a series of CSTR units. The RTD function in this model is dependent on the number of tanks that have been linked. The RTD function of the TIS model converges towards the PFR model as the number of units increases and towards the CSTR distri- bution as the number of units decreases. One tank returns the RTD function of a CSTR and as the number of tanks in series approaches infinity, the RTD function converges towards the PFR model. The characteristic RTD curves of the aforementioned ideal flow models are shown in Fig. 1 [3]. In this figure, t∗ denotes the duration of tracer in- jection and dt represents the time delay when the tracer appears in the outlet. In this investigation, black-box models were used to approximate the behavior of the tank. Black- and white- box modelling are two distinct modelling techniques to describe the traits of various systems. White-box models utilize well-known physical and chemical laws to charac- terize a system. These laws come in the form of differen- tial equations, algebraic equations, etc. The model defines the internal functioning of the system and calculates the specific output from the input by taking these into con- sideration. Generally in the case of complex mixing prob- lems, computational fluid dynamics (CFD) simulators are utilized to solve the white-box model of the system [7,8]. CFD simulations can determine the flow field in various units [9]. These simulators can be used to solve scientific prob- lems from a diverse spectrum of disciplines, e.g. aero- dynamics, hydrodynamics, ocean engineering, chemical engineering, etc. [10]. To solve these problems, various numerical methods have been invented over the years among which the finite element method (FEM) is one of the more frequently applied approaches [11]. However, in cases where the FEM is insufficiently advanced, alterna- tive numerical methods can also be implemented, e.g. the spectral element method (SEM) [12]. Their disadvantage is, however, that since they model the systems so thor- oughly, they require considerable computational power and computation time. Through these methods, it is pos- sible to analyze processes which are more complex or re- quire long experimentation times to be dissected. Despite this, these methods also require experimental data for val- idation. Black-box models, on the other hand, provide gener- alized means for the modelling of various problems [13]. They often employ empirical methods to calculate the output of a system from the input without thoroughly de- scribing the internal properties of the system [14]. The relationship between the input and output can be mod- elled using various techniques with the most popular be- ing neural networks as well as state-space and transfer function models [15]. In this study, the transfer function model was utilized to approximate the behavior of a tank of special geometry. 2. Experimental and the modelling pro- cess The geometry of the investigated tank is shown in Fig. 2. The tank consists of 5 cylindrical units integrated with each other. A removable wall was implemented which makes it possible to create the desired number of cylin- ders from the main tank. Therefore, the number of cylin- drical units can be varied between 1 and 5. During the investigation, geometries consisting of 1, 3 and 5 cylin- drical units were investigated. The residence time, cal- culated from the inlet volume flow rate and volume of the liquid, was 3 h for all geometries. Therefore, the in- let volume flow rate changed as the volume of liquid and number of cylindrical units varied. To obtain the response function of the system, a Heaviside function was used. The utilized tracer was a borax solution of low concentra- tion (6.24 gl−1). The tracer was fed into the tank which Hungarian Journal of Industry and Chemistry MIXING IN TANKS OF A SPECIAL GEOMETRY 65 Figure 2: The investigated tank. was filled with a high concentration of borax solution (29.16 gl−1). Conductivity measurements were used to monitor the change in concentration and obtain the response function. The experimental equipment is shown in Fig. 3. The ex- perimental conditions for the feed and residence time are shown in Table 1. The positions of the inlet and outlet during the exper- iments are shown in Table 2. The (0, 0, 0) coordinates of the system were defined as the bottom of the tank at the position of the outlet. The coordinates x, y and z cor- respond to the longitudinal, vertical and horizontal dis- tances, respectively. The RTD curves of the three different geometries of the tank were investigated. After experimentally obtain- ing the data from the RTD curves, transfer function mod- els were fitted to the experimental results. Two types of fitted functions were investigated. Both were black-box convolution models of the PFR and TIS. However, one contained only serial links between the perfectly mixed units, the other assumed parallel connections as well. Dif- ferential equations describe the change in concentration of the output of the tank as a function of the input. In var- ious ideal flow models, these equations were converted into algebraic equations by using the Laplace transform. After this the transfer function of the ideal flow mod- Figure 3: Setup of the experimental equipment (1– experimental tank, 2–tracer containment tank, 3– peristaltic pump, 4–buffer unit, 5–rotameter, 6–output containment tank). Table 1: Experimental parameters. Number of Cylinders τ (h) B (lh−1) 1 0.4 3 3 1.1 5 1.7 Table 2: Inlet and outlet positions. Number of Cylinders Inlet position Outlet position x (mm) y (mm) z (mm) x (mm) y (mm) z (mm) 1 75 100 0 0 2 0 3 225 100 0 0 2 0 5 375 100 0 0 2 0 els was defined as the ratio of the output (the Laplace transformed response function) to the input of the system (the step function of the tracer injection). If the product of multiple elementary transfer functions is computed, a TIS model can be generated. A general model structure for this network can be seen in Fig. 4. On the other hand, by combining the elementary transfer functions of ideal flow models, a network of flow models with parallel links can be created. A general model for this structure is shown in Fig. 5. In both cases, a PFR unit has been included in the structure, in general this is useful for modelling systems with a time delay. The time delay in the RTD function could point to either PFR tendencies in the flow structure or the presence of dead volumes within the tank. To distinguish which is present, further investigations with a CFD simulator are required. The general form of the transfer function that only uses serial connections is Gfitted = a bnsn + bn−1sn−1 + · · ·+ b0s0 e−tds, (3) where n denotes the order of the denominator and bn, a, and td are parameters of the transfer function. In the case of the model containing both parallel and serial links, the general transfer function is Gfitted= an−1s n−1 + · · ·+ a0s 0 bnsn + bn−1sn−1 + · · ·+ b0s0 e−tds, (4) where n is the order of the denominator and bn, an, and td are parameters of the transfer function. However, due to the parallel connections, the elementary transfer func- tions are not just multiplied but totalled. Therefore, the Figure 4: General model structure involving only serial links [5, 6]. 47(2) pp. 63–70 (2019) 66 TARCSAY, BOBEK, EGEDY, AND RIPPEL-PETHŐ Figure 5: General model structure involving both serial and parallel links [5, 6]. nominator also contains a higher order polynomial. The order of this polynomial can be one at most below the order of the denominator for the system to be suitable. In both cases, the fitted RTD curve Ffitted was cal- culated from the step response of the system which was characterized by the transfer function, Gfitted. The overall parameters of the transfer function were determined by numerically solving a minimization prob- lem. The minimization problem is expressed in min [E (n, a, b, d)] = ∑ [Fexp − Ffitted (n, a, b, d)] 2∑ (Fexp) 2 (5) where E is a function denoting the sum of the squared difference of the experimental and fitted RTD functions divided by the sum of the square of the experimental RTD curve. As such,E is a function similar to the relative error of the fitting which is dependent on the parameters of the transfer function, Gfitted. Inversely, the goodness of the fit, I , can be defined as the complementer of the fitting error: I (n, a, b, d) = 1− E (n, a, b, d) . (6) For the minimization of the function, the interior-point method was applied in MATLAB R2011. Constraints were defined for the lower and upper boundaries of the parameters, namely 0 and 100 for all parameters, respec- tively. The fitting algorithm was modified to favor trans- fer functions with the smallest possible order and smallest relative error. The algorithm can be seen in Fig. 6. The algorithm follows a general route for fitting both parallel and serial models. During the first step, a low- order (first-order) transfer function is created. After that the parameters of this transfer function are identified us- ing the interior point method. After obtaining the param- eters, the E function is calculated to evaluate the error of the fitting. If the error is less than a predefined threshold (Erropt), then the fitting will stop and the current trans- fer function be defined as the best possible fit. If the error of the fit exceeds the defined threshold, then the fitting continues. In this case, the order of the denominator is increased and the error of the higher order transfer func- tion evaluated as well. If the error of fitting is less than the lower order transfer function, then the fitting contin- ues until either a transfer function which satisfies the er- ror boundary is identified or the error of the fitting does not decrease significantly even though the order of the denominator increases. In the case of models which include parallel links, this algorithm contains an additional iteration step. Dur- ing this step, the algorithm investigates the fit of an nth order transfer function (if n ≥ 2) with varying orders of the nominator. To achieve this, the fitting is evaluated by applying increasing orders of the nominator. During this process, either the maximum order of the nominator (n − 1) is reached or a point where the increase in the order of the nominator is identified that does not signif- icantly impact the fitting error. In either case, fitting is stopped and the solution accepted as the best possible fit for a transfer function with order of the denominator, n. To evaluate the proposed fitting algorithm for both models, a custom-defined transfer function was analyzed. This reference transfer function (Gref is defined by Gref = 1 1.6s2 + 2.6s+ 1 e−5s (7) in the case of the model contains only serial links. The step response function for this transfer function was used as reference data. The extrema of function E were numerically determined by using the devised algo- rithm. The fitted transfer function is Gref = 1 1.58s2 + 2.59s+ 1 e−5s. (8) Fig. 7 shows the fitting of the RTD of the identified trans- fer function to the reference RTD function. The goodness of the fit was 98.2 %. The testing was also conducted in cases where paral- lel links are present besides serial links. In this case, the reference system is characterized by the transfer function: Gref = 8.6s+ 1 s4 + 2.5s3 + 11s2 + 10s+ 1 e−2s. (9) The fitting algorithm was augmented to identify both the optimal orders of the nominator and denominator. The fitted transfer function is Gfitted = 8.95s+ 1 s4 + 2.37s3 + 11.14s2 + 10.33s+ 1 e−2s. (10) The results of the fitting are displayed in Fig. 8. In this case, the goodness of the fit exceeded 99 %. Both Fig. 8 and the measurements of the parameters show that the algorithm is capable of precisely fitting both models, therefore, it is suitable for identifying the trans- fer functions of the system. After validating the fitting Hungarian Journal of Industry and Chemistry MIXING IN TANKS OF A SPECIAL GEOMETRY 67 Figure 6: The general applied fitted algorithm for a model containing both serial and parallel links. Figure 7: The results of the algorithm validation for fitting the model containing only serial connections. Figure 8: The results of the algorithm validation for fitting the model containing parallel and serial links. method, transfer functions were identified for the tank by using the experimental RTD curves. 3. Results and Analysis The results of the fittings in the case of the model with only serial connections can be seen in Figs. 9–11 for all three geometries of the tank that were investigated. The geometry consisting of 5 cylinders was origi- nally modeled using a second order transfer function that yielded the best fit for the experimental RTD curve (92 %). However, to compare the parameters of the fitting Figure 9: Results of the fitting of the serial model for 1 cylindrical unit. Figure 10: Results of the fitting of the serial model for 3 cylindrical units. and characterise the system, it was more favorable to fit a third order transfer function here as well since it describes the other two RTD curves best. In this case, the average goodness of the fit was 91 % with the goodness of the fit in the case of the geometry consisting of 5 cylindrical units (86 %) being the lowest. The numerical parameters of the fit are shown in Table 3. The results are as follows: to describe the system in the cases of 1 and 3 cylindrical units, a third order transfer function was used. In the case of 5 cylindrical units, a second order transfer function was optimal. In the case of the model containing both parallel and 47(2) pp. 63–70 (2019) 68 TARCSAY, BOBEK, EGEDY, AND RIPPEL-PETHŐ Figure 11: Results of the fitting of the serial model for 5 cylindrical units. Table 3: Fitted parameters for the model containing only serial links. Parameters N n a0 b3 b2 b1 b0 td (h) I (%) 1 3 1 0.1 0.3 1.0 1 0.9 93.43 3 3 1 0.2 0.6 1.5 1 0.6 94.00 5 3 1 0.3 0.7 2.0 1 0.7 85.76 5 2 1 - 0.9 2.3 1 0.7 92.00 serial links, the fitted curves are displayed in Figs. 12–14. In the case of the model containing both parallel and serial connections, the tendencies are similar to the re- sults of the fitting conducted using the model consisting of only serial links. However, the goodness of the fit was enhanced compared to the simpler models (over 99 % for all investigated geometries). In this case, the behavior of the tank consisting of 1 cylindrical unit was described by a transfer function where the orders of the nominator and denominator were 4 and 2, respectively. In the case of 3 and 5 cylindrical units, the orders of the denominator and nominator were 3 and 1, respectively. The order of the transfer function (the number of CSTR units in the TIS model) increases as the number of cylinders decreases. This is evident of an increase in the PFR/LFR tenden- cies of the functions which can also be seen in the fig- Figure 12: Results of the fitting of the parallel model for 1 cylindrical unit. Figure 13: Results of the fitting of the parallel model for 3 cylindrical units. Figure 14: Results of the fitting of the parallel model for 5 cylindrical units. ures. It should also be noted, however, that the fitting in all models was exceedingly high (over 95 % on average). This might suggest that the model algorithm is too pre- cise and the model was also fitted according to the exper- imental noise remaining in the data. In the absence of the noise, such tendencies might be clearer and the decrease in the order of the nominator between the different ge- ometries more prominent. Further investigations carried out by CFD simulations may clarify these concerns. The parameters of the fit are displayed in Table 4. 4. Conclusion In this paper, the velocity field of a tank of special ge- ometry was investigated. In order to determine the re- lationship between the flow profile and geometry of the tank, step-response experiments were conducted. The RTD curves of the system were measured by applying various geometries of the tank. By analyzing the exper- imental RTD curves, a model was proposed to be fitted onto the results. The transfer function of the system was optimized by minimalizing the squared difference of the experimental and fitted RTD curves. The interior-point method was applied to minimalize the difference. The following conclusions were drawn from the fitted RTD curves: Hungarian Journal of Industry and Chemistry MIXING IN TANKS OF A SPECIAL GEOMETRY 69 Table 4: Fitted parameters of the model containing both parallel and serial links. Parameters N n m a2 a1 a0 b4 b3 b2 b1 b0 td (h) I (%) 1 4 2 6 2.4 1 1 4.3 8 3.12 1 1.2 98.28 3 4 1 − 5.5 1 1 3.4 8.5 7.1 1 0.6 99.77 5 3 1 − 12 1 1 3.5 18 14.4 1 0.5 99.43 In the case of models containing only serial links, 1 and 3 cylindrical units could be optimally described us- ing a third order transfer function. The geometry consist- ing of 5 cylindrical unitscan be described by a second order transfer function. The average goodness of the fit in this case was 91 %. In the case of the model utilizing par- allel links, the RTD function of the tank consisting of 1 cylindrical unit was approximated by a transfer function where the orders of the denominator and nominator were 4 and 2, respectively. In the cases of 3 and 5 cylindrical units, the orders of the denominator and nominator were 3 and 1, respectively. The goodness of fit on average in the case of this model was over 99 %. It can be concluded that all geometries of the tank can be described using fourth order or lower transfer func- tions with increasing model parameters according to the number of cylinders. Further investigations to characterise the reasons for such behavior will be carried out by using CFD meth- ods. Simultaneously, different step-response experiments will be conducted to observe the behavior of the system. During these experiments, instead of feeding the solution of lower concentration into that of higher concentration within the tank, a positive step function will be used. As a result, the orders will be reversed with the solution of low concentration being in the tank and that of high con- centration in the feed. Notations Latin letters a nominator parameter of transfer function b denominator parameter of transfer function dt time delay [h] h order of the nominator of transfer function (1,n-1) n order of the denominator of transfer function [-] t time [h] td time delay of transfer function [h] x longitudinal coordinate [mm] y vertical coordinate [mm] z horizontal coordinate [mm] Capital letters B volume flow rate [lh−1] E error of fitting [%] F residence time distribution function G transfer function I goodness of fit [%] N number of cylinders R response function [gl−1] S step function [gl−1] Greek letters δ Dirac delta function τ mean residence time [h] Indices ∗ time of tracer injection d delay exp experimental fitted fitted ref reference steady steady 0 initial Acknowledgement This research was supported by the project EFOP-3.6.1- 16-2016-00015 “Smart Specialisation Strategy (S3) - comprehensive institutional development program at the University of Pannonia to promote sensible individual ed- ucation and career choices”. REFERENCES [1] Nienow, A. W., Edwards, M. F., Harnby, N. Mixing in the Process Industries. (Elsevier, Oxford, United Kingdom 1997) ISBN: 0- 7506-3760-9 [2] Blais, B., Lassaigne, M., Goniva, C., Fradette, L., Bertrand, F. Development of an unresolved CFD–DEM model for the flow of viscous sus- pensions and its application to solid–liquid mix- ing. J. Comp. Phys., 2016, 318, 201-221 DOI: 10.1016/j.jcp.2016.05.008 [3] Levenspiel, O. Chemical reaction engineering. Ind. Eng. Chem. Res., 1999, 38(11), 4140–4143 DOI: 10.1021/ie990488g [4] Danckwerts, P. V. Continuous flow systems: Distri- bution of residence times. Chem. Eng. Sci., 1995, 50(24), 3857–3866 DOI: 10.1016/0009-2509(96)81811-2 [5] Fazli-Abukheyli, R., Darvishi, P. Combination of axial dispersion and velocity profile in parallel tanks-in-series compartment model for prediction of residence time distribution in a wide range of non-ideal laminar flow regimes. Chem. Eng. Sci., 2019, 195, 531–540 DOI: 10.1016/j.ces.2018.09.052 [6] Haag, J., Gentric, C., Lemaitre, C., Leclerc, J. P. Modelling of Chemical Reactors: From Sys- temic Approach to Compartmental Modelling. Int. J. Chem. Reactor Eng., 2018, 16(8), 1–22 DOI: 10.1515/ijcre-2017-0172 [7] Blazek, J. Computational Fluid Dynamics: Prin- ciples and Applications. (Butterworth-Heinemann, 47(2) pp. 63–70 (2019) https://doi.org/10.1016/j.jcp.2016.05.008 https://doi.org/10.1016/j.jcp.2016.05.008 https://doi.org/10.1021/ie990488g https://doi.org/10.1021/ie990488g https://doi.org/10.1016/0009-2509(96)81811-2 https://doi.org/10.1016/j.ces.2018.09.052 https://doi.org/10.1515/ijcre-2017-0172 https://doi.org/10.1515/ijcre-2017-0172 70 TARCSAY, BOBEK, EGEDY, AND RIPPEL-PETHŐ Oxford, United Kingdom 2015) ISBN: 978-0-08- 099995-1 [8] Toro, E. F. Riemann Solvers and Numerical Meth- ods for Fluid Dynamics. (Springer, Heidelberg, Ger- many 1999) ISBN: 978-3-662-03492-7 [9] Pawlowski, S., Nayak, N., Meireles, M., Portugal, C. A. M., Velizarov, S., Crespo, J. G. CFD mod- elling of flow patterns, tortuosity and residence time distribution in monolithic porous columns recon- structed from X-ray tomography data. Chem. Eng. J., 2018, 350, 757–766 DOI: 10.1016/j.cej.2018.06.017 [10] Delafosse, A., Collignon, M. L., Calvo, S., Delvi- gne, F., Crine, M., Thonart, P., Toye, D.: CFD-based compartment model for description of mixing in bioreactors. Chem. Eng. Sci., 2014, 106, 76–85 DOI: 10.1016/j.ces.2013.11.033 [11] Zienkiewicz, O. C., Taylor, R. L., Zhu, J. Z. The Fi- nite Element Method: Its Basis and Fundamentals. (Elsevier, Oxford, United Kingdom, 2005) ISBN: 0- 7506-6320-0 [12] Solin, P., Segeth, K., Dolezel, I. Higher-Order Fi- nite Element Methods. (Chapman and Hall/CRC, New York, 2003) ISBN: 978-0-42-920527-9 [13] Witt, C., Bux, M., Gusew, W., Leser, U. Predictive performance modeling for distributed batch pro- cessing using black box monitoring and machine learning. Information Systems 2019, 82, 33–52 DOI: 10.1016/j.is.2019.01.006 [14] Suykens, Johan A. K.; Vandewalle, Joos P. L. (ed.) Nonlinear Modelling: Advanced Black-Box Tech- niques. (Springer Science and Business Media, Leu- ven, 2012) ISBN: 978-1-4613-7611-8 [15] Bunge, M. A general black box theory. Philosophy of Science, 1963, 30(4), 346–358 DOI: 10.1086/287954 Hungarian Journal of Industry and Chemistry https://doi.org/10.1016/j.cej.2018.06.017 https://doi.org/10.1016/j.ces.2013.11.033 https://doi.org/10.1016/j.ces.2013.11.033 https://doi.org/10.1016/j.is.2019.01.006 https://doi.org/10.1016/j.is.2019.01.006 https://doi.org/10.1086/287954 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 71–75 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-22 ASSESSMENT OF THE ECOTOXICITY OF NANOPLASTICS NÓRA KOVÁTS*1, BETTINA ECK-VARANKA1, ZSÓFIA BÉKÉSSY1, DORINA DIÓSI1, KATALIN HUBAI1, AND JÁNOS KORPONAI2 1Institute of Environmental Sciences, University of Pannonia, 8200 Veszprém, Egyetem u. 10, HUNGARY 2Department of Environmental Science, Sapientia Hungarian University of Transylvania, 400193 Cluj-Napoca, Calea Turzii 4, ROMANIA The presence of micro- and nanoplastics in aquatic environments (including freshwater and marine ecosystems as well as their sediments) is becoming an increasingly serious problem worldwide. A wide range of studies have addressed the ecological effects these particles pose on biota. The main exposure pathway are food chains, e.g. under laboratory conditions these particles accumulate in the brain tissues of fish that feed on zooplankton causing brain damage. These studies, however, report mainly on the physical effects. In order to establish actual ecotoxicological effects, nanoplastics (50 nm in diameter) were assessed using the Vibrio fischeri bioluminescence inhibition bioassay (VFBIA). Our results showed that even environmentally relevant concentrations might trigger ecotoxicological effects. This study can be con- sidered to be a first screening, however, results indicate the need for more complex testing on a battery of aquatic test organisms. Keywords: nanoplastic; ecotoxicity; Vibrio fischeri ; kinetic assay 1. Introduction Aquatic environments contaminated by plastic litter are an emerging problem. Remote, pristine mountainous ar- eas are even contaminated by atmospheric microplastic deposition [1]. Polymer particles < 5 mm in diameter are defined as microplastics (MP) and may be derived directly from the use of industrial pellets or indirectly from the degradation and fragmentation of plastic parti- cles [2]. Polystyrene was proven to degrade into micro- and nanoplastics under laboratory conditions [3]. High levels of contamination have been reported in both ma- rine and freshwater habitats [4, 5]. Micro- and nanoplas- tics (NP) can float freely in bodies of water or be de- posited as sediments. The highest risk associated with these particles is their ingestion, which occurs at different levels in the aquatic food chain. Jabeen et al. [6], for example, listed approximately 150 different fish species where ingestion and accumulation have been reported. Particles can also progress upwards in the trophic levels of the food chain, i.e. fish can be exposed to the ingestion of zooplankton which is not able to discriminate between different food sources and consumes micro- and nanoplastics [7]. An experimental study showed that in fish exposed to NPs via the food chain, these particles caused brain damage and behavioural disorders as a result of accumulation in *Correspondence: kovats@almos.uni-pannon.hu brain tissues [8]. Biomagnification may also affect food safety and human health, though certain knowledge gaps exist in this field [9]. Ingestion may actually lead to starvation and even- tually the impairment of their physical condition. Un- der laboratory conditions, Daphnia magna exposed to polystyrene nanoparticles (PS-NP) exhibited reduced body size and severe alterations in terms of reproduc- tion [10]. D. magna is a widely studied species due to its key role in the aquatic food chain. It was shown to ingest nano- and microplastics (20 nm to 70 µm in diam- eter) from water [11]. In a laboratory study by Mattsson et al. [8], particles 52 nm in diameter elucidated the most severe effects. Cui et al. [12] exposed D. galeata to PS- NPs (5 mg/l, 52 nm in diameter) and detected a significant mortality rate after 2 days of exposure until the end of the study which lasted 5 days. Although a standard ecotoxi- cological test was conducted in this case, the mechanisms of mortality are still unclear: physical contact might have led to a reduction in the survival rate. In general, most ecotoxicological studies have used relatively high concentrations. Manfra et al. [13] investi- gated the impact of green fluorescently labeled carboxy- lated polystyrene nanoparticles of 40 nm in diameter with various surface charges on the marine rotifer Brachionus plicatilis. It was found that anionic PS-NPs did not elu- cidate mortality within the range of concentrations tested (0−50 µg/ml), while cationic PS-NPs caused mortality at https://doi.org/10.33927/hjic-2019-22 mailto:kovats@almos.uni-pannon.hu 72 KOVÁTS, ECK-VARANKA, BÉKÉSSY, DIÓSI, HUBAI, AND KORPONAI concentrations ≥ 2.5 µg/ml. Changes in oxidative stress enzymes were detected within the concentration range of 10−20 µg/ml in different organisms, e.g. the rotifer Bra- chionus koreanus and the marine copepod Paracyclop- ina nana [14]. The same concentration, 10 µg/ml, was reported to cause 40 % growth inhibition in the green mi- croalga Dunaliella tertiolecta [15]. In order to distinguish real (eco)toxicological effects from physical damage, a test based on the biolumines- cence inhibition of the marine bacterium Vibrio fischeri was selected. The species has been reclassified as Ali- ivibrio fischeri [16], however, as most standards and even recent papers from the literature still use the name V. fis- cheri, it will be used hereinafter. Bioluminescence is regulated by the enzyme system NAD(P)H:FMN oxidoreductase-luciferase. In toxic en- vironments, enzyme inhibition is reflected by a rapid de- crease in the luminous emittance of the bacterium. The reduction in light intensity is easy to measure as it is pro- portional to the strength of the toxicant, therefore, pro- vides a quantifiable endpoint. This test has been used in various environmental matrices [17–19]. Lappalainen et al. [20,21] developed a special version of the test which was later standardised (ISO 21338:2010: Water quality - Kinetic determination of the inhibitory ef- fects of sediment, other solids and coloured samples on the light emission of Vibrio fischeri (kinetic luminescent bacteria test)) in which bacteria are kept in suspension in direct contact with potentially toxic solid particles. Lumi- nescence readings were taken when the test commenced and the light intensity continuously monitored over the first 30 secs after the sample had been mixed with the bacteria. The light output pattern, therefore, might al- ready provide some indication of the expected toxicity of the sample [22]. The light intensity was measured once more after the pre-set exposure time (5, 15 or 30 mins as per standard). Toxicity values are normally expressed as EC50 and EC20, i.e. concentrations causing lumines- cence inhibitions of 50 and 20 % in this assay, respec- tively. 2. Materials and Methods In our experiments, the Ascent luminometer (Flash sys- tem, marketed by Aboatox, Finland) was used. A suspen- sion of the test bacteria (NRRL B-11177) was prepared in accordance with manufacturer instructions (Hach Lange GmbH). Polystyrene particles with a nominal diameter of 50 nm were used as a sample (supplier Thermo Fisher Sci- entific). As no comparative data were available on the potentially toxic concentration, a range-finding concen- tration series was set [23]. Three initial sample concen- trations were selected (1 g/l, 1 µg/l and 1 ng/l), which were further diluted, the number of dilutions was 11 (the number of concentrations the 96-multiwell plate permits) and the dilution ratio 1 : 2. Table 1: The measured EC20 values of the polystyrene nanoparticles. Concentration 1 g/l 1 µg/l 1 ng/l EC20 5.2 17.31 30.51 The Vibrio fischeri strain NRRL B-11177 was recon- stituted by adding the contents of one vial of +4 ◦C 1243- 551 Reagent Diluent. The reconstituted reagent was equi- librated at +4 ◦C for 30 min. Then the reagent was sta- bilised at +15 ◦C for 30 mins before being pipetted into the wells. Luminescence readings were taken when the test commenced, Time0, and after the pre-set exposure time of 30 mins, Time30. The luminescence inhibition of each sample was calculated as follows: CF = IC30/IC0 INH \% = 100 $-$ 100 x (IT30 / CF x IT0) where CF = correction factor IC30 = luminescence intensity of the control sample after the contact time (30 mins) in the RLU IC0 = initial luminescence intensity of the control sample in the RLU IT30 = luminescence intensity of the test sample after the contact time (30 mins) in the RLU IT0 = initial luminescence intensity of the test sample in the RLU EC20 values were calculated using the Ascent software, also developed by Aboatox Oy. 3. Results and Discussion Table 1 shows the ecotoxicity expressed in EC20, i.e. the calculated concentration of the sample that caused 20 % bioluminescence inhibition. Fig. 1 illustrates the biolumi- nescence inhibition during the first 30 secs for the sam- ples of 1 g/l and 1 ng/l in concentration. EC20 (or in some cases, EC10) are considered thresholds for the estimation of the lowest observed effec- tive concentration [24], i.e. the sample is normally con- sidered (eco)toxic if the elucidated effect exceeds 20 %. These results show that the V. fischeri bioassay de- tected a measurable degree of toxicity even at a concen- tration of 1 ng/l. Booth et al. [25] used the non-kinetic version of this bioassay (Microtox®), however, in their study, the calculated toxic concentration exceeded the range of concentrations studied (0.001−1000 mg/l). The same negative effect was reported by Casado et al. [26]. The higher degree of (detectable) toxicity in our study might be explained by the differences in the test system used. While Microtox® is a non-kinetic test, the Flash system (Ascent luminometer) was especially developed to test the toxicity of different suspensions or samples containing solid particles. The Ascent luminometer uses a 96-multiwell microplate. A specific feature of it is that Hungarian Journal of Industry and Chemistry ASSESSMENT OF THE ECOTOXICITY OF NANOPLASTICS 73 (a) (b) Figure 1: Kinetic diagram of the 1 g/l (a) and 1 ng/l (b) samples. The light output is recorded over the first 30 sec- onds. After the peak, toxicity causes a rapid reduction in the light output, on the other hand, it remains constant dur- ing the control. The two columns show the two replicates. E1-F1/G1-H1 (left): control. E2-F2/G2-H2 (right): sam- ple, maximum concentration. during luminescence readings, the microplate is continu- ously shaken by the instrument, resulting in the resuspen- sion of particles. According to our results, environmentally relevant concentrations might already pose ecotoxic effects. Ac- tual environmental concentrations are relatively difficult to compare and assess, mostly due to difficulties in sam- pling and the lack of standardized sampling methodolo- gies [27,28]. Indicative data are available: e.g. microplas- tic concentrations of 0.4 − 34 ng/l in bodies of freshwa- ter in the USA [29] or 0.51 mg/l in marine environments [10]. However, in real-world environments, even higher levels of toxicity can be expected as particles might absorb organic pollutants from the surrounding water [30], including highly toxic pesticides or polychlorinated biphenyls (PCBs) [31]. Though their bioavailability is still questionable [32], Batel et al. [33] conducted a lab- oratory study on microplastics and one polycyclic aro- matic hydrocarbon (PAH), benzo[a]pyrene (BaP). It was demonstrated that BaP adsorbed on microplastics and was transferred via an artificial food chain. These par- ticles might also possess inherent toxicity due to the use of additives during manufacturing processes [34]. 4. Conclusions It is a well-known paradigm in ecotoxicology that the sensitivity of different test organisms to a particular chemical varies, therefore, the V. fischeri test can be re- garded as a first screening. The bioluminescence inhibi- tion assay is an acute test that uses a maximum expo- sure of only 30 minutes. Naturally, chronic effects can- not be extrapolated from these results. However, the fact that the tested nanoplastics have already elucidated eco- toxicological effects in environmentally relevant concen- trations emphasises the need for more complex ecotoxi- cological testing involving a properly selected battery of test organisms. In addition to widely used aquatic test or- ganisms such as the aforementioned Daphnia magna, an ideal candidate could be the Caenorhabditis elegans test. It is a standardised bioassay using a sediment-dwelling, widely distributed nematode. However, in order to distin- guish physical damage from toxic effects, the measure- ment of changes in oxidative stress enzymes can be use- ful no matter which test organism is applied. Acknowledgement This work was funded by the BIONANO_GINOP-2.3.2- 15-2016-00017 project. REFERENCES [1] Allen, S.; Allen, D.; Phoenix, V.R.; Le Roux, G.; Jiménez, P.D.; Simonneau, A.; Binet, S.; Galop, D.: Atmospheric transport and deposition of microplas- tics in a remote mountain catchment, Nature Geo- science 2019, 12(5), 339–344 DOI: 10.1038/s41561-019- 0335-5 [2] Mintenig, S.M.; Dris, R.; Imhof, H.; Sanchez, W.; Gasperi, J.; Galgani, F.; Tassin, B.; Laforsch, C.: Beyond the ocean: Contamination of freshwater ecosystems with (micro-)plastic particles, Environ. Chem. 2015, 12(5) 539–550 DOI: 10.1071/en14172 [3] Lambert, S.; Wagner, M.: Characterisation of nanoplastics during the degradation of polystyrene, Chemosphere 2016, 145 265–268 DOI: 10.1016/j.chemosphere.2015.11.078 [4] Cole, M.; Lindeque, P.; Halsband, C.; Galloway, T.S.: Microplastics as contaminants in the marine environment: a review, Mar. Pollut. Bull. 2011, 62 2588–2597 DOI: 10.1016/j.marpolbul.2011.09.025 [5] Dris, R.; Gasperi, J.; Rocher, V.; Saad, M.; Renault, N.; Tassin, B.: Microplastic contamination in an ur- ban area: a case study in Greater Paris, Environ. Chem. 2015 12(5), 592–599 DOI: 10.1071/en14167 [6] Jabeen, K.; Su, L.; Li, J.; Yang, D.; Tong, C.; Mu, J.; Shi, H.: Microplastics and mesoplastics in fish from coastal and fresh waters of China, Environ. Pollut. 2017, 221, 141–149 DOI: 10.1016/j.envpol.2016.11.055 [7] Chae, Y.; An, Y.-J.: Effects of micro- and nanoplas- tics on aquatic ecosystems: Current research trends and perspectives, Mar. Pollut. Bull. 2017, 124(2), 624–632 DOI: 10.1016/j.marpolbul.2017.01.070 [8] Mattsson, K.; Johnson, E.V.; Malmendal, A.; Linse, S.; Hansson, L.-A.; Cedervall, T.: Brain damage and behavioural disorders in fish induced by plastic nanoparticles delivered through the food chain, Sci. Rep. 2017, 7(1), 11452 DOI: 10.1038/s41598-017-10813-0 47(2) pp. 71–75 (2019) https://doi.org/10.1038/s41561-019-0335-5 https://doi.org/10.1038/s41561-019-0335-5 https://doi.org/10.1071/en14172 https://doi.org/10.1016/j.chemosphere.2015.11.078 https://doi.org/10.1016/j.marpolbul.2011.09.025 https://doi.org/10.1071/en14167 https://doi.org/10.1016/j.envpol.2016.11.055 https://doi.org/10.1016/j.marpolbul.2017.01.070 https://doi.org/10.1038/s41598-017-10813-0 74 KOVÁTS, ECK-VARANKA, BÉKÉSSY, DIÓSI, HUBAI, AND KORPONAI [9] Barboza, L.G.A.; Vethaak, A.D.; Lavorante, B.R.B.O.; Lundebye, A.-K.; Guilhermino, L.: Marine microplastic debris: An emerging is- sue for food security, food safety and human health, Mar. Pollut. Bull. 2018, 133, 336–348 DOI: 10.1016/j.marpolbul.2018.05.047 [10] Besseling, E.; Wang, B.; Lürling, M.; Koel- mans, A.A.: Nanoplastic Affects Growth of S. obliquus and Reproduction of D. magna, Envi- ron. Sci. Technol. 2014, 48(20), 12336-12343 DOI: 10.1021/es503001d [11] Rosenkranz, P.; Chaudhry, Q.; Stone, V; Fernandes, T.F.: A comparison of nanoparticle and fine particle uptake by Daphnia magna, Environ. Toxicol. Chem. 2009, 28(10), 2142–2149 DOI: 10.1897/08-559.1 [12] Cui, R.; Kim, S.W.; An, Y.-J.: Polystyrene nanoplastics inhibit reproduction and induce abnor- mal embryonic development in the freshwater crus- tacean Daphnia galeata, Sci. Rep. 2017, 7(1), 12095 DOI: 10.1038/s41598-017-12299-2 [13] Manfra, L.; Rotini, A.; Bergami, E.; Grassi, G.; Faleri, C.; Corsi, I.: Comparative ecotoxicity of polystyrene nanoparticles in natural seawater and reconstituted seawater using the rotifer Brachionus plicatilis, Ecotoxicol. Environ. Saf. 2017, 145, 557– 563 DOI: 10.1016/j.ecoenv.2017.07.068 [14] Prokić, M.D.; Radovanović, T.B.; Gavrić, J.P.; Fag- gio, C.: Ecotoxicological effects of microplastics: Examination of biomarkers, current state and future perspectives, Trends. Analyt. Chem. 2019, 111, 37– 46 DOI: 10.1016/j.trac.2018.12.001 [15] Gambardella, C.; Morgana, S.; Bramini, M.; Ro- tini, A.; Manfra, L.; Migliore, L.; Piazza, V.; Gar- aventa, F.; Faimali, M.: Ecotoxicological effects of polystyrene microbeads in a battery of ma- rine organisms belonging to different trophic lev- els, Mar. Environ. Res. 2018, 141, 313–321 DOI: 10.1016/j.marenvres.2018.09.023 [16] Urbanczyk, H.; Ast, J.; Higgins, M.J.; Carson, J.; Dunlap, P.V.: Reclassification of Vibrio fischeri, Vibrio logei, Vibrio salmonicida and Vibrio woda- nis as Aliivibrio fischeri gen. nov., comb. nov., Aliivibrio logei comb. nov., Aliivibrio salmonicida comb. nov. and Aliivibrio wodanis comb. nov., Int. J. Syst. Evol. Microbiol. 2007, 57(12), 2823–2829 DOI: 10.1099/ijs.0.65081-0 [17] Girotti, S.; Ferri, E.N.; Fumo, M.G.; Maiolini, E.: Monitoring of environmental pollutants by biolumi- nescent bacteria, Anal. Chim. Acta 2008, 608(1), 2– 29 DOI: 10.1016/j.aca.2007.12.008 [18] Ma, X.Y.; Wang, X.C.; Ngo, H.H.; Guo, W.; Wu, M.N.; Wang, N.: Bioassay based luminescent bac- teria: Interferences, improvements, and applica- tions, Sci. Total Environ. 2014, 468–469, 1–11 DOI: 10.1016/j.scitotenv.2013.08.028 [19] Abbas, M.; Adil, M.; Ehtisham-ul-Haque, S.; Mu- nir, B.; Yameen, B.; Ghaffar, A. et al.: Vibrio fis- cheri bioluminescence inhibition assay for ecotoxi- city assessment: A review, Sci. Total Environ. 2018, 626, 1295–1309 DOI: 10.1016/j.scitotenv.2018.01.066 [20] Lappalainen, J.; Juvonen, R.; Vaajasaari, K.; Karp, M.: A new flash method for measuring the toxicity of solid and colored samples, Chemosphere 1999, 38(5), 1069–1083 DOI: 10.1016/s0045-6535(98)00352-x [21] Lappalainen, J.; Juvonen, R.; Nurmi, J.; Karp, M.: Automated color correction method for Vibrio fis- cheri toxicity test. Comparison of standard and ki- netic assays, Chemosphere 2001, 45(4–5), 635–641 DOI: 10.1016/s0045-6535(00)00579-8 [22] Mortimer, M.; Kasemets, K.; Heinlaan, M.; Kurvet, I.; Kahru, A.: High throughput kinetic biolumines- cence inhibition assay for study of toxic effects of nanoparticles, Toxicol. in Vitro 2008, 22(5), 1412– 1417 DOI: 10.1016/j.tiv.2008.02.011 [23] USEPA (2000) Method Guidance and Recommen- dations for Whole Effluent Toxicity (WET) Testing (40 CFR Part 136). EPA 821-B-00-004. U.S. Envi- ronmental Protection Agency, Office of Water [24] Ventura, S.P.M.; Silva, F.A.; Gonçalves, A.M.M.; Pereira, J.P.; Gonçalves, J.P.; Coutinho, J.A.P.: Eco- toxicity analysis of cholinium-based ionic liquids to Vibrio fischeri marine bacteria, Ecotoxicol. Environ. Saf. 2014, 102, 48–54 DOI: 10.1016/j.ecoenv.2014.01.003 [25] Booth, A.M.; Hansen, B.H.; Frenzel, M.; Johnsen, H.; Altin, D.: Uptake and Toxicity of Methylmethacrylate-Based Nanoplastic Particles in Aquatic Organisms, Environ. Toxicol. Chem. 2016, 35(7), 1641–1649 DOI: 10.1002/etc.3076 [26] Casado, M.P.; Macken, A.; Byrne, H.J.: Ecotox- icological assessment of silica and polystyrene nanoparticles assessed by a multitrophic test battery, Environ. Int. 2013, 51, 97–105 DOI: 10.1016/j.envint.2012.11.001 [27] Rocha-Santos, T.; Duarte, A.C.: A critical overview of the analytical approaches to the occurrence, the fate and the behavior of microplastics in the en- vironment, Trends Analyt. Chem. 2015, 65, 47–53 DOI: 10.1016/j.trac.2014.10.011 [28] Twiss, M.R.: Standardized methods are required to assess and manage microplastic contamination of the Great Lakes system, J. Great Lakes Res. 2016, 42(5), 921–925 DOI: 10.1016/j.jglr.2016.07.032 [29] Eriksen, M.; Mason, S.; Wilson, S.; Box, C.; Zellers, A.; Edwards, W. et al.: Microplastic pollu- tion in the surface waters of the Laurentian Great Lakes, Mar. Pollut. Bull. 2013, 77(1–2), 177–182 DOI: 10.1016/j.marpolbul.2013.10.007 [30] Mato, Y.; Isobe, T.; Takada, H.; Kanehiro, H.; Ohtake, C.; Kaminuma, T.: Plastic resin pellets as a transport medium for toxic chemicals in the marine environment, Environ. Sci. Technol. 2001, 35(2), 318–324 DOI: 10.1021/es0010498 [31] Cedervall, T.; Hansson, L.A.; Mattsson, K.: Nano- plastics in the aquatic environment, Environ. Sci. Process. Impacts 2015, 17(10), 1712–1721 DOI: 10.1039/c5em00227c Hungarian Journal of Industry and Chemistry https://doi.org/10.1016/j.marpolbul.2018.05.047 https://doi.org/10.1016/j.marpolbul.2018.05.047 https://doi.org/10.1021/es503001d https://doi.org/10.1021/es503001d https://doi.org/10.1897/08-559.1 https://doi.org/10.1038/s41598-017-12299-2 https://doi.org/10.1016/j.ecoenv.2017.07.068 https://doi.org/10.1016/j.trac.2018.12.001 https://doi.org/10.1016/j.marenvres.2018.09.023 https://doi.org/10.1016/j.marenvres.2018.09.023 https://doi.org/10.1099/ijs.0.65081-0 https://doi.org/10.1016/j.aca.2007.12.008 https://doi.org/10.1016/j.scitotenv.2013.08.028 https://doi.org/10.1016/j.scitotenv.2013.08.028 https://doi.org/10.1016/j.scitotenv.2018.01.066 https://doi.org/10.1016/s0045-6535(98)00352-x https://doi.org/10.1016/s0045-6535(00)00579-8 https://doi.org/10.1016/j.tiv.2008.02.011 https://doi.org/10.1016/j.ecoenv.2014.01.003 https://doi.org/10.1002/etc.3076 https://doi.org/10.1016/j.envint.2012.11.001 https://doi.org/10.1016/j.envint.2012.11.001 https://doi.org/10.1016/j.trac.2014.10.011 https://doi.org/10.1016/j.jglr.2016.07.032 https://doi.org/10.1016/j.marpolbul.2013.10.007 https://doi.org/10.1021/es0010498 https://doi.org/10.1039/c5em00227c https://doi.org/10.1039/c5em00227c ASSESSMENT OF THE ECOTOXICITY OF NANOPLASTICS 75 [32] Beckingham, B.; Ghosh, U.: Differential bioavail- ability of polychlorinated biphenyls associated with environmental particles: Microplastic in compari- son to wood, coal and biochar, Environ. Poll. 2017, 220, 150–158 DOI: 10.1016/j.envpol.2016.09.033 [33] Batel, A.; Linti, F.; Scherer, M.; Erdinger, L.; Braunbeck, T.: The transfer of benzo[a]pyrene from microplastics to Artemia nauplii and further to ze- brafish via a trophic food web experiment–CYP1A induction and visual tracking of persistent organic pollutants, Environ. Toxicol. Chem. 2016, 35(7), 1656–1666 DOI: 10.1002/etc.3361 [34] Avio, C.G.; Gorbi, S.; Regoli, F.: Plastics and mi- croplastics in the oceans: From emerging pollutants to emerged threat, Mar. Environ. Res. 2017, 128, 2– 11 DOI: 10.1016/j.marenvres.2016.05.012 47(2) pp. 71–75 (2019) https://doi.org/10.1016/j.envpol.2016.09.033 https://doi.org/10.1002/etc.3361 https://doi.org/10.1016/j.marenvres.2016.05.012 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(2) pp. 77–83 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-23 BIOCOMPONENT-BASED QUALITY IMPROVEMENT OPPORTUNITIES OF BINDERS FOR ROAD CONSTRUCTION GRÉTA TÁLOSI1, PÉTER GERGÓ1, AND ANDRÁS HOLLÓ *1,2 1Department of MOL Hydrocarbon and Coal Processing, Institute of Chemical and Process Engineering, University of Pannonia, Egyetem u. 10, Veszprém, H-8200, HUNGARY 2Downstream R&D, MOL Plc., Olajmunkás u. 2, Százhalombatta, H-2440, HUNGARY In our experimental work, stabilised pinewood-based fast pyrolysis bio-oil and the linear block copolymer styrene butadi- ene styrene (SBS) were used as additives in bitumen used for road construction (Penetration Grade 50/70) to produce a higher-performance asphalt binder. Our aim was to investigate the modifying effect of the biocomponent on bituminous binders and prepare a comparative analysis. In order to characterize our samples, conventional and rheological mea- surements were performed. It was concluded that the biocomponent by itself cannot provide a favourable bituminous binder with beneficial mechanical properties, however, a favourable solution may be achieved by combining it with the block copolymer SBS. Based on our test results, stabilised pinewood-based fast pyrolysis bio-oil could be blended with the examined bituminous binders up to concentrations of 10 w/w% resulting in good bitumen quality. Keywords: bitumen 50/70, polymer modified bitumen, biocomponent, quality improvement 1. Introduction The stricter demands of road construction has brought about the application of new and high performance, al- beit expensive, binders (e.g. polymer modified bitumen). Thus, the development of alternative, high-quality but cheaper bitumen for road construction has become an in- tensively researched area [1–5]. Another aspect is that in line with sustainable devel- opment, the use of alternative materials derived from non- fossil sources, e.g. biomass or waste, can be a long-term option for the partial replacement of bitumens derived from crude oil. A further issue is that the production of low carbon fuels in the future might reduce the demand for refining crude oil, which could result in a lack of sup- ply in the bitumen market as well. Many researchers have found that different biomass-derived oils could be suit- able components for blending, substituting for and mod- ifying bituminous binders [6]. The raw materials of the processed biomass that produce bio-oils include microal- gae [7], animal waste [8, 9], maize [10], garden waste [11], tea residue [12], coffee residue [13], rapeseed [14], and soybean [15]. The utilisation of these bio-oils can be beneficial following proper treatment, since these materi- als could increase the rutting performance of bituminous binders. These studies show that the utilisation of renew- able materials can be very promising [16]. Many researchers have studied the modifying effect *Correspondence: ahollo@mol.hu of different wood-based bio-oils. Peralta et al. examined the application of red oak residues as direct alternative binders and to reduce the environmental impact of bitu- minous binders [17–19]. Grilli et al. used pinewood bio- oil [20], Jiménez del Barco-Carrión et al. used pine resin [21, 22] and Lei et al. used bio-oil derived from wood as a rejuvenator [23]. Raouf and Williams investigated the utilisation options of oak wood bio-oil [24]. Cooper et al. [25], Ball et al. [26], and Bearsley and Haverkamp blended bio-oil originated pine derived tall oil as a bitu- men extender [27, 28]. Gondim et al. studied the effect of plant sap from a “Petroleum Plant” [29]. Yang et al. examined bio-oils derived from waste wood resources as bitumen modifiers and extenders [30]. In summary, bio-oils typically have a softening effect, which improves their performance at low and intermedi- ate temperatures but impairs their high-temperature per- formance. This negative impact can be mitigated, e.g. by the addition of crumb rubber or other polymers. The application of stabilised wood-based pyrolysis bio-oil as an additive for bitumen binders has been in- tensively researched but remains a field of future bitumen production that is yet to be clarified. The utilisation of polymer modified bitumen is widespread but the modify- ing effect of stabilised wood-based pyrolysis bio-oil on polymer modified bitumen has yet to be clearly defined. Thus, our target was to investigate the effect of stabilised pinewood-based fast pyrolysis bio-oil on the character- istics of different binders, such as road pavement bitu- https://doi.org/10.33927/hjic-2019-23 mailto:ahollo@mol.hu 78 TÁLOSI, GERGÓ, AND HOLLÓ men (Penetration Grade 50/70) and polymer modified bi- tumen. 2. Experimental The aim of our study was to examine the quality im- provement opportunities of bitumen used in road con- struction (Penetration Grade 50/70). Over the course of our experiments, the modifying effect of a biocompo- nent derived from pinewood-based fast pyrolysis bio- oil on the original bitumen 50/70 (B-series) was investi- gated. Furthermore, experiments that applied the biocom- ponent together with the block copolymer styrene butadi- ene styrene (SBS) (P-series) were conducted. 2.1 Samples The investigated feedstock was a commercially-available bitumen (producer: MOL Plc.) with Penetration Grade 50/70 (softening point: 51 ◦C; penetration: 56 0.1mm). The modifying agent for the polymer modified bitumen samples was the linear block copolymer styrene buta- diene styrene (DST L 30-01, producer: Voronezhsin- tezkauchuk OJSC) applied at a concentration of 4 w/w%. The styrene concentration of the copolymer was 30 w/w%. The applied biocomponent was derived from a commercial pyrolysis liquid of pinewood (producer: BTG-BTL, The Netherlands [31], elemental analysis of the BTG-BTL product: C: 46 %, H: 7 %, O: 47 %; water content: 25 w/w%). Before blending, the sample was sta- bilised. The water content was removed and the volatile organic components that are boiling up to 220 ◦C were removed with the water as well, to adjust the flash point of the additive by distillation. Samples of the biocompo- nent were applied in concentrations of 1, 5, and 10 w/w%. Table 1 shows the composition of the samples. 2.2 Measurements In order to characterise our samples, conventional bitu- men tests were conducted before and after aging. The softening point of the samples was determined according to the standard MSZ EN 1427 using an automatic ring- and-ball softening point tester (Petrotest RKA 5). Pen- etration measurements were conducted according to the standard MSZ EN 1426 using an automatic penetrometer (Petrotest PNR 12). Aging was simulated by the rolling thin-film oven test (RTFOT) according to the standard MSZ EN 12607-1. The change in the rheological properties was investi- gated by using an Anton Paar MCR 301 type dynamic shear rheometer (DSR). The linear viscoelastic range (LVE range) of the samples was investigated by ampli- tude sweep analysis. During the measurements, the fre- quency was constant (10 s−1) and the amplitude varied. The amplitude sweep analysis was measured at 10 and 60 ◦C which resulted in the corresponding γ (deforma- tion) values. The sweep frequency response analysis was Table 1: Composition of the samples Name of the sample Bitumen 50/70, w/w% SBS, w/w% Biocomponent, w/w% Bref 100 - - Bbio01 99 - 1 Bbio05 95 - 5 Bbio10 90 - 10 Pref 96 4 - Pbio01 95 4 1 Pbio05 91 4 5 Pbio10 86 4 10 also investigated at 10 and 60 ◦C, within the frequency range of 0.01− 100 Hz. The multiple stress creep recov- ery (MSCR) test was performed at 60 ◦C over 10 cycles according to the standard ASTM D7405. During the anal- ysis, the samples were subjected to a loading force of 100 N for 1 sec which was then removed for 9 secs. After the 10th cycle, the loading force was changed to 3200 N and another 10 cycles measured. From the experimental val- ues the percentage of the recovery could be calculated. The temperature-dependence of the rheological prop- erties was investigated within the temperature range of 10− 70 ◦C, with a heating rate of 1 ◦C/min at a constant frequency of 1 Hz. During the measurements, the values of the complex modulus, namely the stiffness and com- plex viscosity, were investigated. 3. Results and Analysis 3.1 Basic properties As for the conventional properties, the softening points of the samples were increased by the biocomponent as well as SBS when compared to the original bitumen 50/70 (Table 2). The increasing effect tendentiously escalated as the amount of applied biocomponent increased. In the case of applying two modifiers simultaneously (biocom- ponent and SBS polymer), synergistic effects were expe- rienced. In accordance with the softening points, the penetra- tion of the samples decreased by using every blending component, that is the samples became harder (Table 2). The increased stiffness could be explained by the blend- ing of the rigid biocomponent. From the results measured after RTFOT, it can be stated that in the case of the SBS polymer modified bi- tumen (PmB) samples, the aging effect was less signifi- cant when compared to the sample of reference bitumen. Moreover, in the case of the PmB samples, the biocompo- nent enhanced its resistance to hardening. Every sample complied with the technical specifications for road con- struction currently in force in Hungary (for Bbio samples MSZ EN 12591, for Pbio samples MSZ EN 14023). Table 2 shows the conventional properties of the sam- ples before and after aging. Hungarian Journal of Industry and Chemistry QUALITY IMPROVEMENT OF BINDERS FOR ROAD CONSTRUCTION 79 Table 2: Basic properties of the samples and the requirements of the bitumen standards before and after RTFOT. Before RTFOT After RTFOT Softening point, ◦C Penetration, 0.1 mm Change in mass, w/w% Increase in softening point, ◦C Retained penetration, % Bref 51.0 56 −0.02 5.9 67.9 Bbio01 51.1 54 0.006 2.1 75.9 Bbio05 51.2 53 −0.30 3.9 77.4 Bbio10 51.4 52 −0.50 4.7 71.1 MSZ EN 12591 50/70 46− 54 50− 70 ≤ 0.50 ≤ 8.0 ≥ 55.0 Pref 72.5 35 0.04 2.6 82.9 Pbio01 74.1 32 0.01 3.6 96.9 Pbio05 77.0 31 −0.01 3.1 93.5 Pbio10 79.7 29 −0.02 3.2 100.0 MSZ EN 14023 10/40-65 ≥ 65.0 10− 40 ≤ 0.30 ≤ 8.0 ≥ 50.0 25/55-65 ≥ 65.0 25− 55 ≤ 0.50 ≤ 8.0 ≥ 50.0 3.2 Rheological properties The results of the amplitude sweep analysis are presented in Fig. 1. The deformation values increased by applying the SBS polymer component at 10 ◦C, thus the LVE ranges were widened. In every case, the blending of the biocom- ponent narrowed the LVE range. The addition of the rigid component resulted in a stiffer structural material with a lower resistance to deformation. The extent of the de- crease was significant in the case of the samples that were only modified with biocomponents. In the other cases, the decrease was less notable, the bituminous matrix could positively compensate for the narrowing effect of the bio- component. As for the measurements conducted at 60 ◦C, the per- manent deformation decreased in every case when com- pared to the reference sample, nevertheless, the value of the permanent deformation increased. In the case of the samples that were only modified with a biocomponent, the narrowing of the LVE ranges was minimal when com- pared to the values measured at 10 ◦C. Fig. 2 shows the results of the sweep frequency re- sponse analysis. At 10 ◦C and 60 ◦C, the values of the complex modulus were measured, presented and eval- uated at frequencies of 25 Hz and 30 Hz, respectively. These conditions correlate with similar asphalt tests. At 10 ◦C, the complex modulus values, namely the stiffness, were preferably larger in the case of every mod- ified composition when compared to the reference sample (Bref). However, in the case of the samples modified with only a biocomponent, by adding the biocomponent, the stiffness decreased tendentiously. In view of the Bbio10 sample (10 w/w%), the value of the complex modulus was almost identical to that of the reference sample. As for the PmB samples, similar behaviour was observed. By blending the biocomponent in small amounts (up to 5 w/w%), the complex modulus decreased, but in accor- dance with this effect the fatigue performance improved. The stiffness of the samples containing 10 w/w% of biocomponent (Bbio10 and Pbio10) were comparable to those of the corresponding reference samples (Bref, Pref), thus the stiffness of the asphalt mixture can be suit- able as well. The results measured at 60 ◦C showed that the com- plex modulus positively increased by the blending in of the SBS polymer. In the case of the samples modified with only a biocomponent, the modifying effect of the blending component was negligible. As for the flexible behaviour, the results of the MSCR test are presented in Figs. 3 and 4. The loaded and relax- ation periods are shown in the diagrams on the left-hand side, while the percentage of the elastic recovery at the end of the 1st and 10th cycles can be seen on the right. Every diagram shows the results measured when a load- ing force of 3200 N was applied. In the case of the samples modified with only a bio- component, the component applied in small concentra- tions (up to 5 w/w%) influenced the elastic recovery pos- itively when a loading force of 100 N was applied. How- ever, when applied at a concentration of 10 w/w%, the percentage of the elastic recovery decreased when com- pared to the reference sample, however, the value of the permanent deformation increased. When high loading forces (3200 N) were applied, the permanent deforma- tion was adversely increased by the biocomponent. The elastic recovery of the samples was negligible, the elas- tic property practically ceased under the test conditions, their viscosity dominated. In the case of the PmB samples, when a loading force of 100 N was applied, the biocomponent influenced the flexible behaviour positively at every concentration. The permanent deformation decreased and, in line with this positive effect, the percentage of the elastic recovery in- creased. Thus, the rutting performance of the samples was favourable. The higher the blending concentration of the biocomponent was, the more flexible the behaviour 47(2) pp. 77–83 (2019) 80 TÁLOSI, GERGÓ, AND HOLLÓ Figure 1: LVE range of the samples – a) 10◦C, b) 60◦C. Figure 2: Complex modulus of the samples – a) 10◦C, 25 Hz frequency, b) 60◦C, 30 Hz frequency. Figure 3: The results of the MSCR test of the samples modified with only a biocomponent when subjected to a loading force of 3200N. it exhibited. When high loading forces were applied, the tendency was the same. The favourable values of the PmB samples can be explained by a phenomenon in which the biocomponent did not influence the flexible properties of the PmB matrix, while it could mitigate the liability to deform as a result of its rigid structure. By comparing the values of every sample measured by the MSCR test, the flexible behaviour of the PmB sam- ples was significantly better than in the case of the refer- ence bitumen sample. Nowadays the modification of bitu- men is strictly necessary in order to comply with the more rigid requirements enforced by the technical practices of road construction, namely that the binders have a suf- ficient degree of rutting performance. According to our measurements, the biocomponent by itself cannot provide a good solution to this problem, however, by blending it with SBS polymer the modifying effect of the biocompo- nent could be beneficial. The results of the Temperature Sweep Analysis are presented in Fig. 5. In the case of the samples that were only modified with the biocomponent, the values of the complex modulus decreased at lower concentrations (up to 5 w/w%), but by applying the biocomponent at 10 w/w%, the values of the complex modulus were almost identical to that of the reference sample (Bref). As for the PmB samples, the stiffness increased in every case when compared to the Bref sample. The value of the complex modulus increased by increasing Hungarian Journal of Industry and Chemistry QUALITY IMPROVEMENT OF BINDERS FOR ROAD CONSTRUCTION 81 Figure 4: The results of the MSCR test of the samples modified with SBS polymer when subjected to a loading force of 3200N. the biocomponent concentration. During the Temperature Sweep Analysis, the temperature dependence of the com- plex viscosity was measured as well. The tendency was very similar to that of the temperature dependence of the complex modulus. In particular, the complex viscosity significantly increased in the case of the PmB samples when compared to the reference sample (Bref). By in- creasing the blending concentration of the biocomponent, the complex viscosity increased further. 4. Conclusion The aim of this work was to investigate the quality im- provement opportunities of bitumen used in road con- struction (Penetration Grade 50/70) and to provide a com- Figure 5: Temperature dependency of complex viscosity – a) samples modified only with biocomponent, b) PmB samples. parative analysis of the modifying effect of stabilised pinewood-based fast pyrolysis bio-oil on SBS polymer modified bitumen. According to our measurements, the biocomponent can be blended at a concentration of 10 w/w% with the tested bituminous binders. As for its flex- ible properties, the utilisation of the biocomponent by it- self does not result in a good bitumen quality, but together with the SBS polymer the rutting performance improves significantly. Based on these results, further measure- ments are being planned with other modifying agents to identify the exact and credible binder composition from economical and technical points of view. Furthermore, it is necessary to examine the methods of blending and treating opportunities of the biocomponent that can pro- duce a binder with better mechanical properties and pos- sibly with a greater concentration of the biocomponent in the blend. Acknowledgements This research was supported by the ÚNKP-18-2 New Na- tional Excellence Program of the Ministry of Human Ca- pacities. The authors thank MOL Plc. providing the bitu- men sample. REFERENCES [1] Read, J.; Whiteoak, D.: The shell bitumen handbook (Thomas Telford, London, UK, 2003) ISBN:0 7277 3220 X [2] Geiger, A.; Holló, A.; Fehér, P.: MOL Rubber bitumen–Part 1. MOL Group Scientific Magazine 2012/2, 54–61 DOI: 10.24193/subbchem.2017.2.19 [3] Holló, A.; Geiger, A.; Gergó, P.: (2017). Upgrad- ing options of heavy residues, rubber modified bi- tumen case study. Fuels and Lubricants, 1(1), 18– 23. https://www.yumpu.com/en/document/read/ 59539765/fuels-lubricants-magazine [4] Angyal, A., Miskolczi, N., Bartha, L., Gergo, P. Synthesis and evaluation of modified polyethylene wax applied as dispersant in rubber bitumen com- posites. Hung. J. Ind. Chem., 2009, 37(1), 21– 47(2) pp. 77–83 (2019) https://doi.org/10.24193/subbchem.2017.2.19 https://www.yumpu.com/en/document/read/59539765/fuels-lubricants-magazine https://www.yumpu.com/en/document/read/59539765/fuels-lubricants-magazine 82 TÁLOSI, GERGÓ, AND HOLLÓ 25 https://mk.uni-pannon.hu/hjic/index.php/ hjic/article/view/217 [5] Adorjányi, K., Füleki, P. Performance evaluation of bitumens at high temperature with multiple stress creep recovery test. Hung. J. Ind. Chem., 2011, 39(2), 195–199. https://mk.uni-pannon. hu/hjic/index.php/hjic/article/view/409 [6] Yang, X., You, Z., Dai, Q., Mills-Beale, J. Mechanical-performance of asphalt mixtures modi- fied by bio-oils derived from waste wood resources. Construction and Building Materials, 2014, 51, 424–431 DOI: 10.1016/j.conbuildmat.2013.11.017 [7] Chailleux, E., Audo, M., Bujoli, B., Queffelec, C., Legrand, J., Lepine, O. Alternative Binder from mi- croalgae: Algoroute project. Workshop: Alternative binders for sustainable asphalt pavements, 2012, pp. 7–14. [8] Fini, E. H., Kalberer, E. W., Shahbazi, A., Basti, M., You, Z., Ozer, H., Aurangzeb, Q. Chemical charac- terization of biobinder from swine manure: Sustain- able modifier for asphalt binder. Journal of Materi- als in Civil Engineering, 2011, 23(11), 1506-1513 DOI: 10.1061/(ASCE)MT.1943-5533.0000237 [9] Mills-Beale, J., You, Z., Fini, E., Zada, B., Lee, C. H., Yap, Y. K. Aging influence on rheol- ogy properties of petroleum-based asphalt mod- ified with biobinder. Journal of Materials in Civil Engineering, 2012, 26(2), 358–366 DOI: 10.1061/(ASCE)MT.1943-5533.0000712 [10] Raouf, M. A., Williams, R. C. Rheology of fraction- ated cornstover bio-oil as a pavement material. In- ternational Journal of Pavements, 2010, 9, 58–69 [11] Hill, D. R., Jennings, A. A. (2011). Bioasphalt from Urban Yard Waste Carbonization. A Student Study for the Ohio Department of Transportation Office of Research and Development (State Job Number 134464) https://rosap.ntl.bts.gov/view/dot/ 23358/dot_23358_DS1.pdf? [12] Uzun, B. B., Apaydin-Varol, E., Ateş, F., Özbay, N., Pütün, A. E. Synthetic fuel production from tea waste: characterisation of bio-oil and bio-char. Fuel, 2010, 89(1), 176–184 DOI: 10.1016/j.fuel.2009.08.040 [13] Chaiya, C. Production of bio-oil from cof- fee residue using pyrolysis process. 2010, Proceedings of the World Congress on Engi- neering and Computer Science (Vol. 2), pp. 19–21. http://www.iaeng.org/publication/ WCECS2011/WCECS2011_pp704-707.pdf [14] Onay, O., Koçkar, O. M. Pyrolysis of rape- seed in a free fall reactor for production of bio-oil. Fuel, 2006, 85(12), 1921-1928., DOI: 10.1016/j.fuel.2006.03.009 [15] Ţensöz, S., Kaynar, Ý. Bio-oil production from soy- bean (Glycine max L.); fuel properties of Bio-oil. Industrial Crops and Products, 2006, 23(1), 99–105 DOI: 10.1016/j.indcrop.2005.04.005 [16] Ingrassia, L. P., Lu, X., Ferrotti, G., Canestrari, F. Renewable materials in bituminous binders and mixtures: Speculative pretext or reliable oppor- tunity?. Resources, Conservation and Recycling, 2019, 144, 209–222 DOI: 10.1016/j.resconrec.2019.01.034 [17] Peralta, J., Williams, R. C., Rover, M., Silva, H. M. R. D. Development of rubber modified fractionated bio-oil for use as noncrude petroleum binder in flexible pavements. 2012, Transporta- tion Research Circular E-C165: Alternative Binders, pp. 23–36. https://repositorium.sdum. uminho.pt/bitstream/1822/22357/1/Artigo_ TRB-Circular-ec165-mypaper.pdf [18] Peralta, J., Williams, R. C., Silva, H. M., Machado, A. V. Combining Asphalt T-Rubber (AR) and Fast- Pyrolysis Bio-Oil to Create a Binder for Flexible Pavements. 2013, Proceedings of the 2nd Edition of the International Conference and Exhibition WASTES: Solutions, Treatments and Opportunities, pp. 1–8. https://pdfs.semanticscholar.org/ 9909/205f9c3af88eb53e053aa09dd1ec53fe73e1. pdf?_ga=2.105672513.336972856. 1576827673-1545734105.1576827673 [19] Peralta, J., Williams, R. C., Silva, H. M., Machado, A. V. A. Recombination of asphalt with bio-asphalt: Binder formulation and asphalt mixes application. Journal of the Association of Asphalt Paving Tech- nologists, 2014, 83, 1–25. [20] Grilli, A., Bocci, E., Bocci, M. Hot recycling of reclaimed asphalt using a bio-based additive. In 8th RILEM International Symposium on Testing and Characterization of Sustainable and Innovative Bituminous Materials, (Springer, Dordrecht, 2016), pp. 953–964 DOI: 10.1007/978-94-017-7342-3_76 [21] Jiménez del Barco Carrión, A., Lo Presti, D., Pouget, S., Airey, G., Chailleux, E. Linear vis- coelastic properties of high reclaimed asphalt con- tent mixes with biobinders. Road Materials and Pavement Design, 2017, 18(sup2), 241–251., DOI: 10.1080/14680629.2017.1304253 [22] Jiménez del Barco Carrión, A., Pérez-Martínez, M., Themeli, A., Lo Presti, D., Marsac, P., Pouget, S., ... Airey, G. Evaluation of bio-materials’ rejuvenating effect on binders for high-reclaimed asphalt con- tent mixtures. Materiales de Construcción, 2017, 67(327), 1–11 DOI: 10.3989/mc.2017.04516 [23] Lei, Z., Bahia, H., Yi-qiu, T., Ling, C. Effects of refined waste and bio-based oil modifiers on rheo- logical properties of asphalt binders. Construction and Building Materials, 2017, 148, 504–511 DOI: 10.1016/j.conbuildmat.2017.05.101 [24] Raouf, M. A., Williams, R. C. Temperature and shear susceptibility of a nonpetroleum binder as a pavement material. Transportation Research Record, 2010, 2180(1), 9–18 DOI: 10.3141/2180-02 [25] Cooper III, S. B., Mohammad, L. N., Elseifi, M. Evaluation of asphalt mixtures containing renew- able binder technologies. International Journal of Pavement Research and Technology, 2013, 6(5), 570–575 DOI: 10.6135/ijprt.org.tw/2013.6(5).570 Hungarian Journal of Industry and Chemistry https://mk.uni-pannon.hu/hjic/index.php/hjic/article/view/217 https://mk.uni-pannon.hu/hjic/index.php/hjic/article/view/217 https://mk.uni-pannon.hu/hjic/index.php/hjic/article/view/409 https://mk.uni-pannon.hu/hjic/index.php/hjic/article/view/409 https://doi.org/10.1016/j.conbuildmat.2013.11.017 https://doi.org/10.1061/(ASCE)MT.1943-5533.0000237 https://doi.org/10.1061/(ASCE)MT.1943-5533.0000712 https://doi.org/10.1061/(ASCE)MT.1943-5533.0000712 https://rosap.ntl.bts.gov/view/dot/23358/dot_23358_DS1.pdf? https://rosap.ntl.bts.gov/view/dot/23358/dot_23358_DS1.pdf? https://doi.org/10.1016/j.fuel.2009.08.040 http://www.iaeng.org/publication/WCECS2011/WCECS2011_pp704-707.pdf http://www.iaeng.org/publication/WCECS2011/WCECS2011_pp704-707.pdf https://doi.org/10.1016/j.fuel.2006.03.009 https://doi.org/10.1016/j.fuel.2006.03.009 https://doi.org/10.1016/j.indcrop.2005.04.005 https://doi.org/10.1016/j.resconrec.2019.01.034 https://repositorium.sdum.uminho.pt/bitstream/1822/22357/1/Artigo_TRB-Circular-ec165-mypaper.pdf https://repositorium.sdum.uminho.pt/bitstream/1822/22357/1/Artigo_TRB-Circular-ec165-mypaper.pdf https://repositorium.sdum.uminho.pt/bitstream/1822/22357/1/Artigo_TRB-Circular-ec165-mypaper.pdf https://pdfs.semanticscholar.org/9909/205f9c3af88eb53e053aa09dd1ec53fe73e1.pdf?_ga=2.105672513.336972856.1576827673-1545734105.1576827673 https://pdfs.semanticscholar.org/9909/205f9c3af88eb53e053aa09dd1ec53fe73e1.pdf?_ga=2.105672513.336972856.1576827673-1545734105.1576827673 https://pdfs.semanticscholar.org/9909/205f9c3af88eb53e053aa09dd1ec53fe73e1.pdf?_ga=2.105672513.336972856.1576827673-1545734105.1576827673 https://pdfs.semanticscholar.org/9909/205f9c3af88eb53e053aa09dd1ec53fe73e1.pdf?_ga=2.105672513.336972856.1576827673-1545734105.1576827673 https://doi.org/10.1007/978-94-017-7342-3_76 https://doi.org/10.1080/14680629.2017.1304253 https://doi.org/10.1080/14680629.2017.1304253 https://doi.org/10.3989/mc.2017.04516 https://doi.org/10.1016/j.conbuildmat.2017.05.101 https://doi.org/10.1016/j.conbuildmat.2017.05.101 https://doi.org/10.3141/2180-02 https://doi.org/10.6135/ijprt.org.tw/2013.6(5).570 QUALITY IMPROVEMENT OF BINDERS FOR ROAD CONSTRUCTION 83 [26] Ball, G. F. A., Herrington, P. A., Patrick, J. E.. Tall oil pitch as bitumen extender. New Zealand Journal of Forestry Science, 1993, 23(2), 236–242 https://www.scionresearch. com/__data/assets/pdf_file/0016/17701/ NZJFS2321993BALL236-242.pdf [27] Bearsley, S. R., Haverkamp, R. G. Age hardening potential of tall oil pitch modified bitumen. Road Materials and Pavement Design, 2007, 8(3), 467– 481 DOI: 10.1080/14680629.2007.9690084 [28] Bearsley, S. R., Haverkamp, R. G. Adhesive proper- ties of tall oil pitch modified bitumen. Road Materi- als and Pavement Design, 2007, 8(3), 449–465 DOI: 10.1080/14680629.2007.9690083 [29] Gondim, L. M., Soares, S. A., Barroso, S. H. Petroleum Plant Sap as an Asphalt Modifier for Pavement Applications. International Journal of Engineering & Technology IJCEE-IJENS, 2006, 16(06) [30] Yang, X., Mills-Beale, J., You, Z. Chemical char- acterization and oxidative aging of bio-asphalt and its compatibility with petroleum asphalt. Journal of Cleaner Production, 2017, 142, 1837–1847 DOI: 10.1016/j.jclepro.2016.11.100 [31] BTG-BTL, Pyrolysis oil properties. https://www. btg-btl.com/en/applications/oilproperties (Accessed: 11/06/2019) 47(2) pp. 77–83 (2019) https://www.scionresearch.com/__data/assets/pdf_file/0016/17701/NZJFS2321993BALL236-242.pdf https://www.scionresearch.com/__data/assets/pdf_file/0016/17701/NZJFS2321993BALL236-242.pdf https://www.scionresearch.com/__data/assets/pdf_file/0016/17701/NZJFS2321993BALL236-242.pdf https://doi.org/10.1080/14680629.2007.9690084 https://doi.org/10.1080/14680629.2007.9690083 https://doi.org/10.1080/14680629.2007.9690083 https://doi.org/10.1016/j.jclepro.2016.11.100 https://doi.org/10.1016/j.jclepro.2016.11.100 https://www.btg-btl.com/en/applications/oilproperties https://www.btg-btl.com/en/applications/oilproperties Introduction Materials and Methods Cultivation of the bacteria Analysis Neural networking Results and Discussion Conclusion Introduction Materials and Methods Cultivation of the bacteria The modelling of microbial growth Measurement of CGTase activity Results and Discussion Batch fermentation Fed-batch fermentation Semi-continuous fermentation Comparison of the different fermentation techniques Conclusion Introduction Experimental Emulsification tendency Fouling tendency Neutron activation analysis Compatibility Results and Discussion Conclusion Introduction Materials and methods Materials Measurement of enzyme activity Immobilization of laccase inside an alginate matrix Characterization of the immobilized laccase Effect of temperature on alginate beads Storage stability of the beads Reuse of the alginate beads Determination of kinetic parameters with ABTS Results and Analysis Characterization of the alginate beads Determination of kinetic parameters Effect of temperature on alginate beads Storage stability of the beads Reuse of the beads Discussion Conclusions Introduction MATERIALS AND METHODS The Collection of Samples Sample preparation Stripping of bitumen by leaching whilst being stirred Stripping of bitumen by leaching in the absence of stirring Panning of the tar-free sand Determination of the specific gravity of tar-free sand Shaking tabling Particle size analysis using a sieving method Reflected light microscopy Thin section microscopy X-ray Fluorescence spectroscopy RESULTS AND DISCUSSION Conclusions Introduction Discrete element method Calculation of the normal and tangential displacement of particles Calculation of the slip velocity Calculation of the slide of particles Calculation of the damping coefficient Treatment of particle-wall collisions Conditions of simulation Geometric representation of the silo and rotary valve Application of the GPU Simulation results and discussion The startup process of silo discharge Analysis of the mass outflow rate through the rotary valve Flow patterns of the particles at different rotational speeds of the rotor Conclusions Introduction Experimental work Tracing and Measurement Technique Modelling and Simulation Methodology Results and Analysis Experiments Possible Flow Structures and Parameters Generation of the Simulation Models Local Programming of the Generated Process Structure Evaluation of the RTD Measurements Discussion Conclusions Introduction Applied modeling method: Programmable Process Structures Modelling and analysis of Hydrogen Supply Chains by Programmable Process Structures Description of the example HSC process network Meta-prototypes-based generation and analysis of the example model Meta-prototypes-based functional prototypes of the example model Generation and simulation of the example HSC system Analysis of the simulated results Discussion Conclusion Introduction Experimental and the modelling process Results and Analysis Conclusion Introduction Materials and Methods Results and Discussion Conclusions Introduction Experimental Samples Measurements Results and Analysis Basic properties Rheological properties Conclusion