**Comparison of Di**ff**erent Hydrotalcite Solid Adsorbents on Adsorptive Desulfurization of Liquid Fuel Oil**

**Mozammel Mazumder 1 , Rajib Das 1 , Md Symon Jahan Sajib 2 , Andrew Jewel Gomes 1,**† **, Mohammad Islam <sup>3</sup> , Thinesh Selvaratnam <sup>4</sup> and Ashiqur Rahman 1, \***


Received: 28 February 2020; Accepted: 25 April 2020; Published: 27 April 2020

**Abstract:** With increasingly stringent environmental regulations, desulfurization for gasoline oil production has become an important issue. Nowadays, desulfurization technologies have become an integral part of environmental catalysis studies. It is also important for processing of fuel for fuel-cells, which has a strict requirement for sulfur content for internal combustion engines. In this study, we focused on the preparation and characterization of magnesium hydroxide/aluminum supported NiO, ZnO, ZrO2, NiO-ZnO, NiO-ZrO2, adsorbents for the adsorptive desulfurization of liquid fuels. These hydrotalcite adsorbents were prepared by co-precipitation method and used for adsorption of thiophene (in n-pentane, as model fuel) and dibenzothiophene at ambient temperature and pressure. The physicochemical behaviors of the fresh adsorbents such as structure, composition, and bonding modes were determined using X-ray diffraction (XRD), Raman spectroscopy, Fourier-transform infrared spectroscopy (FTIR), energy dispersive X-Ray analysis (EDAX), scanning electron microscopy (SEM), X-ray photoelectron spectroscopy (XPS) and thermogravimetric analysis (TGA). The sulfur concentration in the mixture (thiophene and n-pentane) was measured by UV-Vis spectrophotometry. The percentages of thiophene removal and the adsorption capacity (mg of sulfur per g of adsorbent) of the five adsorbents were compared. The adsorption performance confirmed that NiO-ZrO<sup>2</sup> and NiO-ZnO adsorbents are more efficient in removing thiophene/dibenzothiophene than that of three other adsorbents. The qualitative studies using XPS confirmed the efficient adsorption nature of modified hydrotalcite adsorbents on dibenzothiophene.

**Keywords:** adsorption; hydrotalcite; thiophene/dibenzothiophene; n-pentane; desulfurization

#### **1. Introduction**

The combustion of fuels containing sulfur produces sulfur dioxide, which is liable to a series of air pollution events. Human exposure to sulfur dioxide in the ambient air has been related to respiratory system diseases and even lung cancer [1]. Therefore, in 2006 the U.S. Environmental Protection Agency (EPA) reduced the allowable sulfur levels in liquid fuels. The gasoline sulfur limit was reduced to 30 ppm, while the diesel fuel sulfur limit was reduced to 15 ppm [2]. With time, EPA regulations

became more stringent. However, the existing desulfurization technologies failed to reduce the sulfur level of gasoline or diesel fuel to less than 10 ppm [1,3].

Various methods such as adsorption desulfurization (ADS) [3], biodesulfurization (BDS) [4], extraction desulfurization (EDS) and oxidative desulfurization (ODS) have been reported for the removal of sulfur compounds from fuels [1,5,6]. Hydrodesulfurization (HDS) is the most popular and effective sulfur removing technology in refineries. However, the production of ultra-low sulfur fuels requires a large volume of catalyst [7]. Adsorptive desulfurization can provide low sulfur fuel for fuel cells and catalyst beds. Many technologies have diverged from HDS to produce low sulfur products; however, sorption, catalytic oxidation, and evaporation show the most potential among them [8]. Oxidative desulfurization (ODS) is an alternative technique for adsorptive sulfur removal. Normally an oxidative reagent is used in combination with the catalyst to oxidize the sulfur. In recent years, several kinds of ODS systems were developed successfully such as H2O<sup>2</sup> organic matrixes, H2O2/Ti-modified zeolites, H2O2/polyoxometalates (POMs), H2O2/ionic liquid, H2O2/polyoxometalate-based ionic liquid, WO3/TiO<sup>2</sup> and CeO2/TiO<sup>2</sup> [9,10]. The oxidation products are generally removed by an extraction process using solvents. Recently Ullah et al. reported on the adsorptive removal of benzothiophene (BT) from liquid fuel using a highly porous metal-organic framework based on a bicomponent zirconium (IV) benzene-tricarboxylate Zr (BTC), and its post-synthetically modified hybrid form with dodeca-tungstophosphoric acid (HPW/Zr(BTC) [11]. Additionally, different microorganisms have been used to remove sulfur. Both aerobic and anaerobic microorganisms prove to be effective desulfurization agents while maintaining aliphatic and aromatic content in the fuel [12]. Raj et al. investigated the effect of temperature, time and mass ratio for extractive desulfurization [13].

The objective of the present study is to identify an adsorbent that selectively removes sulfur from transportation fuels. The candidate adsorbents contain nickel, zinc, and zirconia. Five different types of adsorbents have been examined. Identical support of magnesium and aluminium hydroxides (hydrotalcites) was used for each of these adsorbents in order to determine the influence of the main component (Ni, Zn and ZrO2). Nickel is relatively inexpensive and has previously shown to be more promising for sulfur removal [14]. Sulfur molecules from liquid fuels are adsorbed by direct interaction. Zinc oxide (ZnO), a highly active component used to remove sulfur species from refinery liquids [15]. Zirconium is used as an efficient adsorbent for the desulfurization because of its moderate surface area, bifunctional properties of acid and base [14]. The addition of nickel on ZnO and ZrO<sup>2</sup> is an innovative approach that takes advantage of the selectivity of Ni towards S-species and the high adsorptive capacity of ZnO and ZrO<sup>2</sup> support. It changes the nature of the active metal sites and increases the sulfur removing capacity of the adsorbent. Therefore, the current study compares different types of adsorbents in order to pave the way to remove sulfur to ultra-low levels by the selective adsorption of sulfur from liquid fuels.

#### **2. Materials and Methods**

Five modified hydrotalcite adsorbents were prepared by the combination of Mg, Al, Ni, Zn, and ZrO2. The first adsorbent is the combination of Mg: Al: Ni in the approximate molar ratio of 4:1:3. It is a highly active nickel adsorbent supported on magnesium and aluminium. Increasing the capacity of a nickel adsorbent is highly beneficial, as nickel adsorbents are very effective at removing sulfur compounds. It is believed that the addition of some substance to the fuel can help increase the capacity of the nickel adsorbent. The second adsorbent is the combination of Mg: Al: Zn in the approximate molar ratio of 4:1:3. The third adsorbent is prepared with the combination of Mg: Al: Ni: Zn in the approximate molar ratio of 4:1:3:3. The fourth one is Mg: Al: ZrO<sup>2</sup> in the molar ratio of 4:1:7. The fifth adsorbent is Mg: Al: ZrO2: Ni in the mole ratio 4:1:7:3. All the adsorbents are prepared by the co-precipitation method. They are used for adsorption of thiophene (in n-pentane) and dibenzothiophene at ambient temperature and pressure. A model fuel consisting of 50 ml of n-pentane (C5H12) and thiophene (100 ppm) was prepared, and adsorption experiments were performed under ambient conditions. The physicochemical characterizations of the fresh adsorbents

such as structure, component, bond, elements, and composition were conducted using X-ray diffraction (XRD), Raman spectroscopy, Fourier-transform infrared spectroscopy (FTIR), energy dispersive X-Ray analysis (EDAX), scanning electron microscopy (SEM), X-ray photoelectron spectroscopy (XPS) and thermogravimetric analysis (TGA). The sulfur concentration in the mixture was monitored by UV-Vis spectrophotometry. Finally, the thiophene removal efficiency of the five adsorbents were compared. In addition, qualitative studies using XPS were performed to investigate the sulfur removal efficiency and adsorptive nature of adsorbents on dibenzothiophene.

#### *2.1. Adsorbent Preparation*

Five hydrotalcites with Mg: Al ratio 3:1 were prepared by the co-precipitation method. In this method, two solutions, A and B, were added at the same rate 50 ml h−<sup>1</sup> to a beaker containing 100 ml of deionized water while stirring.

Solution A was prepared by mixing an equimolar solution of Mg and Al metal nitrates (200 ml) in the 3:1 molar ratio [16]. For making the adsorbent, Mg and Al metal nitrates were used as supporting metals and Zn, Ni, Zr were used as active metal. To make Ni adsorbent, Mg: Al: Ni molar ratio were used as 4:1:3. For making Zn adsorbent, Mg: Al: Zn molar ratio were used as 4:1:3. For making Zn and Ni combined adsorbent, Mg: Al: Ni: Zn molar ratio were used as 4:1:3:3. Zr adsorbent was prepared by maintaining Mg: Al: Zr molar ratio 4:1:7. Lastly, Ni and Zr adsorbent were prepared by maintaining Mg: Al: Ni: Zr molar ratio 4:1:3:7.

Solution B was prepared by dissolving 14 g sodium hydroxide (0.35 mol) and 15.9 g sodium carbonate (0.15 mol) in 200 ml deionized water. The pH of the suspensions was around 10. The precipitates were aged at 75 ◦C for 18 h in a dryer. The resulting product was filtered, washed thoroughly with deionized water until the filtrate showed no presence of NaOH, and subsequently dried at 95 ◦C for 24 hours. Part of the samples was heated at 450 ◦C for 12 h in a furnace for further calcinations and catalytic activity study.

#### *2.2. Model Fuels using Thiophene (in Pentane) and Dibenzothiophene*

Model fuels were used in some of the adsorbent testing to determine the selectivity of the adsorbent towards certain compounds. The model fuel consisted of pentane (C5H12) and thiophene, an organosulfur compound with the chemical formula C4H4S. The composition was 50 ml of pentane and 5 microliters of thiophene. It contained 100 ppm of thiophene in 50 ml of pentane. The sulfur compounds found in the model fuel contained highly substituted sulfur compounds, which tended to be more difficult to remove due to the steric hindrance around the sulfur atom. For making the calibration curve, five different concentrations of thiophene and pentane solution (20 ppm, 15 ppm, 10 ppm, 5 ppm, 0 ppm) were prepared. In addition, for the adsorbent performance testing, five different concentrations (100 ppm, 50 ppm, 25 ppm, 12 ppm, 6 ppm of thiophene in pentane solution) were prepared. A five-point calibration curve was made in order to retain the quality of analysis. For more accuracy and to understand the error, the response at each concentration was repeated to obtain the error bar from the responses. Additionally, the adsorbents were mixed with the dibenzothiophene for a qualitative study using the XPS. Percentages of sulfur removal capacity by the five adsorbents were measured by these five different concentrations (100 ppm, 50 ppm, 25 ppm, 12 ppm, 6 ppm of thiophene (in pentane solution) and dibenzothiophene.

#### *2.3. Characterization Techniques*

The powder X-ray diffraction (XRD) experiment of the sample was performed by the Bruker AXS D8 Discover diffractometer with GADDS (General Area Detector Diffraction System, Bruker Corporation, Billerica, MA, USA); operated by a Cu-Kα radiation source and filtered with a graphite monochromator (l = 1.5406 A◦ ). A HI STAR two-dimensional multi-wire area detector was used. All the samples were first ground to make very fine particle powder. The X-ray beam was 40 kV and 40 mA power. The incident ω angle was 5◦ . A laser system was used to ensure the alignment of the sample position

on the instrument. XRD scans were recorded from 5◦ to 77◦ for 2θ with a 0.050◦ step-width and 180 step times. Further, the crystal sizes of different adsorbents were determined by Scherrer equation (L = kλ/(FWHM)Cosθ) with dimensionless shape factor, k = 0.94 and x-ray wavelength, λ = 1.5406 A◦ . The full width at half maximum (FWHM) was determined using OriginPro software (version 9.0, Northampton, MA, USA). The multiple peak fit function was used in the case of nonlinear curve fitting.

Fourier-transform infrared spectroscopy (FTIR) spectra were recorded on a BRUKER TENSOR 27 FTIR Spectrometer (Bruker Corporation, Billerica, MA, USA). This machine works in the range of 400–4000 cm−<sup>1</sup> wave number. In order to minimize the amount of bound water, the samples were kept in an air-tight container at room temperature until measurement, however, the possibility of water absorption from the atmosphere is not entirely excluded.

Hitachi S-3400N scanning electron microscope (SEM), (Hitachi, Ltd., Tokyo, Japan) was used to investigate the morphology and the composition of the adsorbents. Automatic beam axis alignment functions like auto beam setting and auto axial alignment were used. For this experiment, high SE resolution of 10 nm at 3 KV and 5-axis motorized stage with high tilt (−20 to +90 degree) and allowance for tall samples up to 80 mm high were used.

The sample's thermal stability was studied by the thermogravimetric analysis (TGA) by the TGA instrument (Netzsch, STA 449C Jupiter, and TA Instruments SDT Q-600, Erich NETZSCH GmbH & Co. Holding KG, Selb, Germany). All the samples were heated from 25 to 1200 ◦C, and an airflow rate of 60 mL min−<sup>1</sup> was maintained. The heating rate was 20 ◦C min−<sup>1</sup> . Differential scanning calorimetry (DSC) (Netzsch, STA 449C, and TA Instruments SDT Q-600, Erich NETZSCH GmbH & Co. Holding KG, Selb, Germany) measurements were carried out by maintaining the nitrogen flow rate of about 60 mL min−<sup>1</sup> .

The Raman spectra of the sample were measured using Perkin Elmer, RAMAN FLEX 400 Raman spectrometer (Perkin Elmer, Waltham, MA, USA). This machine probes in the spectral range from 230 cm−<sup>1</sup> to 3,500 cm−<sup>1</sup> Raman shift.

The UV visible measurements were performed on a Varian Cary 50 Version 3 UV Visible Spectrometer (Agilent Technologies, Santa Clara, CA, USA) coupled with the potentiostat for applying electrochemical potentials. To study the evolution, the difference between the maximum peak absorbance at a particular wavelength (λmax), and the absorbance at the initial scanning wavelength (λ0) was accounted and plotted vs. time [17].

#### **3. Result and Discussions**

From the UV spectroscopy, pure pentane absorbance value of 0.564 at the 230 nm wavelength was determined. For every 20 ppm, 15 ppm, 10 ppm, 5 ppm concentration of thiophene in 50 ml pentane was measured at the 230 nm wavelength and then pure pentane value was subtracted from the measured value to determine the thiophene concentration. The absolute values for different thiophene concentrations are shown in Figure 1. The calibration curve and equation are used to calculate the thiophene concentration in pentane after using the adsorbent.

*Technologies* **2020**, *8*, x FOR PEER REVIEW 5 of 16

**Figure 1.** Calibration curve for thiophene concentration in pentane. **Figure 1.** Calibration curve for thiophene concentration in pentane.

#### *3.1. XRD Analysis of the Adsorbents 3.1. XRD Analysis of the Adsorbents*

Figure 2a–e shows the X-ray diffractogram of the five adsorbents. From Figure 2a, it is clear that there has been the presence of hydrotalcite in the Ni adsorbent. Hydrotalcite peaks were found at 11.35°, 22.4°, 34.34°, 38.49°, 61.72°, 65.41° 2-theta angles. Figure 2a shows the presence of aluminium oxide, magnesium oxide, nickel oxide, magnesium aluminates, nickel aluminates, and magnesium nickel. Moreover, the formation of the hydrotalcite in the samples indicates the presence of carbonate in the interlayer of hydrotalcite. Hydrotalcite peaks are also observed at the same 2-theta angles for the Zn adsorbent (Figure 2b) indicating the presence of aluminium oxide, magnesium oxide, zinc oxide, magnesium aluminates, zinc aluminates, and magnesium zinc etc. Hydrotalcite peaks were found at 22.4°, 34.34°, 61.72°, 71.6° 2-theta angle (JCPDS file no. 14-0191) for the Ni and Zn adsorbent (Figure 2c). The XRD patterns of Zr, and Zr and Ni (Figure 2d,e) indicate the presence of zirconia majorly in the form of the metastable tetragonal phase and minorly in the form of the monoclinic phase [18]. Generally, the tetragonal phase of zirconia can be stabilized by incorporating the promoters into the zirconia lattice. It is well known that the tetragonal phase of Zr is more active in catalysis [14]. The average crystalline sizes were calculated as 9.1 ± 2.8 nm for Ni adsorbent, 14.0 ± 3.9 nm for Zn adsorbent, 14.9 ± 7.6 nm for Ni and Zn adsorbent, 12.5 ± 1.0 nm for Zr adsorbent, 16.0 ± 7.5 nm for Zr and Ni adsorbent. The peak centers and corresponding FWHM are presented in the Supplementary Materials (Table S1). Figure 2a–e shows the X-ray diffractogram of the five adsorbents. From Figure 2a, it is clear that there has been the presence of hydrotalcite in the Ni adsorbent. Hydrotalcite peaks were found at 11.35◦ , 22.4◦ , 34.34◦ , 38.49◦ , 61.72◦ , 65.41◦ 2-theta angles. Figure 2a shows the presence of aluminium oxide, magnesium oxide, nickel oxide, magnesium aluminates, nickel aluminates, and magnesium nickel. Moreover, the formation of the hydrotalcite in the samples indicates the presence of carbonate in the interlayer of hydrotalcite. Hydrotalcite peaks are also observed at the same 2-theta angles for the Zn adsorbent (Figure 2b) indicating the presence of aluminium oxide, magnesium oxide, zinc oxide, magnesium aluminates, zinc aluminates, and magnesium zinc etc. Hydrotalcite peaks were found at 22.4◦ , 34.34◦ , 61.72◦ , 71.6◦ 2-theta angle (JCPDS file no. 14-0191) for the Ni and Zn adsorbent (Figure 2c). The XRD patterns of Zr, and Zr and Ni (Figure 2d,e) indicate the presence of zirconia majorly in the form of the metastable tetragonal phase and minorly in the form of the monoclinic phase [18]. Generally, the tetragonal phase of zirconia can be stabilized by incorporating the promoters into the zirconia lattice. It is well known that the tetragonal phase of Zr is more active in catalysis [14]. The average crystalline sizes were calculated as 9.1 ± 2.8 nm for Ni adsorbent, 14.0 ± 3.9 nm for Zn adsorbent, 14.9 ± 7.6 nm for Ni and Zn adsorbent, 12.5 ± 1.0 nm for Zr adsorbent, 16.0 ± 7.5 nm for Zr and Ni adsorbent. The peak centers and corresponding FWHM are presented in the Supplementary Materials (Table S1).

Overall, the XRD patterns show both sharp peaks and broad humped peaks, as are presented by the Figure 2a–e, which indicate that the adsorbent materials are partially crystalline and partially amorphous. Overall, the XRD patterns show both sharp peaks and broad humped peaks, as are presented by the Figure 2a–e, which indicate that the adsorbent materials are partially crystalline and partially amorphous.

**Figure 2.** XRD analysis of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent (legends indicated on top of the peaks correspond to the various crystalline and amorphous phases identified from JCPDS file no. 14-0191). **Figure 2.** XRD analysis of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent (legends indicated on top of the peaks correspond to the various crystalline and amorphous phases identified from JCPDS file no. 14-0191).

#### *3.2. SEM Analysis of the Adsorbents*

*3.2. SEM Analysis of the Adsorbents*  The SEM images of the adsorbents are shown in Figure 3a–e. These figures show the surface texture of the adsorbents. The micrographs reveal that the materials have definite crystalline structures and suggest that no deformation of the metal oxides occurred during the preparation. Hydrotalcite is frequently reported as the substituted form of brucite [Mg(OH)2] with the related hexagonal crystal shape. For instance, Brady et al. reported a hexagonal/rhombohedral crystal shape for the hydrotalcite [19]. In addition, our support materials; magnesium and aluminium interlocked with the active materials such as nickel and zinc that lead to a definite crystal shape and an increase in the adsorbent volume. The elemental composition of the five calcined adsorbents from EDAX analyses are presented in wt% in Figure 4. This figure clearly shows the relative abundance of elements in the adsorbent. Additionally, the elements and their corresponding atomic percentage (at. The SEM images of the adsorbents are shown in Figure 3a–e. These figures show the surface texture of the adsorbents. The micrographs reveal that the materials have definite crystalline structures and suggest that no deformation of the metal oxides occurred during the preparation. Hydrotalcite is frequently reported as the substituted form of brucite [Mg(OH)2] with the related hexagonal crystal shape. For instance, Brady et al. reported a hexagonal/rhombohedral crystal shape for the hydrotalcite [19]. In addition, our support materials; magnesium and aluminium interlocked with the active materials such as nickel and zinc that lead to a definite crystal shape and an increase in the adsorbent volume. The elemental composition of the five calcined adsorbents from EDAX analyses are presented in wt% in Figure 4. This figure clearly shows the relative abundance of elements in the adsorbent. Additionally, the elements and their corresponding atomic percentage (at. %) are presented in the Supplementary Materials (Table S2).

%) are presented in the Supplementary Materials (Table S2).

*Technologies* **2020**, *8*, x FOR PEER REVIEW 7 of 16

*Technologies* **2020**, *8*, x FOR PEER REVIEW 7 of 16

**Figure 3.** SEM analysis of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent. **Figure 3.** SEM analysis of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent. **Figure 3.** SEM analysis of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent.

**Figure 4.** Elemental composition of the five calcined adsorbents from energy dispersive X-Ray analysis (EDAX). **Figure 4.** Elemental composition of the five calcined adsorbents from energy dispersive X-Ray analysis (EDAX). **Figure 4.** Elemental composition of the five calcined adsorbents from energy dispersive X-Ray analysis (EDAX).

#### *3.3. FTIR Analysis of the Adsorbents 3.3. FTIR Analysis of the Adsorbents 3.3. FTIR Analysis of the Adsorbents*

The recorded peaks corresponding to various functional groups of the carbonate intercalated hydrotalcite are presented in Table 1. The recorded spectra are presented by the Figure S1 in the Supplementary Materials. The recorded peaks corresponding to various functional groups of the carbonate intercalated hydrotalcite are presented in Table 1. The recorded spectra are presented by the Figure S1 in the Supplementary Materials. The recorded peaks corresponding to various functional groups of the carbonate intercalated hydrotalcite are presented in Table 1. The recorded spectra are presented by the Figure S1 in the Supplementary Materials.

When compared to the infrared spectrum of the sample and hydrotalcite reference materials, a number of similarities, as well as some differences, are observed. For the Ni adsorbent, the metalhydroxide peaks present in the hydroxyl stretching region (3000 cm−1–4000 cm−1) are similar to the bands found in the hydrotalcite materials. The formation of a number of C-O stretches that either correlates with carbonate or oxalate appears at 1700 cm−1 and 1400 cm−1. The expected peaks for carbonate in hydrotalcite are 1640 cm−1, 1365 cm−1, and 1313 cm−1. However, it is difficult to accurately determine the bands of carbonate or oxalate in this range because of the similarities of carbonate and oxalate stretches due to the carbon-oxygen stretching vibrations. The Zn adsorbent shows the metal-oxide and metal-hydroxide vibration between 3000 cm−1 to When compared to the infrared spectrum of the sample and hydrotalcite reference materials, a number of similarities, as well as some differences, are observed. For the Ni adsorbent, the metalhydroxide peaks present in the hydroxyl stretching region (3000 cm−1–4000 cm−1) are similar to the bands found in the hydrotalcite materials. The formation of a number of C-O stretches that either correlates with carbonate or oxalate appears at 1700 cm−1 and 1400 cm−1. The expected peaks for carbonate in hydrotalcite are 1640 cm−1, 1365 cm−1, and 1313 cm−1. However, it is difficult to accurately determine the bands of carbonate or oxalate in this range because of the similarities of carbonate and When compared to the infrared spectrum of the sample and hydrotalcite reference materials, a number of similarities, as well as some differences, are observed. For the Ni adsorbent, the metal-hydroxide peaks present in the hydroxyl stretching region (3000 cm−1–4000 cm−<sup>1</sup> ) are similar to the bands found in the hydrotalcite materials. The formation of a number of C-O stretches that either correlates with carbonate or oxalate appears at 1700 cm−<sup>1</sup> and 1400 cm−<sup>1</sup> . The expected peaks for carbonate in hydrotalcite are 1640 cm−<sup>1</sup> , 1365 cm−<sup>1</sup> , and 1313 cm−<sup>1</sup> . However, it is difficult to accurately determine the bands of carbonate or oxalate in this range because of the similarities of carbonate and oxalate stretches due to the carbon-oxygen stretching vibrations.

3500 cm−1 wavenumbers [20]. There is a close correlation between the sample and reported peaks in oxalate stretches due to the carbon-oxygen stretching vibrations. The Zn adsorbent shows the metal-oxide and metal-hydroxide vibration between 3000 cm−1 to 3500 cm−1 wavenumbers [20]. There is a close correlation between the sample and reported peaks in The Zn adsorbent shows the metal-oxide and metal-hydroxide vibration between 3000 cm−<sup>1</sup> to 3500 cm−<sup>1</sup> wavenumbers [20]. There is a close correlation between the sample and reported peaks in the

1350 cm−1–1380 cm−<sup>1</sup> range which ascribed to the carbonate antisymmetric increase, which indicates the hydrotalcite has carbonate between the layers.

At the Ni-Zn adsorbent the presence of water stretching bands at 1700 cm−<sup>1</sup> is clear. This presence is mostly due to the existence of a water bridging mode at 3401 cm−<sup>1</sup> in the hydroxyl stretching region. The presence of a number of C-O stretches that also associate with carbonate or oxalate comes out at 1700 cm−<sup>1</sup> and 1400 cm−<sup>1</sup> . The probable peaks for carbonate in hydrotalcite are 1640 cm−<sup>1</sup> , 1365 cm−<sup>1</sup> , and 1313 cm−<sup>1</sup> .

Zirconium adsorbent corresponds to the stretching vibration νOH at the broad bands ranging between 3000 and 3400 cm−<sup>1</sup> . This indicates the presence of both free and hydrogen-bonded OH groups on the sample. The peak in the region at 1620 cm−<sup>1</sup> is given to the δHOH of quasi-free H3O<sup>+</sup> group. The presence of the sharp band at 1620 cm−<sup>1</sup> is due to the residual presence at the surface and is typically zirconia and water connected. Bands around 700–500 cm−<sup>1</sup> correspond to Zr–O2–Zr asymmetric and the formation of ZrO<sup>2</sup> phases is confirmed by the Zr–O stretching modes.


**Table 1.** Infrared peak of the five adsorbents.

#### *3.4. Raman Analysis of the Five Adsorbents*

Table 2 presents the peak list of the Raman spectroscopy carried out on the adsorbent samples (the recorded spectra are presented by Figure S2 in the Supplementary Materials). The spectra of carbonate intercalate hydrotalcite with the hydrotalcite show the correlation between the modified hydrotalcite and the reference hydrotalcite. In addition, the spectra show clear proof of the presence of carbonate in the adsorbents. This corroborates with the absence of the peaks that are associated with the free sodium carbonate as presented by powder XRD pattern in Figure 2 (cf. vide supra). These results further confirm that the carbonate is successfully intercalated into the interlayer spaces of the hydrotalcite.

For the Ni adsorbent in the hydroxyl stretching region (3000–4000 cm−<sup>1</sup> ) there is some variation between the carbonate intercalated hydrotalcite and the reference hydrotalcite. The first of these variations are in the metal hydroxide bands, which shifted to lower wavenumbers, which is evocative to a variation in the chemical environment because of the presence of carbonate [22,23].

For the Zn adsorbent, there is some peak at higher wavenumbers for the metal hydroxide stretching vibrations which take place at 3389 cm−<sup>1</sup> and 3200 cm−<sup>1</sup> as reported previously [24]. In addition, there is a peak for the interlayer water/water-carbonate bridging mode at 3239 cm−<sup>1</sup> . For this atom, this is the indication of some control of freedom due to hydrogen bonding (which is seen around 200 cm−<sup>1</sup> ). The aluminium hydroxide deformation at 999 cm−<sup>1</sup> is at the same wavenumber as reported

in the literature. This indicates that the metal-oxygen and metal hydroxide bonds form in a similar way as those in the other hydrotalcites reported before [24].


**Table 2.** Raman peak list of the five adsorbents.

For the Ni-Zn adsorbent some differences are observed between the carbonate intercalated modified hydrotalcite and the hydrotalcite in the hydroxyl stretching region (3000–3500 cm−<sup>1</sup> ). The first difference is in the metal hydroxide bands, which shifted to lower wavenumbers. It indicates the difference in the chemical environment, due to the presence of carbonate. This twist is caused by the presence of the carbonate anions in the structure of the material. There is a band present at 950 cm−<sup>1</sup> , which is a carbon-oxygen bond in the carbonate anions. The expected carbonate peak appeared at 1750 cm−<sup>1</sup> instead of 1660 cm−<sup>1</sup> . This difference in peak position is due to the presence of carbonate anions, which is opposing to the carbonate anions within the interlayer. Furthermore, there is a peak visible at 200 cm−<sup>1</sup> , which disperses as a metal-oxygen band [24], and there is also a band at 100 cm−<sup>1</sup> , which is likely a hydrogen bonding band [24]. The presence of these peaks suggests that there is a little dissimilarity between the metal cation layers of the carbonate and modified hydrotalcites.

For the Zr adsorbent, there is some peak at higher wavenumbers for the metal hydroxide stretching vibrations, which take place at 3289 cm−<sup>1</sup> and 3239 cm−<sup>1</sup> as reported in the literature [24]. Furthermore, there is also a peak for the interlayer water/water-carbonate bridging mode at 3239 cm−<sup>1</sup> . For this atom, this indicates some control of freedom for these atoms due to hydrogen bonding (which is observed around 100 cm−<sup>1</sup> ). This is an indication of the metal-oxygen and metal hydroxide bonds forming in a similar way like the other reported hydrotalcites [24].

#### *3.5. TGA and DSC Analyses of the Five Adsorbents*

The thermogravimetric study of the Ni adsorbent shows a number of features (Figure 5a–e). The weight loss has been observed to take place in four steps; loosely bound water, tightly bound water, dehydroxylation, and decarbonization. The first weight loss is a mass loss step around 110 ◦C, which is due to the loss of loosely bound surface water. The large proportion of the mass loss indicates that this sample also has a large portion of tightly bound water molecules when analyzed, though it is calcined to 450 ◦C. The 6.73% mass loss around 400 ◦C is due to the loss of water molecules that are most tightly bound to the hydrotalcite and in the interlayer spaces. There are a number of simultaneous dehydroxylation steps occurring around 400 ◦C. At around 750 ◦C, weight loss is due to the loss of carbon monoxide and water from the sample, indicating simultaneous dehydroxylation and decarbonization. The mass loss at 1050 ◦C indicates the completion of the decarbonization step. The percentage mass loss and the calculated decomposition steps are shown in Table 3.


950 0.57 5MgO + MgAl2O<sup>4</sup>

**Table 3.** Mass loss step for calcined five adsorbents. Ni supported 110 8.03 Mg6Al2(OH)16(CO3)·2H2O

**Mass Loss (wt%)** 

**Proposed Formula** 

*Technologies* **2020**, *8*, x FOR PEER REVIEW 10 of 16

simultaneous dehydroxylation steps occurring around 400 °C. At around 750 °C, weight loss is due to the loss of carbon monoxide and water from the sample, indicating simultaneous dehydroxylation and decarbonization. The mass loss at 1050 °C indicates the completion of the decarbonization step.

**Table 3.** Mass loss step for calcined five adsorbents.

The percentage mass loss and the calculated decomposition steps are shown in Table 3.

**(° C )** 

**Figure 5.** *Cont*.

*Technologies* **2020**, *8*, x FOR PEER REVIEW 11 of 16

**Figure 5.** *Cont*.

*Technologies* **2020**, *8*, x FOR PEER REVIEW 12 of 16

**Figure 5.** TGA of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent. **Figure 5.** TGA of calcined (**a**) Ni adsorbent; (**b**) Zn adsorbent (**c**) Ni and Zn adsorbent, (**d**) Zr adsorbent (**e**) Zr and Ni adsorbent.

For the Zn adsorbent, from the TGA data in figure 5b, the decomposition steps are the loss of some interlayer water up to 200 °C. The next decomposition step takes place at 280 °C. The first of these steps is the loss of the last interlayer water by the relocation of the brucite-like layer. This is followed by the change in the structure of the material by the partial dehydroxylation of the brucite layer and the partial merging of the carbonate group into the structure. The third decomposition step around 400 °C signals the end of the dehydroxylation, followed by the finishing point of decarbonization, which finishes around 750 °C [20]. The final dehydroxylation of the products to form a mixed metal solid solution occurs at almost 950 °C. For the Zn adsorbent, from the TGA data in Figure 5b, the decomposition steps are the loss of some interlayer water up to 200 ◦C. The next decomposition step takes place at 280 ◦C. The first of these steps is the loss of the last interlayer water by the relocation of the brucite-like layer. This is followed by the change in the structure of the material by the partial dehydroxylation of the brucite layer and the partial merging of the carbonate group into the structure. The third decomposition step around 400 ◦C signals the end of the dehydroxylation, followed by the finishing point of decarbonization, which finishes around 750 ◦C [20]. The final dehydroxylation of the products to form a mixed metal solid solution occurs at almost 950 ◦C.

For the Ni and Zn adsorbent (Figure 5c), the 8.55 % mass loss around 400 °C corresponds to the loss of water molecules that are most tightly bound to the hydrotalcite and formed in the interlayer spaces. The weight loss around 750 °C corresponds to the loss of carbon monoxide and water from the sample and indicates a simultaneous dehydroxylation and decarbonization. The mass loss finishes at 950 °C, which completes the decarbonization step. For the Ni and Zn adsorbent (Figure 5c), the 8.55 % mass loss around 400 ◦C corresponds to the loss of water molecules that are most tightly bound to the hydrotalcite and formed in the interlayer spaces. The weight loss around 750 ◦C corresponds to the loss of carbon monoxide and water from the sample and indicates a simultaneous dehydroxylation and decarbonization. The mass loss finishes at 950 ◦C, which completes the decarbonization step.

At the Zr adsorbent (Figure 5d) the next decomposition step takes place at 600 °C. The first of these is the loss of the last interlayer water through the relocation of the brucite-like layer. This is followed by the change in the structure of the material by the partial dehydroxylation of the brucite layer and the partial merging of the carbonate group into the structure. The third decomposition step occurs around 600 °C, which is the end of the dehydroxylation followed by the finishing point of decarbonization, which finishes around 700 °C [20]. The final dehydroxylation of the products to form a mixed metal solid solution occurs at almost 900 °C. At the Zr adsorbent (Figure 5d) the next decomposition step takes place at 600 ◦C. The first of these is the loss of the last interlayer water through the relocation of the brucite-like layer. This is followed by the change in the structure of the material by the partial dehydroxylation of the brucite layer and the partial merging of the carbonate group into the structure. The third decomposition step occurs around 600 ◦C, which is the end of the dehydroxylation followed by the finishing point of decarbonization, which finishes around 700 ◦C [20]. The final dehydroxylation of the products to form a mixed metal solid solution occurs at almost 900 ◦C.

According to the mass losses recorded from the TGA for the Ni and Zr adsorbent, it is clear that there are major mass losses occurring in the sample. These are for the loosely bound water, tightly bound water, dehydroxylation, and decarbonization. The figure also shows an exothermic reaction, and the sample is gradually degraded. The sample is heated to 1200 °C, whereas the major decomposition occurs up to 600 °C. According to the mass losses recorded from the TGA for the Ni and Zr adsorbent, it is clear that there are major mass losses occurring in the sample. These are for the loosely bound water, tightly bound water, dehydroxylation, and decarbonization. The figure also shows an exothermic reaction, and the sample is gradually degraded. The sample is heated to 1200 ◦C, whereas the major decomposition occurs up to 600 ◦C.

#### *3.6. XPS Analysis of the Three Adsorbents*

The three adsorbents were mixed with the dibenzothiophene. Then the samples were tested by the XPS to check the bonding between the sulfur molecule and the adsorbents. Figure 6 below shows the clear evidence of bonding between sulfur molecule and the adsorbents. *3.6. XPS Analysis of the Three Adsorbents*  The three adsorbents were mixed with the dibenzothiophene. Then the samples were tested by the XPS to check the bonding between the sulfur molecule and the adsorbents. Figure 6 below shows

the clear evidence of bonding between sulfur molecule and the adsorbents.

*Technologies* **2020**, *8*, x FOR PEER REVIEW 13 of 16

**Table 4.** Relative percent abundances of species from XPS.  **Zinc Aluminium Magnesium Nickel Sulfur**  Ni Adsorbent - 34.80 44.17 19.47 1.57 Table 4 shows percentages of relative abundances of species from XPS. It shows clear evidence that sulfur is observed on substrates immersed in dibenzothiophene, irrespective of the formulation of the hydrotalcites preparations. Substrates with nickel showed increased sulfur content over that without nickel: Mg: Al: Ni: ≡ 1.57%, Mg: Al: Ni: Zn ≡ 1.48%, Mg: Al: Zn: ≡ 0.87%.

Ni and Zn Adsorbent 6.19 23.97 53.53 14.83 1.48 **Table 4.** Relative percent abundances of species from XPS.

Zn Adsorbent 6.87 37.25 55.00 - 0.87


without nickel: Mg: Al: Ni: ≡ 1.57%, Mg: Al: Ni: Zn ≡ 1.48%, Mg: Al: Zn: ≡ 0.87%

#### *3.7. Performance Comparison of All Five Types of Adsorbents*

*3.7. Performance Comparison of All Five Types of Adsorbents*  Five different types of adsorbents were examined. Each of these adsorbents has the same support Five different types of adsorbents were examined. Each of these adsorbents has the same support of magnesium and aluminium; therefore, an influence specific to the main component is determined.

of magnesium and aluminium; therefore, an influence specific to the main component is determined. Figure 7 shows that, after an hour, fresh thiophene reaches 71 ppm, and when hydrotalcite without any metal is used. The zirconia adsorbent reacts slowly in comparison to the other four adsorbents. On the other hand, nickel-based adsorbents react more quickly compared to the other three adsorbents. As shown in Figure 7, all five adsorbents work within an hour until the thiophene concentration in the pentane becomes almost constant. It turns out that the deposition of nickel on zirconia is an innovative approach that takes advantage of the selectivity of Ni towards S-species and the high adsorptive capacity of zirconia support. It is observed that adding zirconia into the adsorbent structure increase the sulfur removal capacity. Possibly, the added zirconia changes the nature of the Figure 7 shows that, after an hour, fresh thiophene reaches 71 ppm, and when hydrotalcite without any metal is used. The zirconia adsorbent reacts slowly in comparison to the other four adsorbents. On the other hand, nickel-based adsorbents react more quickly compared to the other three adsorbents. As shown in Figure 7, all five adsorbents work within an hour until the thiophene concentration in the pentane becomes almost constant. It turns out that the deposition of nickel on zirconia is an innovative approach that takes advantage of the selectivity of Ni towards S-species and the high adsorptive capacity of zirconia support. It is observed that adding zirconia into the adsorbent structure increase the sulfur removal capacity. Possibly, the added zirconia changes the nature of the active metal sites

active metal sites into the structure and increases its sulfur removal capacity. Figure 8 shows the

into the structure and increases its sulfur removal capacity. Figure 8 shows the percent removal of thiophene with respect to different time interval. It is evident that for all the adsorbents the maximum percentage of the removal occurred in the 0–10 min interval. It is also worth noting that more than 70% removal for both NiO-ZrO2, and NiO-ZnO was obtained in the first 10 min. This finding can be used to investigate the reaction kinetics in future studies. *Technologies* **2020**, *8*, x FOR PEER REVIEW 14 of 16 adsorbents the maximum percentage of the removal occurred in the 0–10 min interval. It is also worth noting that more than 70% removal for both NiO-ZrO2, and NiO-ZnO was obtained in the first 10 min. This finding can be used to investigate the reaction kinetics in future studies. *Technologies* **2020**, *8*, x FOR PEER REVIEW 14 of 16 adsorbents the maximum percentage of the removal occurred in the 0–10 min interval. It is also worth noting that more than 70% removal for both NiO-ZrO2, and NiO-ZnO was obtained in the first 10

**Figure 7.** Thiophene concentrations in pentane with time. **Figure 7.** Thiophene concentrations in pentane with time. **Figure 7.** Thiophene concentrations in pentane with time.

#### **4. Conclusions 4. Conclusions**

**Figure 8.** Percent removal of thiophene with respect to different time interval. **4. Conclusions**  In the current study, we prepared hydrotalcite by the co-precipitation method and characterized using various analytical techniques. The analysis matches the typical behavior of hydrotalcites of the same composition reported in the literature. From the analysis of the modified hydrotalcites, it can be concluded that carbonate was successfully intercalated into the structure of hydrotalcite using co-In the current study, we prepared hydrotalcite by the co-precipitation method and characterized using various analytical techniques. The analysis matches the typical behavior of hydrotalcites of the same composition reported in the literature. From the analysis of the modified hydrotalcites, it can be concluded that carbonate was successfully intercalated into the structure of hydrotalcite using coprecipitation methods. It was also evidenced that there is an intermediate structure for carbonate intercalated hydrotalcite. X-ray diffractograms showed the formation of hydrotalcites in the synthesized materials, while EDAX analysis confirmed the expected elemental ratio present in the In the current study, we prepared hydrotalcite by the co-precipitation method and characterized using various analytical techniques. The analysis matches the typical behavior of hydrotalcites of the same composition reported in the literature. From the analysis of the modified hydrotalcites, it can be concluded that carbonate was successfully intercalated into the structure of hydrotalcite using co-precipitation methods. It was also evidenced that there is an intermediate structure for carbonate intercalated hydrotalcite. X-ray diffractograms showed the formation of hydrotalcites in the synthesized materials, while EDAX analysis confirmed the expected elemental ratio present in

**% removal of thiophene**

**0% 10% 20% 30% 40% 50% 60% 70% 80%**

precipitation methods. It was also evidenced that there is an intermediate structure for carbonate

the modified hydrotalcite adsorbents. We further confirmed the presence of ZnO, NiO, and zirconia on the respective adsorbents. FTIR and Raman spectroscopy showed the presence of carbonate and water in the adsorbents. The removal efficiency of thiophene (100 ppm) from the model fuel containing n-pentane by the five adsorbents was in between 70%–85%. However, only NiO-ZrO<sup>2</sup> adsorbent could remove thiophene from 100 ppm to a 15-ppm level. The thiophene removing capacity of NiO-ZrO2, and NiO-ZnO adsorbents is superior to the other three absorbents (ZnO, NiO, and ZrO2). NiO-ZnO is a well-known adsorbent, however, the adsorption capacity of NiO-ZrO<sup>2</sup> presented in this work can an have potential advantages over other adsorbents including the former one. Additionally, XPS showed clear evidence of the retention of sulfur from dibenzothiophene when interacted with the NiO-ZnO, NiO, and ZnO containing modified hydrotalcite adsorbents. Our results also indicated that the deposition of nickel on zirconia changed the nature of the active metal sites into the structure that successfully increased the sulfur removal capacity of the adsorbents.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-7080/8/2/22/s1, Figure S1: Infrared peaks of (a) Ni adsorbent; (b) Zn adsorbent (c) Ni and Zn adsorbent, (d) Zr adsorbent, Figure S2: Raman peaks of calcined (a) Ni adsorbent; (b) Zn adsorbent (c) Ni and Zn adsorbent, (d) Zr adsorbent, Table S1: Peak center, FWHM and crystal size from X-ray diffractograms, Table S2: Composition of the five calcined adsorbents from EDAX.

**Author Contributions:** Conceptualization, M.M., R.D., M.I. and A.J.G.; methodology, M.M. and A.R.; software, R.D. and A.R.; validation, M.S.J.S., A.J.G. and A.R.; formal analysis, M.M., R.D. and M.I.; investigation, A.G.; resources, A.G.; data curation, M.M. and R.D.; writing—original draft preparation, M.M., R.D., M.S. and A.R.; writing—review and editing, M.M., R.D., M.S.J.S., T.S. and A.R.; visualization, M.M., R.D., M.S.J.S., T.S. and A.R.; supervision, A.J.G. and A.R.; project administration, A.J.G. All authors have read and agreed to the published version of the manuscript.

**Acknowledgments:** This work has been funded in part with funds from the State of Texas as part of the program of the Texas Hazardous Waste Research Center, Texas Air Research Center, Lamar University, and Prairie View A&M University. The contents do not necessarily reflect the views and policies of the sponsor nor does the mention of trade names or commercial products constitute endorsement or recommendation for use. The authors would like to extend special thanks to Hylton G. McWhinney, David Cocke, Paul Bernazzani, Hylton McWhinney, Tony Graddy and Tracy Benson for their analytical support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Process Development of CO2-Assisted Polymer Compression for High Productivity: Improving Equipment and the Challenge of Numbering-Up**

## **Takafumi Aizawa**

Research Institute for Chemical Process Technology, National Institute of Advanced Industrial Science and Technology, 4-2-1 Nigatake, Miyagino-ku, Sendai 983-8551, Japan; t.aizawa@aist.go.jp; Tel.: +81-22-237-5211

Received: 12 April 2019; Accepted: 7 May 2019; Published: 8 May 2019

**Abstract:** The CO2-assisted polymer compression method is used herein to prepare porous polymer materials by bonding laminated polymer fiber sheets using a piston in the presence of CO2. In this work, the CO<sup>2</sup> flow line connections were moved from the pressure vessel to the piston to increase productivity, which makes the pressure vessel free-moving and the processing time of sample introduction and removal seemingly zero. In addition, a numbering-up method suitable for CO2-assisted polymer compression is proposed and verified based on the variability of the products. The variability of the product was evaluated using porosity, which is one of the most important properties of a porous material. It is found that the CO<sup>2</sup> exhaust process, specific to this method, that uses high-pressure CO2, causes product variation, which can be successfully suppressed by optimizing the CO<sup>2</sup> exhaust process.

**Keywords:** CO2-assisted polymer compression; numbering-up; high productivity; CO2; polymer; porous material; process improvement

#### **1. Introduction**

Polymers are an essential material in everyday life due to their light and durable characteristics [1,2]. In particular, porous polymer materials are very lightweight and exhibit various functionalities [3] such as sound absorption, heat insulation, adsorption, filtering, moisture permeation, water absorption, drug loading, and sustained release [3]; their applications vary from mass production to advanced functional materials [4].

Considering the advantages of such materials, CO2-assisted polymer compression (CAPC), a method that uses CO<sup>2</sup> to adhere fibrous sheets to create a porous polymer material, has been developed [5]. Because the method is a room temperature process and no heater is required, it is a low energy process. As the CO<sup>2</sup> used is released into the atmosphere and does not remain in the polymer, impurities do not remain. In addition, CO<sup>2</sup> is a very low toxic substance; it is also used in the food industry, e.g., carbonated beverages. CO<sup>2</sup> is known as an environmentally friendly solvent; in particular its supercritical state is recognized as a green solvent [6], and its application in various processes has been researched [7,8]. CO<sup>2</sup> has also been practically used as an extraction solvent for biomass in the food industry [9]. The demonstrated CO<sup>2</sup> process is a suitable method for manufacturing parts used for the food industry, medical treatment, and living environments that demand high purity. Among gases, CO<sup>2</sup> is known to dissolve well in polymers and has been reported as a specific interaction between the CO<sup>2</sup> and polymer [10], estimation of the solubility of CO<sup>2</sup> in polymers [11], measurement of solubility of CO<sup>2</sup> in crystalline polymers [12], and measurement of solubility of CO<sup>2</sup> in amorphous polymers [13]. CO<sup>2</sup> is known to be capable of plasticizing polymers and decreasing their viscosity [14], melting point [15], and glass transition temperature [16]. Many CO2-based processes with polymers have been

proposed such as dyeing and polymer particle synthesis [17], foam molding [18], fabrication of carbon nanotube composite foam [19], foam molding of polymer blends [20], molecularly imprinted polymer development in supercritical CO<sup>2</sup> [21], and fabrication of two-dimensional porous polymers [22]. As the CAPC process is performed at room temperature with no temperature controller and the pressure of CO<sup>2</sup> is also introduced at the vapor pressure, which is the pressure inside the liquid CO<sup>2</sup> cylinder, the CO<sup>2</sup> is conveniently introduced in the gas phase by opening a valve. As no pump is required, the process is simplified. This technology reduces equipment cost and processing time. Previous work has also shown that material properties, such as porosity, were easily tuned by changing the process conditions, even when using the same raw material [23]. The adhesion strength of the CAPC sample was also evaluated [24]. Others have demonstrated that a drug can be carried in the porous material very easily by loading the drug in the state of the fiber sheet of the raw material, and the sustained release property of the drug when the porous material is immersed in water can be controlled [25]. There are many applications of drug-loaded porous materials in the advanced materials and medical fields, but it is challenging to insert the drug into the center of a thick porous material. The CAPC process facilitates drug carrying by supporting drug insertion at the raw fiber stage.

In my previous studies, only one sample was fabricated in one batch, which is a disadvantage remedied in this work. The CAPC method is used to prepare a porous polymer material by adhering laminated fibrous sheets through a simple apparatus with a short pressing time. Even if it is an overwhelmingly fast process at the laboratory level, when considering mass production, the improvement in the production per unit time is indispensable. Historically, scaling chemical processes up is often challenging. One of the main research goals in chemical engineering is developing ways to smoothly scale up flask-scale reactions to be successful in large reactors.

For processes that cannot be scaled up, numbering-up was proposed. The concept of numbering-up has been structured with a chemical process using a microreactor as shown by numbering-up for multi-phase (gas–liquid) flow microreactor [26], numbering-up for liquid–liquid two-phase capillary microstructured reactor [27], and external numbering-up (parallelized microdevices) [28]. Using a small reactor (microreactor) and a small mixer (micromixer), the microscale chemical process maximizes the heat exchange rate and mixing speed and realizes a high-speed reaction that suppresses by-products. Because the reactors and the mixers must be small, scaling up cannot be performed. Numbering-up that increases the number of reaction lines is important. In manufacturing, it has been common practice to approximate numbering-up conventionally. Adding production lines to increase production is one possible numbering-up strategy. However, a new concept that has emerged in microchemistry is the introduction of the number of reaction lines in an apparatus, wherein numbering-up is considered inside the apparatus.

The micro chemical process is a flow process, but the concept of numbering-up is also adaptable to a batch-type production process in that there is a limit to the amount that can be produced in one reactor. In a batch production process, if the size of the batch is optimized, it is impossible to scale up, and the idea of numbering-up is necessary. When considering production by the CAPC method, scaling up can be applied when large porous polymer materials are fabricated, and the final products are cut into pieces. When it is desired to mold the polymer in the form of a final product, for example, when preparing a drug-loaded tablet, it is necessary to consider numbering-up to increase the production per unit time. The chemical reaction in the flow process and the compression of the polymer fiber sheets are completely unrelated phenomena. Therefore, only the concept of numbering-up for simultaneous production in the same equipment was considered herein.

This study is the first trial to enhance the productivity of the CAPC process. Herein, increasing the productivity of the CAPC process is investigated through two strategies. The flow line of CO<sup>2</sup> is improved to hide both sample preparation time and sample removal time. The application of numbering-up such that the cost of equipment is suppressed is also investigated. Product uniformity is evaluated using the porosity of the samples when numbering-up is incorporated. Porosity is an important property of porous materials, and previous research has shown that the sustained release

rate of the drug trapped in the porous material is greatly influenced by the fabricated material's porosity [25]. The cause of sample variation is identified and a solution to suppress sample variation is presented.

#### **2. Materials and Methods**

#### *2.1. Materials*

A nonwoven fabric spun with a fiber diameter of 8 µm and a basis weight of approximately 30 g/m<sup>2</sup> using a melt-blown method (Nippon Nozzle Co., Ltd., Kobe, Japan) was used as a raw material; a polyethylene terephthalate (PET) pellet (Model number: TK3) of Bell Polyester Products Inc. (Yamaguchi, Japan) was used. The procedure of the melt-blown process is very simple. The PET pellets were melted, extruded from small nozzles, and stretched using high-speed blowing air to create a nonwoven fabric. Circular cut samples (18 mm diameter) of the nonwoven fabric were selected such that the same weight (0.257 g) was obtained for each set of 32 sheets. μ

#### *2.2. Equipment*

The experimental equipment used in this study is presented in Figure 1. The major difference from previously reported research [5,23–25] is found in the connections of the CO<sup>2</sup> introduction exhaustion tubes. The position of the tubes was moved from the high-pressure vessel (B1) to the piston (P). The piston area includes a CO<sup>2</sup> introduction port and an exhaust port at the top, a flow path for carrying CO<sup>2</sup> at the center, and diagonally positioned CO<sup>2</sup> introduction and exhaustion ports into the gap between the high-pressure vessel and piston at the bottom. Owing to the difficulty in creating such a shape with one piece, the system comprised mainly three pieces as shown in Figure 2. The inner diameter of component (P2) through which CO<sup>2</sup> flows is 4 mm, and the outer diameter of component (P3) is 19.5 mm. Black circles in Figure 2 indicate the O-rings. Component (P1) and component (P2) are sealed with an O-ring. An O-ring between components (P2) and (P3) seals the pressure vessel and component (P2). There is no gas seal between component (P2) and component (P3). The piston was attached to the press machine (Model number: JP-1504), manufactured by Janome Sewing Machine Co., Ltd (Hachioji, Japan). The inner diameter of the pressure vessel is 20.0 mm. Both the piston and the pressure vessel are made of stainless steel (SUS316).

**Figure 1.** Schematic illustration of the cross-section of the high-pressure vessel used for CO<sup>2</sup> -assisted polymer compression. B1: Body of the high-pressure vessel, B2: Base of the high-pressure vessel, C: CO<sup>2</sup> cylinder, P: Piston, PC: Laptop computer, S1: Sample, S2: Separator, V1: Intake valve, V2: Exhaust valve, V3: Exhaust valve (optional), and V4: Metering valve (optional).

**Figure 2.** Cross-sectional diagram of the piston components. The piston comprises three main components, P1–P3.

The liquid CO<sup>2</sup> cylinder was connected to the CO<sup>2</sup> inlet of the piston through the introduction valve by a stainless-steel (SUS316) tube. The tube was slacked between the introduction valve and the piston so that the upward and downward movement of the piston was not affected. The stainless-steel tube connected to the outlet of the piston was also connected to the exhaust valve with slack. After the exhaust valve, a T-type branch was provided; one was connected to an additional exhaust valve (V3), and the other was connected to a metering valve (V4). The flow line after the exhaust valve (V2) (dotted red frame in Figure 1) was not provided in the first experiment, and the exhaust was directly released to the atmosphere. Although the introduction valve (V1) was closed during operation, it was temporarily opened when CO<sup>2</sup> was introduced. V2 was normally open, but it was closed when CO<sup>2</sup> was to be retained in the pressure vessel. In operating V4 and V3, when V2 was open and V3 was closed, CO<sup>2</sup> was slowly released through the metering valve. A metering valve (Model number: SS-SS1, Swagelok Inc., Solon, Ohio, United States) was set to the minimum flow rate and used as V4. When V3 was open, CO<sup>2</sup> was instantaneously exhausted. When the experiment is performed without the process outlined in the dotted red frame (see Figure 1), it appears that the same result will be achieved as that when V3 is open. After performing the experiment without the dotted red frame line, the flow line in the red frame was added. This has been depicted as the experiment without a dotted red frame line such that the performed experiment is represented accurately.

The high-pressure container is a container with a hole at the top. The structure was made so that the body part (B1) and the base part (B2) could be separated if the system became clogged. By fixing the container filled with the sample just beneath the piston, the setup of the device is complete.

The laptop computer (PC) controlled each valve and piston such that the timing of the valves opening and closing was consistent for all experiments.

#### *2.3. Experimental Procedure*

Five sets of 32 laminated samples (S1) were placed in the high-pressure vessel. The laminated samples were separated by stainless-steel spacers (S2), each having a diameter of 19.5 mm and a thickness of 1 mm. As the inner diameter of the pressure vessel is 20.0 mm, there is enough clearance for CO<sup>2</sup> to flow between the pressure vessel and the spacers. The inner diameter of the stainless tube from the CO<sup>2</sup> cylinder to the piston is 0.8 mm, and the clearance between the pressure vessel and the spacers is much larger than that. The cross-sectional view of the samples when set in the high-pressure vessel is shown in Figure 1. The measured thickness of the center portion of each spacer was, starting from the top spacer, 1.024 mm, 1.009 mm, 1.010 mm, and 1.035 mm. Operation was started when V1 closed and V2 opened. Then, the piston was lowered to the CO<sup>2</sup> introduction position, and V2 was closed. The introduction valve and the exhaust valve were alternately opened three times to replace the air in the pressure vessel with CO2. As the internal pressure of the liquid CO<sup>2</sup> cylinder is about 6 MPa, once CO<sup>2</sup> is introduced at this pressure and released into the atmosphere (0.1 MPa), the amount of existing air is reduced to 1/60, and the remaining air is replaced with CO2. The volume of air is estimated to become (1/60)<sup>3</sup> following three trials of the procedure. The CO<sup>2</sup> introduction valve was opened again to introduce CO2, and the introduction valve was closed. Then, the piston was lowered to the press position and set to press for 10 s, after which an exhaustion procedure was performed. The exhaustion operation was conducted in two ways depending on the experiment. One exhaust procedure was an experiment without the flow line enclosed by the dotted red line indicated in Figure 1. In this case, CO<sup>2</sup> was exhausted instantaneously (<1 s). In the other procedure, the flow line enclosed by the dotted red line was attached; V2 was opened when V3 was closed, and CO<sup>2</sup> was slowly exhausted for 30 s by V4. Then, V3 was opened to exhaust CO<sup>2</sup> completely. The distances (when the thickness of the spacer is excluded) between the CO<sup>2</sup> introduction position and the bottom of the pressure vessel and that between the press position and the bottom of the pressure vessel were 9.75 mm and 6.50 mm, 9.00 mm and 6.00 mm, 8.25 mm and 5.50 mm, 7.50 mm and 5.00 mm, and 6.75 mm and 4.50 mm. The experimental results correspond to press positions from 4.50 mm to 6.50 mm. For evaluating the sample, the thickness of the center was measured with a micrometer and the variation due to the sample position was assessed.

Porosity was evaluated from the polymer density, weight of the fabric sheet, and sample thickness. The datasheet for the PET pellet indicated that the density of the used polymer is 1.34 g/mL. For a solid without any voids, the thickness can be easily calculated using the weight, density, and diameter of the fabric sheet. In this experiment, laminated 0.257 g of circular sheets with an 18.0-mm diameter were used, such that the thickness of the solid (*L*solid) was calculated to be 0.754 mm. The difference between *L*solid and the actual thickness (*L*sample) was considered as the total area of the pore, and the porosity α*,* was calculated as α = (*L*sample – *L*solid)/*L*sample.

The experiment was performed seven times under each condition and the average value was used for analysis. Error bars are not shown in the plot of the evaluated porosity because the standard deviation of the data was small (maximum 4.9 %); the standard deviation is provided in Tables 1 and 2.

The fiber surface of the central layer of the compressed sample was observed using a scanning electron microscope (SEM; Model number: TM-1000, Hitachi High-Technologies Co., Minato-ku, Japan). The fibrous sheet before compression was also measured using an SEM for comparison.

#### **3. Results and Discussion**

First, the design of the system, specifically the effect of connecting the CO<sup>2</sup> lines to the piston, is discussed. When they were connected to the high-pressure vessel, the vessel could only be moved within the range that the tubes could extend. To pull the pressure vessel out completely, it was necessary to disconnect the tubes from the pressure vessel, resulting in reduced usability. By connecting the tubes to the piston, the high-pressure vessel became free of this limitation. Moreover, the initial setting (corresponding to the beginning of the operation) was complete by locating the high-pressure vessel including the samples under the piston. If a plurality of high-pressure vessels were set up, it would be possible to prepare and remove samples from the system during press operation, which would greatly improve productivity.

Despite reducing the time for setting and removing the samples, the time required to move the piston, replace the air with CO2, and press remains absolutely necessary time, which is nearly 30 s. If the processing time per sample should be reduced further, multiple samples must be fabricated simultaneously; at this point, the concept of numbering-up is attractive. Figure 3 schematically presents

the preparation of one sample and two kinds of numbering-up. Figure 3b shows the concept of preparing a plurality of pressure vessels and pressing a plurality of samples with one press machine. In this case, numbering-up can be realized under the same conditions as the preparation of one sample and it is possible to increase the production efficiency while suppressing variations in the sample. On the other hand, the maximum load determines the performance of the press machine, and when the press area is increased, the required maximum load for the press machine increases, which increases process cost dramatically and should thus be avoided in the practical process. In Figure 3c, numbering-up by stacking samples without changing the press area is shown. Although the stroke of the piston will be long, if the maximum load does not change, the increased cost for the press machine-related owing to the longer stroke will be small; this will be very advantageous. Porosity is an important characteristic of porous materials. When a sample is compressed to fabricate a porous material, the thickness after compression determines the porosity. Thus, to produce samples of the same porosity, the compression rates must be uniform. In the method shown in Figure 3c, the separators are sandwiched between the samples, but the position of the separator is not controlled. Thus, the compression rates of the samples separated by the separators could differ. Owing to low or no uniformity of the raw material sample, the flow of CO<sup>2</sup> inside the pressure vessel, and the difference of the distance between the piston (moving device) and the samples, the products may vary. Therefore, product variation was evaluated using the porosity of the product as an index. As the objective here is to increase productivity, first, the experiment wherein CO<sup>2</sup> is instantly discarded by opening V2, i.e., without flow line shown within the dotted red line in Figure 1, was conducted.

**Figure 3.** Schematic drawings showing the (**a**) basic process, (**b**) parallel numbering-up, and (**c**) serial numbering-up.

Figure 4 shows the SEM images of the raw material sheets and CAPC samples. In the CAPC samples, dents, which are adhesion traces are observed on the fibers. The high compression sample (Figure 4b) has more dents than the low compression sample (Figure 4c), which is caused by strong compression. The origin of the pores is the void between the fibers, implying that the porosity depends on the strength of the compression. As the mechanism of the CAPC process has been previously described [5,23], only a brief explanation is provided herein. The shape of a thermoplastic polymer is maintained by the friction between the polymer chains. The fluidity of the polymer increases when the polymer is impregnated with CO<sup>2</sup> as the CO<sup>2</sup> molecules reduce the friction between the polymer chains. In the case of PET, CO<sup>2</sup> is considered to impregnate the amorphous part of the polymer. By compressing the polymer fibers in a plasticized state, bonding points are formed at the overlapping

portions of the fibers. The adhesion of the sample is completed by removing the CO<sup>2</sup> and the friction of the polymer chains is recovered.

**Figure 4.** Scanning electron microscope images of (**a**) raw material, (**b**) sample at 4.5 mm press position, and (**c**) sample at 6.5 mm press position.

Figure 5 shows the product porosity for different press positions. When the press position was small (high compression), the variation in the sample porosity was small. For larger press positions, larger variation was observed. The top and second samples have higher porosity than the lower samples. The experiments were performed seven times under each condition, and the average thicknesses and porosity are provided along with the standard deviation in Table 1. As seen in Figure 4 and Table 1, the differences in the sample porosity were significant.

**Figure 5.** Product porosity at each position for different press positions. Each value is the average thickness based on seven experiments.

I considered the time when variation occurs: When CO<sup>2</sup> is introduced, when it is pressed, or when the CO<sup>2</sup> exhaust is vented from the system. Introduced CO<sup>2</sup> may apply pressure throughout the samples from above by transmitting through the gap between the piston and the high-pressure vessel. However, there is a gap between the spacer and the high-pressure vessel, and CO<sup>2</sup> should flow downward instantaneously through this gap. In addition, as the polymer is plasticized only when CO<sup>2</sup> is dissolved in it, the effect on the polymer before CO<sup>2</sup> is dissolved is expected to be extremely low. When the piston is lowered after introducing the CO2, compression occurs. Considering that the piston moves suddenly when the sample and the spacer are stationary, the upper sample must

be compressed more than the lower sample. Another explanation is that the lower sample may be more pressurized by the weight of upper samples and spacers. However, as the pressure of the piston used to compress the sample is about 1 MPa (1 MPa at 4.5 mm press position, 0.9 MPa at 5.0 mm press position, and 0.8 MPa at 5.5 to 6.5 mm press position), the weight of samples and spacers should be negligibly small. The last possibility is to consider how CO<sup>2</sup> venting affects the system. When CO<sup>2</sup> exhaust is vented by releasing it to the atmosphere, the CO<sup>2</sup> trapped in the pores of the porous material and the CO<sup>2</sup> dissolved in the polymer are vigorously removed. If this effect differs for each sample layer, the samples may vary. Specifically, CO<sup>2</sup> explosively blows out from the upper sample near the outlet, which compresses the lower sample. Because the lower sample, in which CO<sup>2</sup> is sufficiently dissolved in the polymer, is still plasticized, the sample is likely to compress. When CO<sup>2</sup> is removed from the lower sample, which might compress the upper sample, it seems that the CO<sup>2</sup> concentration in the upper sample had already decreased; as decreased CO<sup>2</sup> would prohibit plasticization, the upper sample would be difficult to compress, thus leading to it being thicker.


**Table 1.** Product thickness and porosity obtained under rapid CO<sup>2</sup> exhaustion.

Assuming that the sample variation was most likely due to how the CO<sup>2</sup> exhaust was vented from the system, the experiment was conducted with an added metering valve; the metering value should enable even CO<sup>2</sup> discharge from each layer by limiting the exhaustion speed. The experiment was performed under the condition of the press position 6.5 mm, for which the dispersion had been largest. As described in the experimental section (Section 2.3), the evacuation operation was performed in a sequence of exhausting the residual pressure instantaneously after exhausting slowly for 30 s. According to the change in the load applied to the piston, the drop in pressure due to CO<sup>2</sup> exhaustion for 30 s was about 4 MPa. The measured sample thickness and evaluated porosity versus sample position are listed in Table 2. Although the porosity of the top sample was slightly larger, the variation among the samples was greatly reduced. Therefore, these results indicate that the sample variation was

caused by the CO<sup>2</sup> venting process. The time required for slow venting doubled the processing time (from 30 s to 60 s), which halves the effect of the 5-fold improvement in the productivity achieved via numbering-up. As the slow venting procedure performed to decrease the sample variation reduces the productivity improvement due to numbering-up, allowable product variation should be considered carefully before adapting slow venting.


**Table 2.** Product thickness and porosity obtained under slow CO<sup>2</sup> venting conditions.

#### **4. Conclusions**

Two strategies were conducted to improve the productivity of the CAPC method. The CO<sup>2</sup> introduction and exhaustion tubes were connected to the piston instead of the pressure vessel, enabling the free movement of the pressure vessel, which removed the processing time otherwise required for sample introduction and removal. The production per unit time was increased by numbering-up vertically and setting multiple samples without changing the load of the piston. This method is very effective because productivity can be improved without a significant increase in equipment cost. However, sample variation in porosity was detected in samples prepared in the numbering-up scheme, which was found to be caused by the CO<sup>2</sup> venting process. Although controlling the venting rate was one of the solutions, it was accompanied by a decrease in productivity. Therefore, carefully judging the allowable product variation together with the benefits of numbering-up is necessary. Although not discussed in this article, in the future, the flow of CO<sup>2</sup> at the time of exhaustion should be simulated to optimize the high-pressure vessel, spacer, and CO<sup>2</sup> flow path, such that CO<sup>2</sup> is removed at equivalent levels from all samples. This will potentially allow the sample variations to be lowered without reducing the CO<sup>2</sup> venting speed.

**Author Contributions:** T.A. conceived and designed the experiments, performed the experiments, analyzed the data, and wrote the paper.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **E**ff**ects of the Infill Density on the Mechanical Properties of Nylon Specimens Made by Filament Fused Fabrication**

#### **Svetlana Terekhina 1 , Innokentiy Skornyakov 2 , Tatiana Tarasova <sup>2</sup> and Sergei Egorov 2, \***


Received: 30 June 2019; Accepted: 13 August 2019; Published: 16 August 2019

**Abstract:** Additive manufacturing of polymer products over the past decade has become widespread in various areas of industry. Using the fused filament fabrication (FFF) method, one of the most technologically simple methods of additive manufacturing, it is possible to produce parts from a large number of different materials, including wear-resistant nylon. The novelty of the work is properties investigation of ±45 ◦ filling configuration with different filling degree for nylon, as well as calculating the effect of infill on the strength characteristics, excluding the shell. This article reflects the process of manufacturing samples from nylon using FFF technology with various internal topologies, as well as tensile tests. The analysis of the obtained results is performed and the relationship between the structure of the sample and the limit of its strength is established. To calculate real filling degree and the effect of internal filling on the strength characteristics of the specimen, it is proposed to use a method based on the geometric and mass parameters. The FFF method is promising for developing methods for producing a composite material. The results of this article can be useful in choosing the necessary manufacturing parameters.

**Keywords:** additive technologies; additive manufacturing; FFF; 3D printing; nylon

#### **1. Introduction**

Rapidly developing for several decades, additive technologies are gradually replacing the classic ways of making products in many industries. This fact is due to the main principle of new methods, the layer-by-layer creation of objects based on digital three-dimensional model controlled by a computer. This approach minimizes the number of equipment and the number of technological operations.

A wide range of materials available for additive production (from metallic powders to polymer filaments) makes it possible to produce analogs of products obtained by classical methods, often not inferior to them by mechanical characteristics. A bright representative of polymer materials used in this kind of production is nylon. Nylon is a thermoplastic polymer belongs to the group of polyamides. These plastics are characterized by high wear resistance but increased hygroscopicity (the ability to absorb moisture). Nylon is used in the production of elements of friction pairs and in the medical industry for the manufacturing of prostheses. Nylon is used to make some parts of machines because it is inexpensive and durable. It is often used in the electronics industry as a nonconductive and heat-resistant material. This polymer is one of the main materials for additive plants using fused filament fabrication (FFF). FFF technology refers to the simplest methods of 3D printing. However, this technology allows to produce a wide variety of products, from mock-ups and prototypes to robust functional elements.

One of the ways to increase the competitiveness of the product is to reduce its weight and the cost of its production. In relation to additive production, both objectives can be achieved by replacing the solid product with a shell with the same geometry and dimensions. To ensure its strength, a structure of the same material is formed inside the shell. The volume fraction of this structure inside the shell can vary from 0% to 100%. Obviously, with increasing this value, the strength of the structure will increase. The purpose of this study is to identify the relationship between the volume fraction of this structure inside the shell (20–100%) and the strength of the sample. It should be noted that in this article the filling structure is considered as a composite material, one of the components of which is nylon and the other is air.

At present, all technologies of additive production are fully reflected in the literature [1–3]. The main attention is paid to the manufacturing of products from metal powders [4,5] using selective laser melting technology in view of the increased interest in obtaining functional parts of various machines and mechanisms. However, polymer composite materials, filled with fibers of various types and composition [6,7], already constitute a serious competition to metals. In this connection, on the wave of general interest in additive technologies, attempts are made to directly produce composite materials using 3D printing [8].

Nowadays, experimental and theoretical studies of the effect of the filling structure of parts printed using additive technologies on their properties are being conducted worldwide [9]. The materials chosen are metal alloys (steel [10], aluminum [11], and copper alloys, etc.) as well as polymers (ABS (Acrylonitrile butadiene styrene) [12], PLA (polylactide) [13], polyamides [14,15], etc.). In a number of works, the influence on the mechanical properties of various parameters of the printing process (speed, temperature, layer thickness [16–18], etc.) and the internal architecture (filling geometry [19], directionality [20], and density [21]) is investigated. As the mechanical properties under investigation, tensile [22] and flexural [23] ultimate strength and the modulus of elasticity are most often chosen. Often, a comparative analysis of various materials [24,25], technologies [26], structures [27,28], and filling directions [29] is carried out. In some works, the authors conducted a comparative analysis of experimental and theoretical data [30–32].

To understand the rationality of the transition from traditional manufacturing methods to new ones, the accumulation of experimental data is required. The world scientific community studies the properties of samples from unfilled polymers made by additive technology methods (particularly by the FFF method [33]), and mathematical models of their formation are created [34]. Much attention is paid to the most widely used plastics in FFF installations, ABS and PLA [35,36]. Less common nylon is used in the manufacturing of gears and friction pairs, where it experiences mainly cyclic loads. In this regard, it is more often subjected to fatigue testing [37]. However, the prospects for using nylon as a matrix of composite material are growing because of its high compatibility with biodegradable natural fibers. This can significantly expand the scope of the material. Tests of tensile strength of nylon, useful in further studies, are described in this article.

#### **2. Materials and Methods**

Nylon filament (NYLON, manufacturer-Print Product, Saint Petersburg, Russia) was used as a material for the production of samples (Table 1). Due to the high hygroscopicity of the material, the drying of the filament was carried out immediately prior to manufacturing. The drying of the filaments before printing was carried out at 60 ◦C in a vacuum oven for 6 h. All the samples were then stored in the dry atmosphere of a desiccator prior to testing.



To determine the required dependence, strength tests of a number of samples having different volume fraction of the filling structure were carried out.

Samples were manufactured on equipment developed at MSTU "STANKIN" using elements of the Prusa Mendel project. Manufacturing of products on this equipment is carried out using FFF technology (fused filament fabrication). FFF is one of the methods of additive manufacturing, which consists of feeding a threadlike material into the heated chamber, melting it, squeezing it out through the nozzle, and depositing it onto the working surface. Simultaneous processes of extruding the filament, moving the extruder, and/or the working surface provide a layer-by-layer forming of the product. On the equipment of MSTU "STANKIN", these processes are controlled by computer using the software of Repetier Host. Moving of mobile elements is carried out according to the algorithm, previously formed by the program based on the three-dimensional model of the product and the characteristics set by the operator (Table 2).


**Table 2.** Characteristics of the fused filament fabrication (FFF) process.

The shape and dimensions of the samples were selected in accordance with GOST 11262-80 "Plastics. Method of tensile testing" (is an analog of ISO 527-2:2012) (Figure 1a). As a filling scheme, a perpendicular scheme was set. Each layer is filled with parallel tracks (strips of cooled polymer) positioned at an angle of 45◦ to the axis of the sample and perpendicular to the direction of the tracks of the previous layer (Figure 1b). Varying the distance between adjacent tracks results in a change in the volume fraction of the fill pattern. Several groups of samples were made, for each of which a theoretical volume fraction of the filling structure was set: 20%, 40%, 60%, 80%, and 100%.

(**a**)

**Figure 1.** *Cont.*

**Figure 1.** Sample in accordance with GOST (**a**), the layer infill scheme (**b**), real specimen (**c**).

The actual geometric and mass characteristics of the samples were previously determined by caliper measurement (accurate to 0.1 mm) and weighing on the Mettler Toledo XPE analytical scales (accurate to 0.001 g).

To determine the tensile strength of the manufactured samples, the electrodynamic testing system ElectroPuls Model E10000 for axial loading with twisting was used.

The tests of each sample were carried out in accordance with GOST 11262-80 "Plastics. The tensile test method" at a speed of 25 mm/min. Specimen after test is shown in Figure 2.

**Figure 2.** Specimen after test.

#### **3. Results**

After the tests, stretch diagrams of the samples were obtained (some of them are shown in Figure 3) and the strength limits were calculated.

**Figure 3.** Tensile diagrams of some samples.

As the degree of filling increases, the ability of the sample to plastic deformation increases.

The obtained data cannot fully characterize the relationship between strength and volume fraction of the sample filling pattern. To identify this relationship, the actual fraction of filling of each sample and the strength of the filling (without taking into account the strength of the shell) were calculated. Calculation of the actual fraction of sample filling was made on the basis of geometric and mass parameters.

∑ The fabricated sample consists of *N*<sup>P</sup> layers of thickness *t*, of which *N<sup>L</sup>* lower ones and *N<sup>U</sup>* upper ones have a continuous filling (Figures 4 and 5a). The number of lower layers varies from sample to sample due to different conditions for separating samples from the bed surface after manufacturing. Each of the *NINT* intermediate layers consists of a solid contour and filling structure with a theoretical volume fraction *AT%* (Figure 5b). Thus, the total number of layers: ∑

<sup>∑</sup> = + + ூே்

$$N\_{\Sigma} = N\_L + N\_{U} + N\_{\text{INT}} \tag{1}$$

**Figure 4.** Working cross-section of the sample.

**Figure 5.** Infill scheme of upper and lower layers (**a**); infill scheme of intermediate layers (**b**). ∑

Taking into account the actual dimensions of each manufactured sample using the AutoCAD 2016 software, the contour was drawn and the values of the areas of the different zones of the layer were determined (Figure 6):

**Figure 6.** Layer zones.

*SC*—solid contour area, *SF*—filling area,

$$\mathcal{S}\sum = \mathcal{S}\_{\mathbb{C}} + \mathcal{S}\_{F} \tag{2}$$

ி = ூே் ∙∙ி —total area of the continuous layer.

ௌெ = ( + )∙∙<sup>∑</sup> + ூே் ∙∙

The difference in these values in the samples is due to changes in the temperature state of the working zone during manufacturing.

=ௌெ ∙ ே + ி ∙ *ρ ρ* According to obtained values, the volumes of solid material and filling were calculated: *VSM* = (*N<sup>L</sup>* + *NU*)·*t*·*S* <sup>P</sup> + *NINT*·*t*·*SC*—volume of solid material, *V<sup>F</sup>* = *NINT*·*t*·*SF*—volume of filling.

ௌெ = ( + )∙∙<sup>∑</sup> + ூே் ∙∙ The measured mass *m* of the sample:

−ௌெ ∙ ே ி

=

$$
\mathcal{m} = V\_{\text{SM}} \cdot \rho\_N + V\_{\text{F}} \cdot \rho\_A
$$

=ௌெ ∙ ே + ி ∙

ி = ே + ூோ ∙ ி = ே ∙ ே + ூோ ∙ ூோ

ρ*N*—density of solid material,

=

ρ*A*—actual filling density.

ி = ூே் ∙∙ி

*ρ ρ* Thus, the actual filling density:

$$\rho\_A = \frac{m - V\_{SM} \cdot \rho\_N}{V\_F} = \frac{m - \left[ (N\_L + N\_U) \cdot t \cdot \mathbf{S\_\Sigma} + N\_{\rm INT} \cdot t \cdot \mathbf{S\_C} \right] \cdot \rho\_N}{N\_{\rm INT} \cdot t \cdot \mathbf{S\_F}} \tag{3}$$

 = ி = ூே் ∙∙ி The mass of the sample is composed of the weight of solid material and the weight of air:

$$
\boldsymbol{m}\_{\boldsymbol{F}} = \boldsymbol{m}\_{\boldsymbol{N}} + \boldsymbol{m}\_{\boldsymbol{A}\boldsymbol{I}\boldsymbol{R}}
$$

$$
\rho\_A \cdot V\_F = \rho\_N \cdot V\_N + \rho\_{AIR} \cdot V\_{AIR}
$$

ூோ ρ*AIR*—air density,

ே

ூோ ே

*VN*—volume of solid material in the filling,

*VAIR*—volume of air in the filling.

Taking into account *V<sup>F</sup>* = *V<sup>N</sup>* + *VAIR* and *V<sup>N</sup>* = *VF*· *AACT*% 100% , the formula was obtained:

$$
\rho\_A \cdot V\_F = \rho\_N \cdot V\_F \cdot \frac{A\_{ACT\%}}{100\%} + \rho\_{AIR} \cdot V\_F \cdot \left(1 - \frac{A\_{ACT\%}}{100\%}\right)
$$

*AACT*% = ρ*A*−ρ*AIR* ρ*N*−ρ*AIR* ·100%—actual volume fraction of the filling structure. Neglecting the density of air, the following dependence was obtained:

$$A\_{\rm ACT\%} = \frac{m - \left[ (N\_L + N\_{\rm II}) \cdot t \cdot \mathcal{S}\_{\Sigma} + N\_{\rm INT} \cdot t \cdot \mathcal{S}\_{\mathbb{C}} \right] \cdot \rho\_N}{N\_{\rm INT} \cdot t \cdot \mathcal{S}\_{\Gamma} \cdot \rho\_N} \cdot 100\% \tag{4}$$

Calculation of the tensile strength of the filling structure (without taking into account the strength of the shell) was done on the basis of the following formula:

$$
\sigma\_B = \sigma\_{\rm BN} \cdot A\_N + \sigma\_{\rm BF} \cdot A\_F \tag{5}
$$

σ*B*—tensile strength of the entire sample,

σ*BN*—tensile strength of solid nylon (due to the change in material properties because of thermal stresses, the maximum tensile strength of 100%-filled sample was used),

σ*BF*—tensile strength of the filling structure,

*AN*—percentage of solid nylon contour in cross-section,

*AF*—percentage of filling structure in cross-section.

Taking into account the geometrical parameters of each sample (number of upper, intermediate, and lower layers, working cross-section width, track width (Figure 3)), the values of *A<sup>N</sup>* and *A<sup>F</sup>* were calculated:

$$A\_N = \frac{(N\_{II} + N\_L) \cdot b + \mathbf{2}b\_{TR} \cdot \mathbf{N\_{INT}}}{(N\_{II} + N\_L + N\_{INT}) \cdot b} \tag{6}$$

$$A\_F = \frac{\text{N}\_{\text{INT}} \cdot (b - 2b\_{\text{TR}})}{(\text{N}\_{\text{U}} + \text{N}\_{\text{L}} + \text{N}\_{\text{INT}}) \cdot b} \tag{7}$$

The obtained dependence for the tensile strength of the filling structure:

$$
\sigma\_{\rm BF} = (\sigma\_{\rm B} - \sigma\_{\rm BN} \cdot \frac{(N\_{\rm II} + N\_{\rm L}) \cdot b + 2b\_{\rm TR} \cdot N\_{\rm INT}}{(N\_{\rm II} + N\_{\rm L} + N\_{\rm INT}) \cdot b} \cdot \frac{(N\_{\rm II} + N\_{\rm L} + N\_{\rm INT}) \cdot b}{N\_{\rm INT} \cdot (b - 2b\_{\rm TR})} \tag{8}
$$

Based on the obtained values of the volume fractions (theoretical and actual) and the strength of the filling structure of each sample, the arithmetic mean values were determined (Table 3) and the dependence was obtained (Figure 7).

**Table 3.** Average arithmetic values of the volume fraction and the tensile strength of the infill structure.


**Figure 7.** Dependence of the tensile strength on the actual volume fraction of the infill structure.

#### **4. Discussion**

The obtained results show that with the increase of the volume fraction of the infill structure, the discrepancy between the preproduction (theoretical) and actual values is growing. This can be explained by the fact that the material has time to partially cool before reaching the bed, and is stacked as an oval-shaped track (Figure 8). As a result, voids are formed, leading to a decrease in the volume fraction of plastic. The increase of the volume fraction and, as a consequence, the supply of material, leads not so much to a decrease of these voids as to an increase of the cross-section of the sample. Thus, the conducted experiment did not allow to establish a dependence in the area from 72% to 100% of the infill, which does not make it possible to compare samples fabricated by the FFF technology and cast samples.

**Figure 8.** Tracks interaction of different theoretical volume fractions of the infill structure: 40% (**top**), 80% (**left**), 100% (**right**).

Nevertheless, this dependence can be useful in choosing the volume fraction of the infill structure in the obtained range of values depending on the required strength characteristics.

Carrying out a qualitative assessment of the results, it can be noted that after increasing the volume fraction of the infill structure above 60%, a significant increase in strength occurs (Figure 7). The ultimate strength of the fabricated samples is determined not only by the amount of material but also by the contact between the parallel tracks. Observation of the manufacturing process showed that when setting a theoretical volume fraction of filling in the range of 20–40%, neighboring tracks of the same layer do not touch each other. When the parameter is increased to 60% (which corresponds to the actual value of 54%), the parallel tracks contact (Figure 8), which leads to the formation of a continuous layer and increases the strength of the entire sample.

Thus, the observation of the production of samples from nylon and the tensile tests made it possible to establish a relation between the internal structure and the strength properties of the samples, and to characterize the effects that occur during the printing process.

**Author Contributions:** Conceptualization, S.T.; data curation, S.T.; formal analysis, S.T.; funding acquisition, T.T.; investigation, I.S. and S.E.; methodology, S.T. and I.S.; project administration, T.T.; resources, I.S. and S.E.; software, I.S.; supervision, T.T.; validation, S.T. and T.T.; visualization, S.E.; writing—original draft, I.S.; writing—review and editing, T.T.

**Funding:** The work was carried out with the financial support of the Ministry of Education and Science of the Russian Federation in the framework of state task No. 11.1267.2017/4.6.

**Acknowledgments:** The studies were carried out using the equipment of the Center for Collective Use of MSTU "STANKIN". The authors express their gratitude to the anonymous reviewers, whose remarks made it possible to improve the article.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


= Mechanical Engineering: Traditions and Innovations. In Proceedings of the Scientific and Technical Conference, Moscow, Russia, 25–26 October 2016; pp. 86–88. (In Russian)


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Compensation for Geometrical Deviations in Additive Manufacturing**

**Christoph Hartmann 1, \* , Philipp Lechner 1 , Benjamin Himmel 1 , Yannick Krieger 2 , Tim C. Lueth <sup>2</sup> and Wolfram Volk 1**


Received: 28 October 2019; Accepted: 28 November 2019; Published: 2 December 2019

**Abstract:** The design of additive manufacturing processes, especially for batch production in industrial practice, is of high importance for the propagation of new additive manufacturing technology. Manual redesign procedures of the additive manufactured parts based on discrete measurement data or numerical meshes are error prone and hardly automatable. To achieve the required final accuracy of the parts, often, various iterations are necessary. To address these issues, a data-driven geometrical compensation approach is proposed that adapts concepts from forming technology. The measurement information of a first calibration cycle of manufactured parts is the basis of the approach. Through non-rigid transformations of the part geometry, a new shape for the subsequent additive manufacturing process was derived in a systematic way. Based on a purely geometrical approach, the systematic portion of part deviations can be compensated. The proposed concept is presented first and was applied to a sample fin-shaped part. The deviation data of three manufacturing cycles was utilised for validation and verification.

**Keywords:** additive manufacturing; 3D-printing; compensation; accuracy; precision

#### **1. Introduction**

A central, key technology in the future will be additive manufacturing (AM) [1]. Opportunities are arising in combination with ongoing digitisation and cross-linking between production and design in line with industry 4.0. For certain components this will allow one to avoid costly and time-consuming production processes, such as the manufacturing of moulds and tools. An increase in flexibility and the production of individualised products containing complex structures without additional effort go hand in hand. In this way, AM can make a significant contribution to mastering the complexity in the production process associated with the rapidly increasing number of variants [2].

However, along with all of the opportunities AM processes offer, a comprehensive transfer to mass production is nowhere in sight. Central hurdles for the application in serial production are above all, the cost-related general conditions (high material costs in combination with an insufficient ratio of plant costs to productivity) and the high manual effort in the process chain due to the lack of physical and digital line integration. This current lack of the integration of AM processes in conventional production environments is due to the specifics of AM processes, such as batch production, and the currently low degree of automation in the peripheral systems, such as component handling, quality measurement and transport chains. On the component level, one major critical point in mass production is accuracy and reproducibility. In AM processes, deviations between the nominal and the actual parts are primarily caused by the highly complex thermo-mechanical conditions leading to high gradients and high temperatures, for example. This significantly limits the use and the associated potential of AM processes for in series production [3].

#### **2. State of the Art**

When considering the dimensional quality of manufactured components, the systematic and stochastic portions of the deviations have to be carefully separated. Varying environmental conditions, material fluctuations, tribology, and changing thermo-mechanical boundary conditions are typical sources of stochastic errors. It is impossible to completely eliminate the stochastic portion in real applications and direct elimination is hardly possible. To cope with stochastic phenomena, process robustness, reasonable component design and the suitable definitions of tolerances are necessary. The systematic deviations can be related to identifiable causes and possess unique, unidirectional characteristics. For that reason, systematic deviations can be addressed by geometric compensation. For a separation of both deviation types, a sufficient number of measurements is necessary.

#### *2.1. Causes for Geometrical Deviations in Additive Manufacturing*

Due to the complexity of the process chain from the CAD model to the finished component, geometric deviations can occur at different points in the manufacturing process. These can be divided into the following error categories [4]:


In laser powder bed fusion (L-PBF), each layer is two-dimensional. The three-dimensional component is, thus, created step-by-step by joining individual layers of the same thickness along the Z coordinate. For this purpose, mathematical layer information is extracted from the digital 3D data set during pre-processing. This results in deviations during the creation of a discrete mesh (usually .STL). In addition, layer-by-layer construction can result in steps in a component, since the smallest possible resolution in the Z-direction is the layer thickness [4,5].

During processing, the powder is heated by a laser to a temperature above its melting temperature in order to bond the particles through phase transformation. There is a risk of thermal stresses due to an inhomogeneous temperature distribution in the build chamber. As a result, distortions occur in the component, as the material shrinks inhomogeneously during cooling, resulting in thermally induced stresses. Dimensional deviations resulting from material shrinkage are dependent on the geometry of the respective component [6,7].

Based on temperature measurements in an SLS machine, a significant temperature gradient in the x-y plane of a layer was determined. From this it can be concluded that the material shrinkage of a geometry additionally depends on its x and y-positions in the build chamber [6]. The dimensional accuracy in the x-y plane of a layer depends directly on the scanning system used by the AM system and its ability to move the laser beam along the desired path [8]. Such machine errors, on the other hand, are measurable and can be compensated to a certain extent in order to minimise their effect on the end product to an acceptable minimum. In addition to the factors already mentioned, despite an identical processing environment, quality fluctuations can occur between the end products, the exact cause of which cannot be defined, so they are, therefore, regarded as stochastic errors. Of all these factors, the temperature-induced shrinkage of the material is the main cause of inaccuracies in the end product [4,8].

#### *2.2. Classical Methods to Compensate for Deviations in Additive Manufacturing*

In the literature, different methods have been found to counter geometry deviations in the component by conventional means. However, these are usually limited to one specific error. In their work, Pham et al. analysed various factors that influence the accuracy of additive manufactured components. This began with the generation of the STL file. The chord height and the angle control have a big influence on the accuracy. Some of these error sources can also be avoided by using alternative data formats [4,8,9].

During further processing from 3D to layer information (slicing), discretisation errors occur, which are mainly due to the z-direction resolution, which is limited to one layer in thickness [10]. These effects can, however, be reduced by various algorithms; for example, by adjusting the layer thickness [11]. In addition, the orientation of the parts in the build chamber can be adjusted. Masood et al. describe in their research, the possibility of minimising the error caused by the step effect by adjusting the orientation in the build chamber. The best possible orientation can be calculated on the basis of the minimum volumetric error. This leads to an improvement of the accuracy and the surface quality of the component by minimising the effects of the stair tread effect [12].

In order to prevent geometric deviations due to residual stresses in the component, heat treatments are used to reduce mechanical stresses in the component. This prevents stresses from being released and dimensional deviations from occurring in one of the subsequent processing stages, such as separation from the construction platform. Reference [13] In addition, thermal distortion is influenced by the scanning strategy and scanning time and can be improved by adjusting them. References [14–16] Kamat and Pei analytically calculated deformations due to residual stresses in parts with internal channels, which can be used to compensate for deformations [17]. To compensate for shrinkage effects, different scaling factors are applied to the individual 2D layers. Furthermore, a compensation of the laser diameter can be achieved by an offset in the 2D planes [4,18,19].

#### *2.3. Geometrical Deviation Compensation in Additive Manufacturing*

Due to interacting factors, such as material shrinkage, geometric complexity, quality fluctuations and machine faults and inaccuracies, it is challenging to predict and control the extent of dimensional deviations. The straightforward way to compensate for deviations is currently to use a constant scaling factor across a part to calculate material shrinkage. Alternatively, different scaling factors can be applied to individual sections in the CAD model. It is assumed that the dimensional deviations occur evenly in the individual sections. However, this assumption is often not legitimate due to the complexity of some geometries, and this can lead to transfer effects, among other things [20].

Huang et al. deal with the compensation of dimensional deviations in the x-y plane and define various quality standards that minimise volume and surface deviations [20,21].

In addition, initial work has been done to use pre-deformation to compensate for dimensional deviations in selective laser beam melting. This involves simulating thermal distortion in the component using a finite element calculation, and then pre-deforming the CAD model against the expected deviations [22–24]. Zhu et al. proposed a statistical shape analysis, which is trained with data obtained by simulation and can predict geometrical deviations for similar parts [25]. Furthermore, they proposed a compensation method for the layers of the sliced part [26]. Schmutzler et al. compared a simulation-based pre-deformation to a deviation compensation in theory, which is based on previous manufactured iterations, for a cantilever geometry. They showed that the simulation-based compensation converged to the desired geometry within three iterations [27].

#### **3. Hypothesis**

Geometrical accuracy has always been one of the main objectives in component manufacturing. The established manufacturing processes, therefore, have a variety of approaches that address this problem. In forming technology, two systematic approaches to deviation compensation are the

displacement adjustment method and the stress-based adjustment method [28,29]. The stress-based method uses measured deviations to calculate a certain stress state on the target geometry. This stress state is reversed and is used as a boundary condition to deform the forming tool. The result of the deformation is the compensated forming tool. The displacement adjustment method works similarly. The measured deviations are directly reversed and used for the geometrical adaptation of the forming tool. Both methods have been used successfully in different applications. The stress-based adjustment methods seem to perform better for processes dominated by a global deformation pattern—sheet metal forming, for example, where the displacement adjustment method tends to be suitable for local compensation, as with bulk forming.

In AM, no tools are needed for manufacturing, but a manufacturing geometry has to be provided to the AM system. We proposed that this manufacturing geometry could be seen as a forming tool in a geometrical compensation sense. We further proposed that the displacement adjustment method, which directly incorporates the deviation data, is a suitable method for AM, because of its incremental characteristics.

In this paper, the objective was to prove those hypotheses. Therefore, a general geometrical compensation procedure based on the displacement adjustment method is suggested for AM to capture the problem of accuracy and reproducibility also with regard to an implementation in production lines. The compensation scheme is shown using a SLS process.

#### **4. Materials and Methods**

For the investigation in this paper, a SLS process is used to apply the conducted general geometrical compensation procedure. In the following subsections, the parts created, manufacturing steps and further processing are presented in detail.

#### *4.1. Specimen and Manufacturing*

A sample part was designed to show the effectiveness of the compensation algorithm. Figure 1 depicts the whole job. The fin-shaped geometry of the component was chosen specifically to boost thermal deviations. A total of 27 parts were positioned in the build chamber in three rows with nine parts each. We measured and compensated the blue coloured parts *ai*–*e<sup>i</sup>* . The light blue parts were measured but not compensated and served as a reference to their compensated direct neighbour to study the reproducibility of the process. For the digitisation of the components, the structured white light scanning device GOM ATOS II 400 (GOM, Braunschweig, Germany) was utilised. For each component, approximately 15 frames had to be acquired. With the setup we used, a resolution of 0.05 mm was achieved. The test specimens were fabricated using an EOS Formiga P100 (EOS, Krailing, Germany) in combination with polyamide 12 (EOS PA2200, Krailing, Germany). The powder we used had a 50/50 ratio of new and recycled powder. The chosen layer thickness was 0.1 mm. The job was positioned in the centre of the build platform with the standard offset of six millimetres in the z plane. The temperature of the process chamber was set to 172 degrees Celsius. The material-dependent scaling for the jobs was (X:3.086%; Y:3.019%; Z(0):2.2%; Z(300):1.6%). These machine parameters were kept constant for all jobs.

#### *4.2. Compensation Algorithm*

The compensation strategy was based on a framework proposed for bulk forming processes [30,31]. The compensation approach was fully driven by geometrical data, which, in principle, ensures general validity and scale invariance. In Figure 2 the flow chart for the compensation procedure is shown (forming in the upper half and the adaptation of AM in the lower half). Each process originates from the part design and the target geometry, which could be provided in any form; for example, as a CAD-design or reverse engineering measurement object. In case of forming processes a forming tool geometry is derived, which is close to the negative of the target geometry, but usually is divided into an upper and a lower tool part due to the uni-axial working direction of the presses. The active surfaces

of the forming tool in the closed state show a certain transformation of the target geometry's surface with the same topology. In the surface generation process (centre of Figure 2), this transformation, also called compensation, is deduced through the iterative tool redesign. Based on the target geometry for AM manufacturing, a discrete AM manufacturing geometry was determined based on the part design. This step corresponds to the generation of active tool surfaces in forming processes. Since AM is an incremental process, no physical tool is necessary, but an active closed surface is required, which is a building instruction for the AM system. It is transferred to the AM system for initial manufacturing. As an alternative, the framework is also capable of handling data in a numerical simulation framework. Whether part manufacturing, numerical analysis or both should used, depends on different process and part specific characteristics and attendant circumstances. Both from real manufacturing or numerical results, a final workpiece geometry was achieved and evaluated. The deviation data of the component compared to the target geometry served as input data for the systematic compensation. As a result, an updated AM manufacturing geometry is provided. This surface, in the sense of forming, "forms" the part, in the AM manufacturing environment. More abstract, in both manufacturing processes, surfaces have to be generated based on the target geometry, from which the actual manufacturing process result is controlled.

**Figure 1.** Stacking of the components in the build chamber. Blue components were measured and compensated; cyan components were measured for reference and not compensated.

The purely geometric approach assumes that all physical causes for deviations occurring during the manufacturing process are included in the systematic portion of the component shape. Hence, different physical boundary conditions have no effects on the general procedure. For example, different materials or processing temperatures can be treated in the same way. The approach also allows one, for example, to incorporate additional processing steps, such as heat treatments or machining. Material properties and environmental conditions are implicitly contained in the approach through the measurement data. The three key aspects pointed out in [31] in the context of bulk forming, reduce to two key aspects that also arise in AM for an efficient and effective implementation of the compensation procedure:


**Figure 2.** Flow chart of the compensation process for determining a suitable manufacturing geometry for AM in comparison to the forming tool compensation approach.

The problem of finding a suitable manufacturing geometry presents as an inverse problem under stochastic conditions. Hence, the design process is built in an iterative manner. A certain number of iterations *i* need to be passed, as shown in Figure 2. The manufacturing geometry *M<sup>i</sup>* leads to a workpiece geometry *W<sup>i</sup>* . With a predefined reference coordinate system and the target geometry *Wtarget*, the deviations *D<sup>i</sup>* can be derived. It is assumed that the deviation values *D<sup>i</sup>* are small compared to the component dimension. Hence, a non-rigid transformation *φW<sup>i</sup>* : *Wtarget* −→ *W<sup>i</sup>* between the target shape *Wtarget* and the workpiece shape *W<sup>i</sup>* can be deduced. Additionally the manufacturing geometry is supposed to be a non-rigid transformation *φM<sup>i</sup>* : *Wtarget* −→ *M<sup>i</sup>* of the the target geometry *Wtarget* . The idea of the geometrical compensation is to use the information of the non-rigid transformations to derive a new manufacturing geometry *Mi*+<sup>1</sup> in a recursive fashion by using

$$M\_{\bar{i}+1} = M\_{\bar{i}} - a\_{\bar{i}+1} (\mathcal{W}\_{\bar{i}} - \mathcal{W}\_{\text{target}}) = \phi\_{M\_{\bar{i}}}(\mathcal{W}\_{\text{target}}) - a\_{\bar{i}+1} (\phi\_{\mathcal{W}\_{\bar{i}}}(\mathcal{W}\_{\text{target}}) - \mathcal{W}\_{\text{target}}).\tag{1}$$

The factor *αi*+<sup>1</sup> weighs the compensation term. In its easiest form, *αi*+<sup>1</sup> is scalar, but in principle, depending on the topological structure of *Wtarget* , tensor-valued *αi*+<sup>1</sup> values are possible, leading to anisotropic compensations. Working with transformations of the target geometry is beneficial for computational reasons, since the same topology is ensured for straight forward calculations. In the present paper, the transformations are based on the target surface normals, which means that the measurement information can be directly used to build the transformations. The transformation *φW<sup>i</sup>* in this case reads

$$
\phi\_{\mathcal{W}\_l}(\mathcal{W}\_{target}) = \mathcal{W}\_{target} + D\_{\mathcal{i}} = \mathcal{W}\_{\mathcal{i}}.\tag{2}
$$

At the beginning of the compensation process, an initial manufacturing geometry has to be defined, either based on experience, numerical analysis or simply the target geometry. After the first run, the geometric adjustment and stochastic phenomena in turn influence the entire AM process again, in the sense of the inverse problem, which could call for another iteration to achieve the desired tolerances.

In this paper, the target geometry *Wtarget* is represented by a surface triangulation, as in common in AM. Based on the node coordinates and the connectivity list, a straightforward application of the transformations *φW<sup>i</sup>* and *φM<sup>i</sup>* is possible.

#### **5. Results and Discussion**

#### *5.1. Geometrical Deviation in the Calibration and Compensation Cycles*

In the calibration cycle (index 0), all 27 components in the build chamber were produced with the target geometry. In Figures 3 and 4, the first rows with index 0 show seven measured parts with their respective positions in the build chamber in Figure 1. Two main geometrical errors can be observed in the measured parts. There is a volume error due to the shrinkage of the material, and the upper right edge of the fin is distorted. It is evident that no optimal, universal compensation for all components is possible, due to the influence of the position in the build chamber. Hence, a unique compensation was calculated for each of the compensated parts *ai*–*e<sup>i</sup>* .

The results of the first compensation cycle are shown in the second row of Figure 3 with part numbers *a*1–*e*1. For part numbers *a*1–*d*1, the compensation reduced the geometrical deviation significantly. *e*<sup>1</sup> was not improved. Analysing the calibration cycle (*a*<sup>0</sup> − *e*0), *e*<sup>0</sup> already is an outlier in the area of the upper-right corner. The compensation was calculated to this deviation and did not improve the accuracy. After further investigation, a finishing error was detected. There was still loose powder on the part in this area, which emulated a thicker part to the measuring system. Hence, the compensation led to a thinner part in *e*1.

The second compensation cycle (index 2) showed minor geometrical improvements for *a*2–*c*2. For *d*2, a minor deviation increase was measured, while *e*<sup>2</sup> was improved significantly. This shows that this compensation process is quite robust regarding process outliers, since in the next compensation cycle they can be corrected without manual interaction. The limit of the compensation is reached when the statistical deviations are in the same order of magnitude as the deterministic deviations. The minor increase in deviation for *d*<sup>2</sup> shows that this limit is possibly reached after two compensation cycles. In Figure 4, the results of the reference components, which were not compensated, can be analysed. The parts show minor changes from one cycle to the next as well. This is due to the statistical deviations.

Comparing cycles 0 and 2, the compensated parts were significantly improved in terms of geometrical accuracy. Furthermore, the parts were more homogeneous and the influence of the position in the build chamber was diminished. In order to quantify this improvement and to analyse the limits of this method, the deviations were numerically analysed; see the next section.

#### *5.2. Evolution of the Part Deviations*

In Figure 5, the evolution of the mean absolute error and the evolution of the standard deviation are shown over the three iterations. For the reference parts (*b re f i* , *c re f i* and *d re f i* in Figure 1) the mean absolute error and the standard deviation remains almost constant, which is in accordance with the impression of Figure 4. For all compensated parts, an enhancement in the mean absolute error and the standard deviation was achieved. The process outlier *e<sup>i</sup>* is also clearly visible in the evolution of the mean absolute error, which stayed higher after cycle 1 but improved significantly in cycle 2. In Figure 5, another outlier with respect to standard deviation in cycle 1 can be seen, which is part *d<sup>i</sup>* . The high standard deviation in cycle 1 finally led to the increased mean absolute error of part *d*<sup>2</sup> in cycle 2.

*a*<sup>2</sup> *b*<sup>2</sup> *c*<sup>2</sup> *d*<sup>2</sup> *e*2 Deviation in mm

**Figure 3.** Compensated components for the calibration cycle (0) and two compensation cycles (1 and 2). The positions on the building platform of the parts *a*-*e* are detailed in Figure 1.

**Figure 4.** Reference components without compensation for the calibration cycle (0) and two compensation cycles (1 and 2). The positions on the building platform of the parts *b re f* -*d re f* are detailed in Figure 1.

**Figure 5.** Evolution of the mean absolute error and the standard deviation for the three cycles.

#### *5.3. Stochastic Deviations*

For the separation of systematic and stochastic errors, three reference parts *b re f i* , *c re f i* and *d re f i* are considered in each manufacturing job. In order to determine the stochastic portion of the procedure, the deviations are normalised with respect to the mean of the surface deviations according to Equation (3) for each measurement point.

$$
\mathfrak{x}\_{normalized} = \mathfrak{x}\_{\mathrm{i}} - \overline{\mathfrak{x}}.\tag{3}
$$

Based on the normalised values, the stochastic mean absolute error and the stochastic standard deviation can be calculated and are summarised in Table 1. The stochastic error measures give the values for the whole procedure shown in Figure 2, including manufacturing (or numerical analysis), measurement (evaluation) and the redesign process of the manufacturing geometry.

The average values of the stochastic errors over all reference parts and all iterations are suitable estimates to characterise the stochastic behaviour of the procedure. With regard to the approach, which is capable of dealing with systematic deviations only, these values could be seen as a process limit for this specific geometry and job in this configuration using the procedure described.

**Table 1.** Separated stochastic errors of the compensation procedure conducted for the three reference parts. (mae: mean absolute error; std: standard deviation).


#### **6. Conclusions**

In this study, a systematic compensation approach from forming technology was adapted and applied for AM. The Am manufacturing geometry was successfully used as a forming tool in the sense of geometrical compensation. The displacement adjustment method, which directly incorporates deviation data for compensation, presented itself as suitable for AM. Compared to conventional manufacturing technologies, the presented procedure is even more reasonable for AM processes, since no expensive physical forming tools or moulds are necessary.

The presented approach is a data-based method, which relies on an initial calibration job. We think that with regard to high engineering and computational costs of AM numerical analysis, this approach is reasonable, especially with regard to build chamber dependencies. The necessary effort and experiments for a highly complex material model have to be taken into account as well.

By using non-rigid transformations of the initial geometry topology, a direct application of measurement data for the compensation is realised. This way, each part is treated individually, which ensures a compensation of spatial dependencies in the build chamber. The systematic procedure is suitable for series production especially, where the first iteration can be seen as a kind of calibration step. We conclude that the proposed hypothesis could be proven for the setup conducted. Of course, future work has to be done to further verify the approach.

During production, the approach can also be used for process control, in the sense of a closed-loop control, and, for example, to compensate for other raw materials, system characteristics or changed environmental conditions. In future work, the approach will also be tested for other AM processes and different materials. Furthermore, the incorporation of simulation results for the initial geometry is a promising way to further enhance the procedure.

**Author Contributions:** Conceptualisation, C.H., P.L. and B.H.; methodology, C.H., B.H. and P.L.; investigation, C.H., P.L. and Y.K.; writing—original draft preparation, C.H. and P.L.; writing—review and editing, B.H., T.C.L. and W.V.; supervision, T.C.L. and W.V.

**Funding:** This research was funded by the Federal Ministry of Education and Research of Germany (BMBF) under grant number: 13N15085, Industrialisierung und Digitalisierung von Additive Manufacturing (AM) für automobile Serienprozesse, acronym IDAM.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Analysis, Optimization, and Characterization of Magnetic Photonic Crystal Structures and Thin-Film Material Layers**

## **Mikhail Vasiliev \* , Kamal Alameh and Mohammad Nur-E-Alam**

Electron Science Research Institute, School of Science, Edith Cowan University, 270 Joondalup Drive, Joondalup, WA 6027, Australia

**\*** Correspondence: m.vasiliev@ecu.edu.au; Tel.: +61-8-6304-5809

Received: 20 May 2019; Accepted: 4 July 2019; Published: 5 July 2019

**Abstract:** The development of magnetic photonic crystals (MPC) has been a rapidly evolving research area since the late 1990s. Magneto-optic (MO) materials and the techniques for their characterization have also continually undergone functional and property-related improvements. MPC optimization is a feature-rich Windows software application designed to enable researchers to analyze the optical and magneto-optical spectral properties of multilayers containing gyrotropic constituents. We report on a set of computational algorithms which aim to optimize the design and the optical or magneto-optical spectral analysis of 1D MPC, together with a Windows software implementation. Relevant material property datasets (e.g., the spectral dispersion data for the refractive index, absorption, and gyration) of several important optical and MO materials are included, enabling easy reproduction of the previously published results from the field of MPC-based Faraday rotator development, and an effective demonstration-quality introduction of future users to the multiple features of this package. We also report on the methods and algorithms used to obtain the absorption coefficient spectral dispersion datasets for new materials, where the film thickness, transmission spectrum, and refractive index dispersion function are known.

**Keywords:** 1D magnetic photonic crystals; multilayer film modeling; modeling of Faraday rotation spectra; MPC optimization; exhaustive computation; materials characterization

#### **1. Introduction and Background**

In recent years, there has been some resurgence in the research interest dedicated to engineering and characterization of magneto-optic iron garnet materials [1–4]. Thin-film magnetic garnets are semitransparent magnetic dielectrics possessing record specific Faraday rotations of up to several thousand ◦ /cm, in the near-infrared spectral range, if containing bismuth substitution [5]. Generically, the chemical composition of garnet materials of this type is described by the formula (BixRe3-x)FeyM5-yO12, where Re is rare-earth metal (e.g., Dy, Sm, Lu, Nd, or Ce), and M is transition metal such as Ga or Al [5]. The exploration of this important subclass of functional materials has a decades-long history, starting from the days of bubble-domain magnetic memory development, and more recently, continued with renewed research activities, in application areas ranging from on-chip nonreciprocal components (waveguide isolators, [2,6]), to magnetic recording [7]. Magnetic garnet materials synthesized by a range of thin-film deposition techniques and crystal-growth methods have also attracted a significant research and development momentum since the 1990s, and throughout the last two decades, in areas ranging from photonic crystals to spintronics [8–15].

Various approaches to the design and manufacture of magnetic photonic crystals (MPC) with tunable properties, and potentially suitable for the fabrication of nonreciprocal optical components have been reported [8–18]. The one-dimensional MPC are structured (periodic or quasi-periodic) sequences of the magnetic and nonmagnetic material layers deposited on optical substrates, which are defined by the layer order and thickness. Many research groups have focused on optimizing the thickness of one-dimensional (1D) MPC structures to simultaneously achieve 45◦ of Faraday rotation angle and maximal possible transmission at optical telecommunication wavelengths. The MPC based on quarter-wavelength thin-film stacks, which are sequences of magnetic and nonmagnetic layers with multiple embedded phase shifts (termed defects, or missing layers), have been shown to possess a significant potential for practical implementation of integrated optical isolators. This is due to the necessity of approaching Faraday rotations as large as 45◦ , which has been shown to be attainable, due to the resonant enhancement of Faraday rotation observed in such structures. Complex 1D MPC designs featuring a "flat-top" spectral response, with almost 100% transmission within a large bandwidth (several nanometers), and close to 45◦ of Faraday rotation at 1550 nm can contain in excess of 200 layers [10,11], limiting their applications to the near-infrared range, where the magnetic garnets possess very low (negligible) optical absorption. Since the original development of the MPC optimization program in 2005 [18], multiple garnet material development efforts have been undertaken within our group [19–21], all of which have relied substantially on material characterization options featured within the same software package. In particular, the option of deriving the data for the spectral dispersion of the absorption coefficient using the measured transmission spectrum and refractive index dispersion data has been very useful to characterize new nanocomposite-type garnet materials synthesized by cosputtering using an extra oxide material source [19]. The present-day performance limits of 1D MPC in the visible spectral range have been evaluated using the same software [22] and using the available spectral data on the optical properties of the best-performing magnetic garnet compositions. Optical constants data of multiple recently-synthesized magnetic garnet compositions have been evaluated using the measured transmission spectra in conjunction with MPC Optimization software and Swanepoel envelope method [23]. Our group's preferred method for calibrating the quartz microbalance sensor's tooling factors of various deposition sources also involves fitting the actual film layer thickness using MPC Optimization software in conjunction with measured transmission spectra.

A graphical snapshot of the controls and features implemented within MPC Optimization software is shown in Figure 1, which also presents a sample optimization result. The computation time to obtain this result is typically around one min (after analyzing almost 5000 MPC designs out of a possible 10,000 defined by the preset film structure features). Some of calculated designs exceed the maximum limits set for thickness or layer number, and thus are not fully evaluated. The default spectral range of calculations is extremely broad, slowing down the optimization algorithm significantly, since the default wavelength settings and resolution have been entered to assist in fitting film layer thickness conveniently, which is one of the most common everyday applications of this software. Running the actual MPC optimization algorithm is best performed by zooming onto the spectral range of interest, which is usually represented by a wavelength region surrounding a narrow peak of transmission or reflection, where it is possible to engineer the enhancement of Faraday rotation.

The default-entered materials-related data used to optimize MPC designs (within the parameter space and constraints also entered as default values not relevant to any particular application) relate to a MPC design which is based on magnetic garnet composition Bi2Dy1Fe4Ga1O<sup>12</sup> and SiO<sup>2</sup> L-type layers, deposited on a glass substrate. The optical constants and gyration data for wavelengths near 630 nm have been obtained from thick garnet layers of this composition, and the corresponding index dispersion dataset is also preloaded into the menu item "Extra data | Account for refractive index dispersion".

**Figure 1.** Front-panel controls of MPC Optimization software and a sample optimized MPC design obtained by clicking the "optimize" button without changing the default-entered data. The result is a MPC structure which is designed to operate at 630 nm, reliant on a magnetic garne<sup>t</sup> material of specific Faraday rotation near 2◦/µm at that wavelength, however, the structure enhances the Faraday rotation to 4.98◦/µm within the spectral transmission peak.

The contents of the compiled html (.chm) help file accessible from the "Help" menu are sufficient to enable most beginner-level MPC (or thin-film) designers to quickly learn the main features of program and its data representation formats. In the following sections of this article, these main features are described in detail, with examples provided to enable the productive and convenient use of this feature-rich software package. The main aim of this present work has been to provide a set of computational tools and algorithms for use by the developers of MPC and other nanostructured material systems (e.g., thin-films and multilayers containing MO materials). These tools will enable both the characterization of functional material layers and the application-specific design of Faraday rotators. Other contributions of this work include a number of experimentally-validated optical and MO property datasets of several garnet material types possessing giant Faraday rotation and good MO quality, synthesized by our group using RF sputtering followed by high-temperature annealing crystallization processes. Examples of particularly important MPC-based Faraday rotator design types are reviewed, focusing on the ways that strong peaks of Faraday rotation are engineered to spectrally coincide with the peaks of transmission, while paying special attention to the role of limiting factors such as absorption.

#### **2. Overview of Package Operation and Key Examples**

Since the structure of 1D MPC is essentially represented by a sequence of magnetic and nonmagnetic material layers on a substrate, most of the core terminology, design structure description conventions, and analysis techniques are derived from the field of multilayer thin-film design. The input datasets necessary to define the layer sequence and the optical properties of each material type and individual layer are entered into the relevant text boxes, starting from the top-left corner of the Windows Form. The design formula field defines the layer sequence, starting necessarily from the capital letter S defining the substrate. The layer sequence must be entered in an alphanumeric string format containing special symbols such as round, square, or curled brackets, corresponding to one of the three main design-string representation types. These types are: (i) quarter-wave stack-based notation, e.g., SLH(ML)2 which is perhaps the most common system of layer structure abbreviation in thin-film design; (ii) physical thickness-based notation, e.g., SM [1000]L[5 0], and (iii) advanced designs notation, e.g., S{1.0}(L/2HL/2)1{1.06}(L/2HL/2)3 further described within the help documentation. The optical properties of the substrate material (which is presumed to be dielectric, semi-infinite, and non-absorbing) are defined using only the real part of refractive index. Up to three different optical materials are allowed, denoted by the letters H, L, and M, however, H-type layers are not restricted to mean "high refractive index". The layer-naming conventions are derived from thin-film terminology for convenience only. In the preceding thin-film abbreviation example (SLH(ML)2), the symbols mean a substrate (S) with the following sequence of layers deposited, in the following order: L, H, M, L, M, and L.

The M layers can optionally be modeled as magnetic dielectrics, where a gyration value at the design wavelength (the imaginary part of the nondiagonal dielectric tensor component, or Im(εxy)) needs to be entered into its relevant text box. This imaginary component of the off-diagonal dielectric tensor element describes magnetic circular birefringence and manifestes as Faraday rotation of the polarization plane in the transmitted or reflected light waves. The real part of this tensor component is treated as zero, in most applications, unless the experimentally-measured value is known, which is related to characterizing magnetic circular dichroism and polarization state ellipticity effects. For applications requiring good numerical accuracy over broad spectral ranges, the spectral dispersion of both the gyration and the refractive indices of all relevant materials need to be loaded from .txt format data files, using the options within the "Extra Data" menu. Optionally, H-type layers are selected to also represent a magnetic dielectric, with the same off-diagonal dielectric tensor components as in M-type material (but optionally with different refractive index and absorption), in order to model a special physical situation where an MPC has layer-specific magnetization reversal possibilities in some individual magnetic layers within the structure. If M layers are modeled as nonmagnetic (e.g., just implying medium-index dielectric layers), then both the real and imaginary parts of εxy are entered as zero values, however, if a small gyration is still entered by error, it will not measurably distort the transmission or reflection spectral calculations, and therefore, in this case, the results regarding the Faraday rotation spectra should be ignored.

The details of the physical situation being modeled, in terms of incidence geometry, and the ways in which the transmitted (or reflected) light intensity is normalized with respect to the incident wave intensity, are defined using sets of checkboxes within the menu entry "Incidence Geometry", where the relevant descriptions are given. The incidence geometry settings defined as default are the ones used most often and generic in general, since these enable the convenient and accurate fitting of lab-measured transmission spectra in deposited film-substrate systems, to the corresponding theoretically-modeled spectra, where the effects of each interface (including the back surface of substrate) are correctly accounted for in the model. The only parameter not accounted for is the physical thickness of substrate, which is modeled as non-absorbing. Alternative settings for the incidence geometry are useful for considering more theoretical situations, such as calculations of the optical intensity transmitted into a semi-infinite substrate medium, or in reflection-mode calculations, where it is often necessary to compare reflected-wave intensities, which vary with the direction of incidence.

One of the most important material parameters in all layer types is the optical absorption coefficient at the design wavelength (entered in cm−<sup>1</sup> into relevant text boxes; the corresponding extinction coefficients will then be displayed after the film design is characterized by pressing the "compute formula" button), especially for materials possessing significant spectral dependency in their optical absorption. For accurate characterization of thin films or MPC, ideally, every material should have its optical constants dataset available for loading from the "Extra Data" menu option "Account for refractive index dispersion" and loaded into the specialized form (shown in Figure 2) prior to calculations. The material-specific optical constants data files are prepared using zeros entered in place of an unknown absorption coefficient, for the purposes of physical layer thickness fitting, based on the accurately measured transmission spectrum data (as discussed more in detail in subsequent sections). The end-of-file (EOF) marker in these .txt material data files prepared using editor applications such as Notepad must be placed immediately after the last figure in the last column, by way of making sure to delete any possible symbols or empty row spaces below.

Once the refractive index dispersion information is loaded from data file(s), the data in text boxes that corresponds to the index and absorption at the design wavelength are no longer used in the main spectral calculations, but used only for calculating the quarter-wave stack physical thickness in nm. For all wavelength points between the 27-point data grid, the values of refractive index and absorption are linearly interpolated from the nearest grid-located points. This can erroneously cause small spectral shifts to appear in the spectral locations of the transmission, reflection, or Faraday rotation peaks, seen away from the precise design wavelength, if the interpolated refractive index at that wavelength does not coincide with the value entered into front-panel text box. This is not a significant issue for the experienced designers of MPC, once the origin of these possible small data errors is known or eliminated by entering precise (same as interpolated from the index dataset for the design wavelength) index data into front-panel text boxes. It is a known a priori that the actual spectral response peak locations in quarter-wave stack-type designs will appear precisely at the design wavelengths, due to the nature of optical interference-related phenomena. In situations when the refractive index dispersion data are only available within a limited spectral interval of interest, rather than for all wavelengths in the 350–1600 nm range, it is recommended that the available refractive index data are entered into the newly-generated data file. The refractive index and absorption coefficient values at other wavelengths still need to be entered. The recommended practice is as follows: for example, if the available data starts from 500 nm, enter the same values as are known at 500 nm into the wavelength grid positions for all shorter wavelengths; alternatively, use any available models to predict the unknown values (e.g., Cauchy dispersion model). If the refractive indices or absorption are only available up to 800 nm, it is best to enter (for all larger wavelengths), the same (n, A) data as at the last data point (800 nm).

Whilst the optical properties cannot be extrapolated outside a known range by the values measured at the range limits, the designers of thin film and MPC structures are encouraged to make numerical experiments with the package to evaluate the errors expected to result from the possible refractive index data inaccuracies, in each particular design case. The refractive index of the exit medium surrounding the substrate-film system is defined only by its real part (typical value is 1.0 for air), and the same exit medium is presumed to precede the (thickness-undefined) substrate, and to exist beyond the last deposited film layer, regardless of the direction of light incidence. The checkbox "reflection mode" sets up the calculations of the reflection spectra, and also the Faraday rotation spectrum for the reflected light, if checked. Optionally, the "show absorption spectrum" checkbox is checked, after which the "compute formula" button will initiate new calculations, resulting in the display of the calculated absorption spectrum.

**Figure 2.** Form dedicated to loading the optical constants data from text files prepared as shown, using the software-specific 27-row wavelength grid and containing the columns data for the refractive index and optical absorption coefficient (in cm−<sup>1</sup> ) at each wavelength point within the spectral grid.

The theoretical analysis of the spectral properties of MPC has been implemented using the 4 × 4 complex-valued transfer matrix approach [16]. The method is reliant on finding a transfer matrix that relates the complex electric and magnetic field amplitudes of the light in front of the MPC and behind it. The transfer matrix for the whole structure is represented by the product of the transfer matrices calculated for each layer of the MPC. The transfer matrices of layers are determined by the thickness of the layers and their dielectric tensor components. Nonmagnetic layers are treated as isotropic media having a diagonal dielectric tensor, whereas, magnetic layers are described by tensors containing the magnetization-dependent off-diagonal dielectric tensor components containing wavelength-dependent gyration information, which are related to specific Faraday rotation spectra. When running the optimization algorithm, a C++ code implemented using Microsoft Visual Studio.NET 2003 generates a look-up table of all possible design structure variations (within the design constraints specified). From this table, an optimum MPC structure is obtained by successive tightening of the desired spectral response specifications, guided by the user inputs. The optimization approach used in our algorithm is

based on the "exhaustive computation" of the entire parameter space defined by an arbitrary selected type of a structural formula composed of "design substrings" representing user-selected elementary building blocks of the MPC design.

#### *2.1. Multi-Defect Multilayer MPC Characterization and Optimisation Examples Suitable for Validating Calculations*

To illustrate the suitability of software to correctly calculate the spectral responses of complex, multi-defect MPC designs, it is easiest to use the design or optimization examples described in the published literature, e.g., [11] and [22]. One of the common goals of optimizing the MPC structures has been to achieve strong spectral peaks in either the transmission or reflection (ideally approaching 100%) that coincident spectrally with peaks of enhanced Faraday rotation, in either the transmission or reflection-mode operation, and ideally approach 45◦ , if efficient modulation of light intensity is required. It is important to note that Faraday rotation in the reflected wave is different in its physical nature from Kerr rotation [5] (although there may exist some terminological misinterpretations, even in the published literature). This program calculates only the angles of polarization-plane rotation due to Faraday effect, in either geometry, and does not account for Kerr effect. Figure 3 presents a graphical summary of the input parameters necessary to enter in the relevant text boxes to evaluate one remarkable MPC design from the published literature [11], as well as the results of calculations. The practical implementation of this MPC design is expected to be difficult, due to factors such as the extremely high number of layers, large total thickness, the expected nonzero (but possibly well below about 0.1 cm−<sup>1</sup> ) optical absorption coefficient in garnets at 1550 nm, as well as the scattering effects expected to occur at multiple layer boundaries. From the theoretical insight perspective, this high-performance MPC design is still an outstanding example of MPC application potential, especially in systems using optimized surrounding media, index-matched to the mean refractive index of structure. A .mpc file (MPC Design from JLT Vol. 19 No. 12 p. 1964.mpc) is included in the subfolder "Optimization results files" of the program installation directory, and can be loaded from the "File" menu, enabling easy recalculation of the contents of Figure 3, using preloaded design data. The way this MPC has been modeled also involves the customized incidence geometry settings, where film-side incidence is modeled, without normalizing the transmitted intensity after the substrate. Rather, transmission into the substrate material is modeled, which explains why the modeled transmission within 1550 nm peak closely approaches 100%. Running the MPC optimization algorithm, on the other hand, requires using substrate-side incidence geometry settings only.

All materials-related data values were used as per the description in [11], and the data in the calculated graphs reproduce the results presented in Figure 3b of [11] with precision. Due to materials-related constraints, such as the spectral dependencies of the absorption coefficient and gyration, achieving strong enhancements in Faraday rotation simultaneously with low optical losses becomes progressively more difficult with the reducing design wavelength. Across the visible spectral range, the optical absorption in bismuth-substituted iron garnet materials becomes the limiting factor, placing stringent limits on the achievable MPC performance characteristics, regardless of their intended application area or the design type. Therefore, the ability to generate and compare multiple and differently-structured optimized MPC designs is crucial for achieving the best possible performance characteristics, limited only by the fundamental, materials-related constraints.

In order to directly reproduce the optimization result reported in [11], by way of running a constrained optimization algorithm using the MPC Optimization software, the materials-related datasets and a set of optimization constraints, as shown in Figure 4, must be used prior to clicking the "optimize" button. It is also necessary to set the substrate-side incidence geometry, and uncheck the second checkbox related to the way the transmitted intensity is normalized. The "flat-top" optimized MPC design will be retrieved from the database of calculated MPC designs within the set of defined criteria (60 different designs will be found to fit these overall criteria and constraints, as set per data of Figure 4), by selecting the design with the maximized spectral response bandwidth. This is done by pressing the "MAX (B)" button. Importantly, all design substrings must be entered as per Figure 4 to reproduce this optimization result, with precisely five substrings and a maximum of 30 entered for the substring-repetition index. Another constraint which needs to be entered relates to the maximum allowed layer number which is 203. The checkbox "calculate symmetric designs only" needs to be checked, and it is best to set the spectral resolution to 1 nm during optimization, followed by changing it to 0.1 nm for calculating the spectral properties more accurately. The optimized design equivalent to the design of Figure 3 (and [11]) is retrieved from 60 possible designs found within the criteria, by pressing the "MAX (B)" button. This retrieves the design with maximized full fidth at half maximum (FWHM) bandwidth, according to the data entered into the TMAX(%) text box.

**Figure 3.** MPC optimization software used to reproduce the flat-top MPC transmission and Faraday rotation spectral properties for a four-defect design reported in [11]. The calculated graphs (obtained after entering the design data and pressing the "compute formula" button) reproduce the data originally presented graphically in Figure 3b of [11].

The settings applied to the normalization of the transmitted intensity within the Incidence Geometry menu corresponded to the physical situation equivalent to that applied during the calculation of the MPC properties in [11], stipulating the transmitted intensity normalization procedure involving transmission "from within" the substrate material, into the index-matched exit medium. Figure 5 shows a graphic summary of the MPC optimization results first reported in [22], which illustrate the current performance limits of MPC designs aimed at developing transmission-mode magneto-optic light intensity modulators working at near 650 nm. The optical and MO material parameters relevant to Bi2Dy1Fe4Ga1O<sup>12</sup> and similar highly bismuth-substituted nanocrystalline garnet materials (synthesized by techniques such as RF sputtering) were used in the calculations.

The data of Figure 5 can be reproduced by running optimization of four-defect structures with up to 40 layers and thickness up to 5 µm, composed of five sequenced (LM) and (ML) substrings, as shown, using the default entered n(L) value and n(M) = 2.376. The gyration value corresponding to 2◦ /µm needs to be entered as –0.02 for wavelengths near 650 nm, accounting for the composition-specific sign of Faraday rotation, which is defined by the convention reported in [5] and other sources. The maximum repetition index value is set to either 10 or, e.g., 12, prior to running the optimization with either absorption coefficient.

**Figure 4.** Optimized "flat-top" high-performance MPC design reproduced by running the MPC optimization algorithm, using checkbox "Calculate symmetric designs only" for speeding up calculations. The exact four-defect MPC design reported in [11] (Figure 3b of [11]) is shown after selecting the design with maximum spectral response bandwidth from the 60 possible MPC designs found to be within the optimization criteria and constrains, as shown.

**Figure 5.** Calculated spectral performance parameters for multi-defect (four-defect structures, having up to 40 total layers) MPC optimized by exhaustive computation to achieve a maximum light intensity modulation capability (maximized value of parameter (T·sin 2 (2**Φ**MPC)), for n(M) = 2.376, n(L) = 1.458, and using two different optical absorption coefficients for magnetic material at 650 nm (a) α = 1000 cm−<sup>1</sup> and (b) α = 2000 cm−<sup>1</sup> . These refractive index and absorption coefficients were considered constant within the design-specific wavelength region, as well as gyration (–0.02), corresponding to approximately 2 ◦ /µm near 650 nm. The graphical information is reproduced from [22].

The results shown in Figure 5 illustrate clearly that the optical absorption is the limiting factor in the visible-range MPC design, even at long-wavelength red wavelengths, where the absorption is still moderate, and the thin constituent garnet layers used within MPC would have been almost visually clear. This is demonstrated by entering a design string, such as SM[68.39], into design formula text box, and running the absorption-mode calculation (e.g., using α = 2000 cm−<sup>1</sup> ) to reveal the graph area by zooming in with mouse, that individual MO layers absorb only about 1.2% of the incident power on each single-pass transmission).

#### *2.2. Generating Optimized Antireflector Film Designs Using Spectral Target Points*

It is possible to apply additional optimization constraints at up to three selected wavelength points, to force the algorithm to output the designs with specific spectral features, in either the transmission or reflection mode. An example of obtaining the optimized five-layer thin-film broadband antireflector coating designs is stored in the file Try\_optimise\_5-layer antireflector.mpc which is placed into the subdirectory of "Optimization results files" in the installation directory. The menu item "Extra Data | Multi-wavelength spectral targets" is used to enter the additional optimization constraints, regardless of whether the Faraday rotation features are being optimized or not. Figure 6 shows the required inputs within the two submenus related to the multi-wavelength spectral targets and the incidence geometry options, which will result in reproducing the five-layer antireflector design shown also in Figure 6. Selecting the design with maximum transmission at the design wavelength (either after running the optimization, or simply after loading the relevant example file) will reveal the reflection spectrum as shown.

This example also demonstrates the use of scaled QWOT layers for use in thin substrings, the thickness of which is then being optimized by the algorithm by adjusting the repetition indices. The calculation of more than 750,000 designs should still take only a few minutes. A number of antireflector-type designs are revealed by using the button "display design N=" with any corresponding number not exceeding the number of designs found within criteria.

Since the optimization algorithm presumes the substrate-side incidence, it is convenient to define air as the subtrate, and set the exit medium to glass. The obtained design S(LL)10(HH)14(LL)4(HH)1(LL)14 needs to have its deposition sequence reversed and be re-evaluated for the film-side incidence case from air. Using the physical thickness notation is preferable in this case, i.e., the design needs to be changed into S(LL)14(HH)1(LL)4(HH)14(LL)10. In physical thickness notation, for a practical deposition-ready design description, this is equivalent to SL[139.49]H[8.19]L[39.85]H[114.68]L[99.64]. It is important to not forget to set n(S) to 1.5 and n(exit) to 1.0(air) in this case. Now, the incidence geometry settings are checked to correspond to the film-side incidence. Note, the way the reflectivity of the back side of the substrate is accounted for in the detailed incidence geometry settings.

#### *2.3. Fitting of the Measured Thin-Film Transmission Spectra to Theoretical Models*

One of the most frequently used applications of MPC Optimization software, apart from optimizing the MPC structures, is expected to be the fitting of actual deposited layer thickness, for thin-film materials with known refractive index dispersion function. The option of loading the measured transmission spectrum for immediate comparison with the modeled transmission spectrum of the same substrate-film system is available from the menu item "Inverse problem | Load T spectrum for fitting". To enable accurate modeling, both of the two checkboxes that correspond to the substrate-side incidence within the menu "Incidence Geometry | Define geometry" must be checked, which is done by default. A measured reflection spectrum (if available for the normal-incidence reflection, which is rare with most instruments) can also be fitted in the same way as transmission, by running the calculations in reflection mode, with the corresponding checkbox checked. During the calibration of in-situ thickness control systems, e.g., quartz crystal microbalances or reflectometer-type systems, the task of determining the actual physical thickness of deposited thin-film layers is very common, where MPC Optimization software can be used effectively, in conjunction with other methods such as SEM or profilometry characterization.


**Figure 6.** *Cont.*

**Figure 6.** Example of menu and algorithm settings required to generate a number of optimized five-layer antireflector film designs on a glass substrate. In this example, the optical materials (ZnS and MgF<sup>2</sup> ) are presumed to possess constant refractive indices and zero absorption across the entire visible spectrum.

Figure 6 illustrates graphically the results of fitting the loaded (measured) transmission spectra of two thin films of composition type Bi0.9Lu1.85Y0.25Fe4.0Ga1O<sup>12</sup> (the refractive index dispersion data file for this composition is supplied within the appropriate sub-folder in the program installation directory of the target machines). A thinner (684 nm) film modeled as deposited on a glass substrate (n(S) = 1.5) was fitted using the refractive index and absorption coefficient data file related to the as-deposited (amorphous-phase) garnet films of this composition. Since the transmission of an annealed (garnet-phase) film was actually loaded, its transmission at shorter wavelengths was in excess of that modeled. The absorption of garnet-phase films is practically always less in the crystallized state as compared with amorphous garnet-precursor films. A thicker film of fitted thickness (transmission shown in Figure 7b) was measured in the amorphous phase, and therefore the quality of fit is better. Some systematic transmission discrepancies across a wide spectral range can be attributed to a combination of possible light scattering on the film surface features and film layer refractive-index nonuniformity.

If the spectral dependency of the film material absorption coefficient is completely unknown (or the data are not reliable), but the dispersion of its refractive index is well known (e.g., from variable-angle spectroscopic ellipsometry data), then the index dispersion data files need to be prepared with zero values entered for all absorption coefficients. These data files still enable, in many cases, very reliable fitting of the physical thickness. It is highly recommended, then, to back up these thickness fitting results by also applying the Swanepoel method-based techniques, e.g., methods reported in [23].

**Figure 7.** Magneto-photonic crystal (MPC) software fitted transmission spectra of different thin film garnet layers. (**a**) Annealed garnet sample of composition type Bi0.9Lu1.85Y0.25Fe4.0Ga1O<sup>12</sup> and thickness 684 nm and (**b**) as-deposited garnet sample of the same composition but from another batch and thickness 1310 nm.

#### *2.4. Fitting of the Absorption Coe*ffi*cient Spectra in Single-Layer Films of Known Refractive Index Dispersion Function, Transmission Spectrum, and Thickness*

In situations when only the refractive index spectral distribution and the measured transmission spectrum of a semitransparent material layer are available, it is possible to use the custom-prepared "zero absorption" refractive index dispersion data file, and then first fit the physical thickness (typical results are shown in Figure 8a), followed by a derivation of fitted absorption coefficient spectrum using a dedicated algorithm from the program menu item "Inverse problem | Derive absorption spectrum". In order to minimize errors, this combination of fitting procedures is only recommended if several deposited films of the same material are available which have a significantly different physical thickness, and the fitting procedures applied to more than one film consistently lead to obtaining well correlated datasets as a result of the absorption fitting. It is also recommended to apply the Swaneopol method-based fitting procedures [23], to reconfirm the validity and accuracy of the physical thickness data.

**Figure 8.** Iterative (bisection algorithm-assisted) fitting of the absorption coefficient spectral dependency and the required pre-fitting of film thickness through matching transmission peak features. (**a**) Peak-aligned transmission spectrum pre-fitting result, from which the physical film thickness information and the measured dataset on the refractive index dispersion are then used, within the option available in the "Inverse Problem" menu, to derive the absorption coefficient data; (**b**) the algorithm-derived (fitted) absorption coefficient spectrum (α, cm−<sup>1</sup> ) for 481 nm thick as-deposited Bi0.9Lu1.85Y0.25Fe4.0Ga1O<sup>12</sup> garnet-precursor film sample on a glass substrate (its measured transmission spectrum is shown as blue curve in a).

Figure 8 illustrates the results of physical thickness and absorption fitting procedures obtained using a transmission file of a 481 nm thick Bi0.9Lu1.85Y0.25Fe4.0Ga1O<sup>12</sup> as-deposited film on a glass substrate, and the corresponding material's "zero-absorption" refractive index data file, both of which are included with program distribution in the corresponding subdirectories of the installation folder. The data file named 701 nm Bi0.9Lu1.85Y0.25Fe4Ga1O<sup>12</sup> garnet layer\_n (SWEM) at A = 0 .txt must be loaded using the "load M data" button and the subfolder corresponding to this material type (we previously characterized the refractive index of a 701 nm thick film of this material). After the best fitted thickness value (481 nm) is obtained by comparing different models of design string, such as SM[450], . . . , SM [481], the same design string must be re-entered using a quarter-wave thickness notation, e.g., S(MM)3, where the QWOT multiplier for M thickness is set to 1.125 (for n(M) = 2.21, corresponding to the nearest index data point to the default 630 nm), to match the physical thickness of 481 nm in this notation. The text box "structure thickness (microns)" must be used to check the physical thickness changes in response to changing either the (MM)N repetition index, or the QWOT multiplier for the M layers. The index n(M) = 2.21 should be entered into its corresponding text box, after looking up the 630 nm (default design wavelength) data for the actual refractive index, to avoid possible misrepresentations of the QWOT data. This notation conversion is required for running the absorption coefficient fitting algorithm, as well as using M-type layers only, regardless of whether the material possesses any magnetic properties or not.

Once within the "derive absorption spectrum" submenu, the same measured transmission file should again be reloaded on the graph from file, followed by a relatively self-explanatory procedure for deriving the graph of absorption spectrum. No changes are usually required in any other text boxes. The result of the fitting procedure with the data files (as described above) is shown in Figure 8b.

In the cases where the fitting algorithm produces exceptions such as described within a dialog text box, e.g., "at some or all wavelengths, even at zero absorption in M layers, the transmission of this structure should be less than specified", adjustments are made to either check that the incidence geometry settings are correct, or the refractive index of the substrate is increased in the model, removing these fitting procedure errors. In cases when these issues persist in the low-absorption wavelength ranges, the wavelength range of the model is reduced to include only the regions where the fit results can be obtained. Nonuniformities in real thin films can lead to reduced refractive index values, leading then to reductions in the modeled reflection, thus showing increased transmission at some wavelengths as compared with theory models. After the absorption spectrum fitting procedure has been completed, the plotted data points are exported into other formats or saved in the data files using the options provided.

#### **3. Installation and Code-Related Information**

The installation of MPC Optimization software is easy, and it is enabled by running the installer (.msi) file supplied within the .zip archived folder used for program redistribution. A necessary prerequisite to program installation is the Microsoft .NET Framework 1.1, which must be installed on any Windows machine prior to running the MPC Optimization installer. The .NET 1.1 Framework installation file (dotnetfx.exe) is supplied within the archived folder file used to redistribute MPC optimization. So far, no known problems have been identified to exist in relation to installing this older version of .NET Framework on modern computer systems. A number of example data files containing the samples of measured garnet thin film transmission spectra and files containing the data on the refractive index and absorption coefficient spectra of various garnets and other optical materials are placed into the selected program installation folder on installation. These data enable the users to recreate the example calculations presented within this article, and therefore are useful in mastering the software operation.

The program has been written as Microsoft Visual Studio .NET 2003 Solution and it is composed of three projects: (i) class library project written using managed C++ code, and built as .dll class library providing matrixrelated computation functions; (ii) main program project, written using visual C# and implemented as Windows forms project; and (iii) deployment project, used for including the necessary reference assemblies, class libraries, license files, organizing the file system structure within the installation directory, and generating the Installer files for program distribution. The installed MPC optimization program runs on any Windows system architecture, whether 32-bit or 64-bit, regardless of the processor type.

Two main third-party software components have been embedded in the program structure and the necessary developer licenses have been purchased commercially. Bluebit Matrix Library 2.2 has been used to provide classes that enable efficient complex-valued 4 × 4 matrix operations functionality. The class library .dll assembly referencing the Bluebit Matrix Library assembly has been required, through the terms of a developer-grade matrix software license, to be compiled as a signed assembly, using a strong-name assembly key file, which is not allowed for redistribution by the developers. The second embedded component was GigaSoft ProEssentials 5 .NET package, which has been installed on the developer's computer to provide ActiveX controls that enable scientific-type graphing output functions. The assemblies GigaSoft.ProEssentials.dll and Pegrpcl.dll are referenced by the C# project, embedding the scientific graphing controls as Microsoft component object model (COM) objects into the structure of Windows forms project-related software assemblies. The use of licensed third-party software components, with the associated restrictions, is the primary reason why the MPC optimization program code is not intended for open-source distribution. Additionally, newer versions

of Bluebit Matrix Library designed to work with present-day versions of .NET Framework, as well as significant code syntax changes applied throughout the projects, would have been required to successfully port the solution into current versions of Visual Studio, e.g., VS2015 or VS2017. The program installation pre-requisites, installation files, and related documentation are available for download from the Supplementary Material section. Appendix A provides a concise summary of the program functionality, code-related technical details, and outlines some limitations in relation to using the embedded graphing controls.

#### **4. Conclusions**

In summary, a set of computational approaches, and a custom software package have been described, designed to enable the design and optimization of 1D magnetic photonic crystals in terms of the achievable combinations of Faraday rotation, transmission, and reflection spectra. A number of important MPC design examples have been introduced, illustrating both the desirable optical property combinations, and the materials-related performance limits of MPC-based Faraday rotators. The same package allows computational modeling of the optical spectral properties of various dielectrics-based generic single- and multilayer thin films. Additional program features include the tools for fitting of the experimentally-measured transmission or reflection spectra to theoretical models, allowing the film physical thickness data recovery, if detailed refractive index information is available. Fitting of the absorption coefficient spectra in absorbing material layers is possible, using an automated algorithm reliant on the data for the measured transmission spectrum, refractive index spectrum, and physical thickness. A number of experimentally- and computationally-derived optical constant-related datasets of different magneto-optic garnet compositions possessing giant Faraday rotation have also been reported.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-7080/7/3/49/s1, Zip archive containing program installation files. A number of MO garnet material-related datasets are available from the program installation directory after installation.

**Author Contributions:** M.V. wrote the solution code for the program using Microsoft Visual Studio 2003. NET Professional and multiple code changes and feature additions have been applied between 2005–2019. During this time frame, the development of the algorithm and code features of the program have been greatly influenced by the academic discussions within our group (K.A. and M.N.-E.-A.), and on-going work on multiple materials-related projects. Both K.A. and M.N.-E.A. contributed to the preparation of the final manuscript and its pre-submission review. All co-authors have performed multiple laboratory and computational experiments with a range of magneto-optic material types using the reported software tools.

**Funding:** This research received no external funding.

**Acknowledgments:** This work was supported by the Edith Cowan University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A : Program Summary and Technical Details**

*Program title:* Optimization of 1D Magnetic Photonic Crystals (alternatively, MPC Optimization)

*Program installation files doi:* available from the Supplementary section of manuscript.

*Licensing provisions:* Creative Commons Attribution-NonCommercial-3.0 Unported (CC BY-NC-3.0)

*Programming language:* Visual C#, compiled using Microsoft Visual Studio .NET 2003

*Nature of computational problem:* Calculation of the optical transmission, reflection, and Faraday rotation spectra in multilayer thin films containing gyrotropic constituents (magnetized material layers possessing magneto-optic properties); optimization of magnetic photonic crystal (MPC) designs aimed at achieving maximized transmission or reflection coincident with maximized Faraday rotation, according to sets of defined criteria; fitting of the experimentally measured transmission or reflection

spectra to theory models; and fitting of the absorption coefficient spectra of single-layer thin film materials using the data for optical transmission, film thickness, and refractive index spectra.

*Implementation:* The program exhaustively calculates multiple possible multilayer structure designs, based on the design structure type(s) defined prior to running optimization. complex-valued 4 × 4 transfer matrix method (accounting for all dielectric tensor components, including the non-diagonal terms responsible for gyrotropic effects) is implemented to compute the complex field amplitudes and optical intensities in either the transmitted or reflected left-hand and right-hand circularly-polarized eigenwaves propagated through the thin-film substrate structure.

*Restrictions:* The program is designed for use in conjunction with reliable optical constant datasets for up to three component dielectric materials, one of which can be modeled as magnetic dielectric possessing Faraday rotation; metallic layers are not implemented. The embedded ActiveX controls, which enable graphical data output, accept only up to 1000 data points per graphing control, whether plotting a single curve or several.

## **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
