Marcin Lefik

Krzysztof Komeza
email: krzysztof.komeza@p.lodz.pl

Ewa Napieralska-Juszczak

Poitr Daniel Roger

Piotr Napieralski

Andrej Napieralskik

Comparison of the reluctance laminated and solid rotor synchronous machine operating at high temperatures

Keywords: Electrical machine, Thermal analysis, Finite element method

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

High-temperature (HT) electrical machines are intended for specific applications. They are designed to work at higher temperatures than those typically used. Of course, in such an environment, the temperature within the machine also increases. For example, the synchronous reluctance motor (SynRM) designed for aircraft applications works at ambient temperature of around 200°C, whereas the temperature within the motor's frame reaches 450-500°C at the hottest spot. With electrical machines capable of running in such an environment, engine manufacturers have been able to reduce the reliance on mechanically driven accessories and, according to the More Electrical Engine concept, provide a better product in terms of efficiency and reliability. Continuous operation in such conditions requires the use of inorganic Electrical Insulation Systems (EIS), as traditional organic materials are not capable of withstanding such a high temperature, and an architecture adjusted to the inorganic insulation technology (Iosif et al., 2016a(Iosif et al., , 2016b;;Cozonac et al., 2014). This topology allows two main options: permanent magnet synchronous machine or synchronous reluctance machine, in which the torque is generated through the phenomenon of magnetic reluctance (Pyrhoonen et al., 2008;Belhadi et al., 2015;Lipo, 1991). The HT engine built of rigid inorganic coils requires a special stator construction with rectangular teeth shape and open slots. Therefore, design calculations of these electromechanical transducers must be preceded by a specific thermal and electromagnetic analysis, aimed to estimate the thermal effects of power loss. This paper presents a 3D FEM-based electromagnetic and thermal analysis of the HT synchronous reluctance machine (SynRM) working with the rotational speed 4800 rpm. The analysis relies on 3D computational fluid dynamics (CFD). Due to the specific stator design linked to the coil winding technology and a large number of poles, the choice of the rotor construction was fairly limited. The present study is based on a salient pole type of rotor. To increase the torque, the poles were notched lengthwise and a round opening at the base of the rotor tooth was added to serve as a flux barrier. Figure 1 shows the main parts of the reluctance motor. Both electromagnetic and CFD solutions take into account changes in parameters and characteristics of materials associated with high operating temperature (Cozonac et al., 2014;Finger and Rubertus, 2000;Ababsa et al., 2018;Komeza et al., 2016Komeza et al., , 2017)).

The rotor can be made as a solid or as a laminated. Solid construction is mechanically more robust allowing higher rotation speed; moreover, induced currents can generate torque during starting. Unfortunately, stator construction, even with used stator slot wedges, generates higher flux harmonics which produced high losses in rotor at synchronous speed. Induced in solid rotor poles currents weakened the effect of reluctance changes. Figure 2 shows comparison of torque versus load angle for both: laminated and solid design. Evidently, a motor with a massive rotor has a much smaller synchronous torque as well as a smaller range of stable operation. This is due to the effect of shielding the rotor through currents induced in its structure.

Electromagnetic field simulation

To obtain heat sources for the thermal model, calculations of the electromagnetic field were made using the Opera 3D program. As commonly known the 3D time-stepping FEM is very time-consuming (Ionel and Popescu, 2010;Bianchi et al., 2013;Ponomarev et al., 2014); thus, the use of 2D analysis is preferred (El-Refaie, 2010;Komeza et al., 2016;El-Refaie and Jahns, 2005;Tang et al., 2015). The time-stepping analysis was performed over approximately 20 cycles until the steady state was reached. First, 2D and after 3D simulation were realized. For the purpose of the 3D simulation, a quarter of the machine has been considered to benefit from the magnetic symmetry. The time-stepping analysis was performed assuming constant speed of rotation of the rotor. Magnetic vector potential formulation with linear edge elements was used in that analysis. The analysis relies upon a remeshing technique. During the transient analysis, before the solution at the next time-step, the position of each moving part is updated. The 3D mesh contains 4.264778 million elements in total (including 700,000 in gap region). In the conductive areas, the elements size was chosen to be less than one-third of the penetration depth. Moreover, the distribution of the density of the mesh is greater than to correctly reproduce the effect of rotation of the rotor relative to the stator. To reduce the computation times, the simulations were performed on a virtual machine with eight cores. Figure 3 shows finite element mesh for 3D model of the quarter of the machine. Figure 4(a) and 4(b) shows the comparison of torque for laminated and solid rotor (made from steel C10) calculated by 2D and 3D models. While for laminated rotor, 2D and 3D results were almost the same, for solid rotor, they differ significantly.

The reason for such a large difference between a 3D and 2D solution for a motor with a massive rotor is the complex distribution of currents [show in Figure 5 Even for synchronous speed, higher harmonics cause significant losses in the conductive material of the rotor (despite the elevated temperature the conductivity is 2.24 MS/m). Figure 7 shows rotor losses versus time for one period. The total losses in the rotor are practically independent of the load angle and amount to 1030 W for the maximum torque.

3D coupled thermal fluid model

To analyse the thermal phenomenon, the 3D coupled thermal-fluid (CFD) model is used (Unai et al., 2014;Arbab et al., 2015). The CFD model allows to calculate the distribution of the temperature field and average temperature of the windings and the other motor's parts. In addition, this model can be used to localize points with the highest value of the temperature (hot points). The CFD model is based on the Reynolds-averaged Navier-Stockes equations (RANS equations) for the steady flow. RANS equations can be written as follows (Grybos, 1998): where V is an averaged values of fluid velocity, P is averaged values of fluid pressure, T is averaged values of fluid temperature, Q is a heat source density, t is a Reynolds-stress tensor, D is an intensity of dissipation and k; c; r; m are thermal conductivity, specific heat, mass density and dynamic viscosity, respectively. Equations ( 1)-( 3) can be solved only with turbulence models. The most commonly used types of turbulence models are k-« models. The k-« models include two additional differential equations which are written for turbulent kinetic energy k and dissipation «. In analysed case k-« Realizable turbulent model is used. In rotation machines, it is necessary to take into account the rotor rotation and thermal radiation. Both phenomena affect the heat flow in analysed machines. Rotation of the machines rotor can be modelled using multiple rotating reference frames (MRF model), where the equations for rotor are written in moving reference frame and for the other machine's parts in stationary reference frame. Radiation inside the machine's housing can be included by radiation model, for example discrete ordinates, and for external walls by combined boundary condition for energy equation with the surface emissivity coefficient « r . Combined boundary condition is given by following formula [18]:

where h c is a convective heat transfer coefficient, T 1 is temperature of the ambient air, T w is surface temperature of the wall, T ext is the ambient temperature and s is the Stefan-Boltzmann constant. The convective heat transfer coefficient h c can be determined from empirical formula which is based on dimensionless Rayleigh and Nusselt numbers. Heat sources in the analyzed machine are the power losses. The power losses in the core and windings of HT° reluctance synchronous machine were calculated. In equation ( 3), the heat source is presented as heat source density Q. It can be calculated from following formula:

Summary

In this study, we analyzed the influence of the design of the rotor (solid or laminated) on efficiency and temperature of the motor active parts. The presented results show clearly that laminated construction is much better from a point of view of efficiency and temperature. However, solid construction can be interesting for highspeed machines due to their mechanical robustness.

The influence of ambient temperature on the winding's coils temperature was also analysed. The results have shown that the optimal temperature for the machine with laminated rotor is 200°C. In this case, the average temperature of the coils is 393.8°C, while the estimated operating temperature of the coils is 450-500°C. The coils are overheated for a higher ambient temperature.

Figure 1 .

 1 Figure 1. Reluctance motor.

  (a)], and the resulting changes in the magnetic flux density distribution [Figure 5(b)]. Only in the stator area as shown in Figure 6(a), the field distribution is analogous as for the motor with laminated rotor [Figure 6(b)].

Figure 2 .

 2 Figure 2. The comparison of torque for laminated and solid rotor design.

Figure 3 .

 3 Figure 3. Finite element mesh for 3D model of the quarter of the machine contains 4.264778 million elements.

Figure 4 .

 4 Figure 4. Comparison of torque 1,500 for laminated rotor design calculated by 1,000 2D and 3D models (a), Comparison of torque or solid rotor design calculated by 2D and 3D models (b).

Figure 5 .

 5 Figure 5. The distribution of current induced in solid rotor [A/mm 2 ] (a); the distribution of magnetic flux density [T] (b).

Figure 6 .Figure 7 .

 67 Figure 6. The distribution of magnetic flux density in stator: for the motor with solid rotor (a); for the motor with laminated rotor (b)

Figure 9 .

 9 Figure 9. Temperature distribution in machine with laminated (a), (b), (c) and solid (d), (e), (f) rotor for the ambient temperature 200°C.

Figure 10 .

 10 Figure 10. Temperature distribution in machine with laminated (a), (b), (c) and solid (d), (e), (f) rotor for the ambient temperature 500°C.

where DP is power losses taken from electromagnetic calculations, whereas V is a volume of motor's part where these power losses were generated.

Results of the thermal analysis

Thermal analysis allowed to compare thermal phenomena occurring in both variants of the designed machine. Thermal field distributions in machines and average temperature values of individual machines parts for two values of ambient temperature (200°C and 500°C) were calculated. In these calculations, the effect of temperature on the change of winding resistance, solid rotor conductivity and the associated change in power losses was taken into account. The results of thermal analysis in Table I, Figures 9 and 10 were shown. The obtained results show that machine with solid rotor has higher values of temperature than machine with laminated rotor. Additional losses in solid rotor cause that other elements of this machine they overheat, in particular the winding coils. The average value of the winding temperature for a laminated rotor at an ambient temperature of 200°C is 393.8°C, while for a solid rotor it is 558.6°C. As can be seen in the case of the massive rotor, the estimated operating temperature (450-500°C) of the coils is exceeded. A similar situation exists for the machine operating at ambient temperature 500°C. The average temperature is 612.9°C for machine with laminated rotor and 697.2°C for solid rotor.


Shaohui Xin

Qiang Wang
email: qiangwang@wipm.ac.cn

Jun Xu

Yueying Chu

Pengfei Wang

Ningdong Feng

Guodong Qi

Julien Trébosc

Olivier Lafon

Weibin Fan

Feng Deng
email: dengf@wipm.ac.cn

Acidic Nature of "NMR-invisible" Tri-coordinated Framework Aluminum Species in Zeolite †

Unambiguous characterization of different acid sites in zeolites is of greatimportance for understanding their catalytic performance and rational design of highly efficient zeolite catalysts. In addition to various well-characterized extra-framework Al species, tri-coordinated framework aluminum species can also serve as Lewis acid site in zeolites, which is however "NMR-invisible" due to its extremely distorted local environment. Here we provide a feasible and reliable approach to elucidate the acidic nature of the tri-coordinated framework Al in dehydrated H-ZSM-5 zeolites via sensitivity-enhanced two-dimensional multiple nuclear correlation NMR experiments coupled with trimethylphosphine oxide (TMPO) probe molecules. Two types of tri-coordinated framework Al sites have been unambiguously identified, which amount to 11.6% of the total Brønsted and Lewis acid sites. Furthermore, it is found that the synergistic effect arising from the close spatial proximity between tri-coordinated framework Al site and Brønsted acid site leads to generation of superacidity (with acid strength stronger than 100% H2SO4) in the zeolite.

Introduction

Heterogeneous catalysts with acid-base prosperities, such as zeolites have been widely used in chemical and petrochemical industry for catalytic cracking, 1 isomerization, 2 alkylation, 3 and disproportionation. 4 In aluminosilicates-type zeolites, Brønsted acid site (BAS) and Lewis acid site (LAS) are responsible for the active sites in numerous catalytic reactions. Therefore, understanding the structures and properties of the acidic sites in zeolites is essential to explore the reaction mechanism and to optimize their catalytic performance. [5][6][7] It is well-known that the BAS is associated with the tetrahedral framework aluminum (FAL) in the form of bridging hydroxyl (Si-OH-Al) in zeolites. And the LAS is generally recognized to extra-framework or framework aluminum species, formed upon calcination or steaming of zeolites. 8,9 The properties of BAS associated with tetrahedral framework Al and LAS associated with various extra-framework Al (EFAL) species have been extensively studied, from their structure characteristics [10][11][12][13][14] to their functions in catalytic reactions. [15][16][17][18] Especially, the LAS in zeolite has proved to play an important role in catalytic transformations of hydrocarbons such as hydrogen transfer process, 7,19,20 and also generate the Brønsted/Lewis acid synergy owing to the vicinity of corresponding Al species in the local architecture of dealuminated zeolites. 11,[21][22][23] Although the tri-coordinated FAL is commonly considered as LAS in zeolites, only a few characterizations [24][25][26] have been attempted to confirm the exsitence of "moisture sensitive" tricoordinated Al species in zeolite framework. Up to now, no experimental approach has been reported to probe the intrinsic acidic property of such important framework Al species in zeolite.

Solid-state NMR has emerged as an important technique for the characterization of heterogenous catalysts at the atomic scale. 22,[27][28][29][30][31] The coordination state of both FAL and EFAL species in zeolites can be determined on the basis of onedimensional (1D) 27 Al magic-angle spinning (MAS) and 2D 27 Al multiple-quantum magic-angle spinning (MQMAS) NMR spectroscopy. 2D 1 H-1 H and 27 Al- 27 Al double-quantum single-quantum (DQ-SQ) homonuclear correlation MAS NMR methods can provide in-depth structural information on the spatial proximity of different acid sites in dealuminated zeolites, especially under high magnetic field ( > 18.8 T). 11,22,32 Recent report has shown that the perturbed aluminum species was attributed to terminal Al-OH in hydrated zeolites, which was considered to be derived from hydrous tri-coordination of Al species. 33 However, most of these measurements were performed under hydration condition. In the case of non-hydrated zeolites, it is difficult to distinguish FAL species due to the resolution/sensitivity issues from the quadrupolar nature of the 27 Al isotope. The stretch of acidic proton in BAS induces a larger distortion of the corresponding tetrahedral FAL in H-form zeolites, leading to its quadrupolar coupling constants (CQ) up to 14 ~ 18 MHz. [34][35][36] Generally, mild calcination and dehydration treatments of zeolites would lead to partial breaking (or hydrolysis) of framework ≡Al-O-bond and generation of tri-coordinated FAL and framework silanol group as illustrated in Scheme 1. [37][38][39] Although the neighboring framework Si-OH would interact with the tri-coordinated FAL to restore tetrahedral FAL after cooling down the sample to room temperature, partial Si-OH derived from the breaking of aluminum oxygen bond that bears the proton may be further dehydroxylized with a neighboring defect site (Si-OH group), thus causing an irreversible formation of three-fold coordinated aluminum. 26 However, reducing the number of neighboring oxygen bound to FAL yields a dramatic distortion of its local environment, and the extreme asymmetry of its surrounding electric field results in a considerable line broadening of its 27 Al NMR resonance (usually with CQ > 30 MHz). 33,40 Therefore, the direct observation of tri-coordinated FAL species in dehydrated sample is generally impossible by conventional 27 Al MAS NMR due to its huge quadrupolar broadening and relative low concentration, and the tri-coordinated FAL was suggested to be the "NMRinvisible" species in zeolites. 26,33 After rehydration, the tri-coordinated FAL could be easily transformed into symmetric tetra-coordinated FAL, 33,38 which is usually indistinguishable from the FAL of BAS. Further hydrolysis of framework Al-O bond leads to removal of Al from zeolite framework, forming EFAL species such as Al(OH)3, Al(OH)2 + , AlOH 2+ , Al 3+ and AlO + etc., which have been well characterized by solid-state NMR spectroscopy. 11,22,32,36,41,42 The adsorption of basic probe molecules such as 2-13 C-acetone and trimethylphosphine oxide (TMPO) is a feasible approach to investigate the acidity of zeolites. 22,[43][44][45] The great majority of studies utilized the chemical shift of certain elements (e.g. 13 C or 31 P) in probe molecules to determine the acidic features of zeolites, which were usually combined with density functional theory (DFT) calculations to reveal the relationship between chemical shifts and the corresponding acidic properties (e.g. type, distribution, and strength) on zeolites. 22,45 However, hitherto it is a great challenge to provide precise and localized information on the structure of the interface between the "guest" probe molecules and the specific acid sites of "host" zeolite via analytical or spectroscopic techniques.

In this contribution, we provide a unique insight into the acidic nature of tri-coordinated FAL species on H-ZSM-5 zeolite by solid-state NMR spectroscopy, and present a clear adsorption picture of TMPO probe molecules on different framework Al species in dehydrated H-ZSM-5 zeolite. In particular, the sensitivity-enhanced 2D 31 P- 27 Al heteronuclear correlation (HETCOR) MAS NMR techniuqe 46 is employed to successfully discriminate the interactions between distinct adsorbed TMPO molecules and Al species in the zeolite framework. Thus, it is straightforward to ascertain the presence of LAS originating from tri-coordinated FAL. Two tri-coordinated FAL species are unambiguously identified by the advanced solid-state NMR technique. Further quantitative NMR analysis reveals that a considerable amount (ca. 11.6%) of this type LAS exists in the framework of dehydrated H-ZSM-5 zeolite even with moderate thermal treatments. Moreover, we demonstrate that the presence of the tri-coordinated FAL is responsible for the formation of superacidity (characterized by a down-field 31 P chemical shift of adsorbed TMPO up to 85-88 ppm) due to the Brønsted/Lewis acid synergy in the zeolite framework which is further confirmed by 2D 27 Al- 27 Al DQ-SQ homonuclear correlation spectroscopy and DFT calculations.

Results and discussion

Structure characterization of dehydrated zeolites

The 27 Al MAS NMR spectra of hydrated parent NH4-ZSM-5 acquired at two different magnetic fields (18.8 and 11.4 T) in Fig. 1a and 1d both displayed a typical narrow line around 54 ppm, respectively, corresponding to relatively symmetric tetrahedral Al sites in zeolite framework. Usually much less quadrupolar broadening of the 27 Al NMR lineshape can be obtained at high field (B0 = 18.8 T) and fast MAS speed (νR = 40 kHz) than that at medium field (B0 = 11.7 T) and moderate MAS speed (νR = 10 kHz). While after calcination and dehydration, the 27 Al spectra of dehydrated H-ZSM-5 zeolite exhibited a very broad resonance (Fig. 1b and Fig. 1e) with a full width of ca. 20 kHz even at 18.8 T, indicative of a dominant contribution from the large second-order quadrupolar broadening. Previous study has demonstrated that this broad signal was mainly attributed to the more-distorted tetrahedral FAL species with a CQ ≈ 16 MHz from Brønsted acidic site on dehydrated zeolites. 35,36 In fact, dealumination is inevitable in most sample preparation procedures, while the mild vacuum or flow temperature-programmed methods used here couldmake framework dealumination rarely occur, which would largely reduce the interference of EFAL species in the study of tri-coordinated FAL. The XRD patterns (Fig. S1) of ZSM-5 zeolites indicated no apparent change of crystallinity after these treatments. To estimate the amount of extra-framework Al in our calcined sample, the parent and calcined forms of these zeolites were treated by 1 M aqueous NaNO3 at 353 K for three times. Since the Na + cations can balance the charge of framework AlO4 -tetrahedra, the Na/Al ratio is a reliable parameter to assess the variation of fraction of tetrahedral framework Al. 47 The Na/Al ratio measured by ICP was 1.01 and 0.98 respectively for parent NH4-ZSM-5 and dehydrated H-ZSM-5 zeolites treated with Na + exchange, suggesting that most of the Al species were still located in the zeolite framework. This is also supported by a slight loss (1%) of framework ≡Si-O-Al moieties in the dehydrated zeolite via 1D 29 Si MAS and 2D 29 Si{ 1 H} HETCOR NMR analysis (Fig. S2,S3). The quantitative 29 Si MAS spectra in Fig. S3 indicated an increase of Si-OH groups in dehydrated H-ZSM-5 zeolite compared with those in parent NH4-form one. In addition, further 1 H MAS NMR analysis (Fig. S4) showed that the amount of Si-OH groups decreased when raising the dehyrated temperature from 573 K to 673 K, suggesting that the dehydroxylation of Si-OH groups may occur as depicted in Scheme 1.

Framework aluminum studied by 2D 27 Al 3QMAS

In our measurements, TMPO molecules were adsorbed on H-ZSM-5 zeolite with their content much lower than that of framework acid sites, by which the interference between TMPO molecules themselves can be largely avoided. Interestingly, besides the narrow 27 Al signal at ca. 54 ppm, a new broad shoulder feature (at ca. 45 ppm) was visible in the 1D 27 Al MAS spectra of dehydrated zeolite loaded with TMPO (P/Al = 0.42, Fig. 1c,f). Since the tri-coordinated FAL with Lewis acidity (if present) can directly interact with the oxygen atom of TMPO, which would significantly decrease the asymmetry of the local environment of the Lewis acidic FAL, yielding its "NMR-detectable" 27 Al resonance on non-hydrated zeolites. Then, 2D 27 Al 3QMAS experiment at 18.8 T was employed to discriminate specific aluminum sites in the dehydrated and TMPO-loaded sample (Fig. 2a). Three signals, denoted as Ala, Alb and Alc, were well-resolved with the isotropic chemical shifts being 54.4, 51.6 and 58.8 ppm, respectively. For the dehydrated zeolite with a lower TMPO concentration (P/Al= 0.18), an obvious decline of the relative proportion of Alc and Alb signals (Fig. 2b,S10) was evident. The NMR parameters deduced from the 2D spectrum were listed in Table 1. The Ala signal with PQ = 2.3 MHz could be attributed to the regular tetrahedral FAL species from Brønsted acidic sites analogous to that of NH4-form or hydrated zeolites. The local structure of tetrahedral FAL in dehydrated H-ZSM-5 became more symmetrical stemming from a shrinkage of the Al-O(H) bond length in BAS after TMPO adsorptions. According to their isotropic chemical shifts, Alb and Alc should be both associated to tetrahedral FAL species as well, while exhibiting a larger inhomogeneous second-order quadrupolar broadening (PQ > 7 MHz). However, their assignments were ambiguous, since they could be from either the FAL in another type of BAS, or the tri-coordinated FAL of LAS directly bound to a TMPO molecule, yielding a new distorted tetrahedral FAL state. Especially, the Alb and Alc resonances contributed to considerable fraction in the quantitative 27 Al NMR spectra of dehydrated zeolites with different TMPO loadings (Fig. 1c). Therefore, further experiments are required to identify the host-guest correlations/interactions between the adsorbed TMPO probe molecules and the framework of zeolites.

Table 1. NMR parameters obtained from 27 Al MQMAS spectra. a δF2 and δF1 represent the center of gravity of each particular species measured along direct dimension F2 and isotropic dimension F1 on sheared spectra, respectively. 1b Isotropic chemical shift δCS, quadrupole interaction product PQ of each aluminum species were deduced from δF2 and δF1. c The quadrupolar coupling constant (CQ) were extracted from the slice of 3QMAS spectra by fitting the corresponding second-order quadrupolar lineshape (η = 0.5) with the DMFIT program.

Site

Acid properties analyzed by 2D 31 P { 1 H} HETCOR

For the acidity characterization by TMPO probe molecules, the dichloromethane (CH2Cl2) is usually used as the solvent in order to well disperse the TMPO molecules into zeolite. 45,48 Although the evacuation at 323 K seems to be insufficient to dispose of all solvent, in the 2D 13 C { 1 H} and 31 P{ 1 H} HETCOR spectra of TMPO/H-ZSM-5 zeolite (Fig. S5) acquired with relative short CP contact time (4 ms or 1ms) there is no correlation between CH2Cl2 and TMPO molecules, indicating that the residual solvent and adsorbed TMPO molecules are well separated. Furthermore, the 27 Al MAS NMR spectrum of H-ZSM-5 zeolite with solely adsorbed CH2Cl2 (Fig. S6) still exhibits a very broad resonance (△δ ≈ 90 ppm), reflecting that the residual CH2Cl2 hardly interacts with acidic sites on zeolite, thus almost has no influence on the acidity characterization by TMPO probe molecules. Note that, after removing the residual solvent on TMPO/H-ZSM-5 zeolite at a higher temperature (363 K, Fig. S7), the 27 Al MAS NMR spectrum of TMPO/H-ZSM-5 was the same as that treated at 323 K (Fig. S8), which clearly indicated that the influence of solvent (CH2Cl2) residues on the framework Al species can be negligible. 1D 31 P MAS NMR spectra of H-ZSM-5 zeolites with different TMPO loadings (P/Al = 0.18 and 0.42) were given in Fig. 3a and3c, respectively. Up to six characteristic peaks at 88, 85, 76, 69, 65 and 51 ppm were identified from the 31 P MAS NMR spectra by using Gaussian deconvolutions (also see Table S1). Indeed, quantitative analyses of 27 Al MAS NMR are difficult due to either different quadrupolar broadenings of distinct FAL species or only partial acid sites interacting with TMPO probe molecules, however, here the 31 P (I =1/2) MAS NMR spectra can provide quantitative information on the framework BAS and LAS that react with TMPO molecules. The 2D 31 P { 1 H} CP-HETCOR spectrum of the TMPO/H-ZSM-5 zeolite (P/Al = 0.18) acquired with a CP contact time of 8 ms in Fig. 3e exhibited two broad distributions of correlation peaks in "sloped" oval (dotted line) centered at (76, 12.2) ppm and (88, 6.4) ppm, respectively, which could be assigned to distinct protonated TMPOH + ionic pair complexes. These correlations were also observable when a short CP contact time (1 ms) was used as (see Fig. S5). The interaction of adsorbed TMPO molecules with BAS in H-ZSM-5 zeolite will result in a downfield shift of 1 H resonance of the acidic proton, and the electron cloud density surrounding the 31 P nucleus in TMPO molecule decreases with increasing acid strength of BAS, leading to the 31 P resonance moved toward higher chemical shift (downfield). 22,44,45 Furthermore, it was also found that the 31 P chemical shift of adsorbed TMPO molecules linearly increases with the increase of Brønsted acid strength, and a 31 P chemical shift of ca. 86 ppm was demonstrated for the threshold of superacidity (with acid strength stronger than 100% H2SO4). 45 Therefore, the 31 P peak at δ31P = 88 ppm reflected the strongest acid site with strength over the threshold of superacidity in the zeolite. In addition, the correlation peaks between all 31 P signals and methyl groups (δ1H = 1.8 ppm) of adsorbed TMPO molecules can be observable all the time. Since the 1 H-1 H spin diffusion will promote the detection of long range interactions in the case of a long CP contact time (8 ms), weak correlations between the protons of the residual CH2Cl2 solvent (δ1H ≈ 5.4 ppm) and 31 P of TMPO molecules were also detectable in the 2D spectrum. However, as illustrated in the 2D 31 P { 1 H} HETCOR spectrum of Fig. 3e, even though 1 H-1 H spin diffusion also contributed the 1 H-31 P correlations when the CP contact time was set to 8 ms, no correlation can be observed for the 31 P signals around 65 ppm with the 1 H signals of BAS, which suggested that the 31 P signals should not originate from the protonated TMPOH + complexes. It is generally recognized that H2O could competitively react with Lewis acid sites, such as tri-coordinated Al species, to form weak BAS. 48 After exposing our samples to humidity for 2 h, the corresponding 31 P MAS NMR spectra in Fig. 3b, 3d exhibited an apparent decline in relative intensity of the signals at 65 and 69 ppm (also see Table S1), and the signal at 88 ppm (or 85 ppm for the sample with P/Al = 0.42) disappeared completely. In addition, the peak at 76 ppm was slightly shifted to upfield with a broader distribution, accompanying a notable increase of the peak at 51 ppm from the TMPO adsorbed on weak acid sites. Therefore, we can reasonably consider that the 31 P signals at 65 ppm and 88 ppm (or 85 ppm) are both associated to Lewis acid sites directly or indirectly. Although the peaks at ca. 88 ppm have proved to be from TMPO molecules adsorbed on Brønsted acid sites, it seems that BAS together with neighboring LAS should contribute to the formation of the superacidity (e.g. δ31P = 88 ppm), indicative of the presence of synergy effect. Accordingly, after humidity exposure of the dehydrated zeolite, the LAS was preferentially hydrated, thus the influence on its neighboring BAS (the superacid feature) was vanished, which in turn caused the BAS back to its original acidity. The two 31 P signals at 65 and 69 ppm were probably from TMPO molecule adsorbed on Lewis acidic FAL species. Concerning the abundance of Alb and Alc in Fig. 2 (also in Fig. S10), we supposed that the aforementioned 31 P and 27 Al species should belong to interacted phosphorusaluminum pairs in the Lewis acid complexes, (CH3)3PO••Al≡. In addition, we also acquired the 2D 31 P{ 1 H} HETCOR spectrum (Fig. S9) on the TMPO/H-ZSM-5 zeolite with the complete removal of solvent at 363 K, which exhibited the similar 1 H-31 P correlations between TMPO molecules and framework acid sites in zeolite compared with that treated at 323 K, which further excluded the influence of residual solvent. Tri-coordinated framework Al characterized by 31 P { 27 Al} HMQC and 31 P { 27 Al} S-RESPDOR To verify the assumption, the Lewis complex (CH3)3PO••Al≡ was identified directly by using our recently developed 31 P- 27 Al population-transfer dipolar-mediated HMQC (PT-D-HMQC) experiment which can establish the direct correlation/ connectivity between TPMO and various Al species. As shown in Fig. S5, the 1D 31 P { 27 Al} PT-D-HMQC MAS NMR spectra via 31 P→ 27 Al→ 31 P filtering on TMPO-loaded H-ZSM-5 zeolite (P/Al = 0.42) revealed that the two main 31 P resonances around 76 and 65 ppm were correlated with 27 Al species in zeolite. When using a short recoupling time (τmix = 1.33 ms), the peaks at ca. 65 ppm became more intense compared to that at ca. 76 ppm (Fig. S11). This indicated the stronger 31 P- 27 Al dipolar interaction of corresponding spin pairs, further suggesting the presence of the Lewis complex [(CH3)3PO••Al≡].

Notably, the 2D 31 P { 27 Al} PT-D-HMQC spectrum acquired at 11.7 T of TMPO-adsorbed H-ZSM-5 zeolite in Fig. 4b clearly illustrated two correlation groups. The correlation peak at around (76, 54) ppm was associated with TMPO adsorbed on the bridging hydroxyl (Si-OH-Al) proton (i.e. Brønsted acid sites). Importantly, other correlations indicated that the 31 P resonance peaks at 65 and 69 ppm are both correlated with broad 27 Al resonances from distorted four-coordinated aluminum, reflecting the presence of two types of Lewis acid sites. We also conducted 2D 27 Al { 31 P} D-HMQC experiments at 18.8 T. Due to the short transverse relaxation times (T2) of 27 Al resonances, we could only observe 27 Al-31 P correlations from the Lewis complex [(CH3)3PO••Al≡] in Fig. 4c, in which 31 P signals at 65 and 69 ppm were correlated with Alb and Alc, respectively (Fig. 4d,4e). Meanwhile, a reduction of the TMPO concentration led to a decrease in the proportion of 31 P signal at 69 ppm relative to that at 65 ppm (Table S1), in consistent with the relative proportion variation of Alc and Alb (Fig. 2). Therefore, we can confirm that Alc and Alb are both from tri-coordinated FAL species bound to the O atom of TMPO molecules, which were well discriminated from their different Lewis acid strength.

Due to its low content and Brønsted acid feature, we didn't observe the 31 P- 27 Al correlations from the 31 P signal at 88 ppm (or 85 ppm) in Fig. 4. Fortunately, Symmetry-based Resonance-Echo Saturation-Pulse DOuble-Resonance (S-RESPDOR) method 49 can be employed here to estimate 31 P- 27 Al heteronuclear dipolar interactions. The results are usually analyzed by plotting the signal fraction ΔS/S0 = (S0 -S')/S0 (S' and S0 represent the signal intensity with and without dipolar dephasing, respectively) as function of the recoupling time τ. The difference 31 P spectrum (ΔS = S0 -S', shown in Fig. S12) obtained by 1D 31 P { 27 Al} S-RESPDOR experiments still suggested the spatial interactions/proximities between all 31 P species of adsorbed TMPO molecules (including that giving rise to the resonance at 88 ppm) and 27 Al species in the framework of H-ZSM-5. Furthermore, the lower concentration of TMPO molecules (P/Al = 0.18) allowed us to estimate distances between 31 P atoms and 27 Al atoms, since we could approximately consider the Lewis complex [(CH3)3PO••Al≡] with 31 P signals at 65 ppm or Brønsted acid complex [(CH3)3PO••H-OAl≡] with 31 P signals at 76 ppm as "isolated" phosphor-aluminum spin pairs in this case. Therefore, the P-Al distance of the former was determined to be 2.9 ± 0.2 Å by fitting the dephasing curve (Fig. 5a), which was rationalized in terms of the structure model of TMPO molecule adsorbed on the framework tri-coordinated aluminum species (Fig. 5c). While the later complex had a reasonable P-Al distance of 4.1± 0.4 Å (Fig. 5b), in which 31 P atom was about four bonds away from the Al site (Fig. 5d). We also considered the dipolar coupling networks of multiple-spin systems (shown in Fig. S13). The representative three spin system are depicted in Fig. S13b-d. In the normalized dipolar dephasing curve with respect to first extremum (or maximum) in each spin system, it is found that the dephasing curves of the observed 31 P signals in multi-spin systems have similar upward trend to that of the isolated 31 P- 27 Al spin pair. This result suggests that the upward trend of the dipolar dephasing curve is dominated by the short-range (strongest) 31 P- 27 Al dipolar interaction. The long-range couplings have little impact on the measurement of the short distance, which mainly affect the oscillation of the dephasing curves after reaching its first extremum. Moreover, the 31 P- 27 Al distances extracted by the dipolar dephasing curves were in good agreement with DFT calculations, where the theoretically predicted 31 P- 27 Al distances of the Lewis acid complex [(CH3)3PO••Al≡] and Brønsted acid complex [(CH3)3PO••H-OAl≡] were 3.1 Å and 4.1 Å (Fig. 5c and5d), respectively.

Although it is difficult to calculate the accurate amount of the tri-coordinated FAL in our sample, an estimation of the lower limit of their concentration via the quantitative analysis of 31 P MAS NMR spectrum (Fig. 3) is available with the following equation:

Ω(FAL Ⅲ ) = χ(P/Al) * Ω( 31 PLAS)
Where Ω(FAL Ⅲ ) is the fraction of tri-coordinated FAL species interacted with TMPO in the total Al content, Ω( 31 PLAS) = 27.0% is the proportion of 31 P signals at 65 and 69 ppm in the quantitative 31 P MAS NMR spectrum (see Table S1), and χ(P/Al) = 0.42 is the atomic ratio of P to Al of the sample obtained by ICP. Therefore, the calculated Ω(FAL Ⅲ ) = 11.6% suggested that a considerable amount of LAS formed and located in the framework of dehydrated H-ZSM-5 zeolite, in agreement with the relative content of Lewis acid sites (13.3%) determined by our FT-IR analysis of pyridine adsorption measurements (Table 2 and Fig. S14).

Brønsted/Lewis synergy effect studied by 2D 27 Al DQ-SQ NMR and DFT calculations

We have successfully obtained the high-resolution 27 Al- 27 Al DQ-SQ homonuclear correlation spectra on hydrated zeolites (HY, MOR, ZSM-5, etc.) by using our developed BR2 1 2 recoupling scheme at high magnetic field, 32,41 whereas only Al species from framework BAS and extra-framework LAS were considered. In the case of non-hydrated zeolite catalysts, the broadening and overlap of the 27 Al quadrupolar patterns usually result in considerable barriers to implementation of the 2D 27 Al DQ-SQ experiment. Fortunately, the improved resolution of 27 Al MAS NMR spectra of dehydrated H-ZSM-5 zeolite after TMPO adsorptions offers an opportunity to investigate the intrinsic structure of different framework Al species under dehydrated state, whose acidic natures have been clearly characterized. Fig. 6a displayed the 27 Al- 27 Al DQ-SQ MAS NMR spectrum of TMPO-adsorbed H-ZSM-5 with P/Al = 0.42 at 18.8 T. Two auto-correlation peaks were observable, in which the intense blue part indicated that tetra-coordinated framework Al species (Ala) was in close proximity one another (Ala), and the weak orange part was ascribed to the spatial proximity of Alb-Alb or Alc-Alc pairs in a low content. Note that, it is still hard to discriminate the Alb and Alc in the 27 Al DQ-SQ MAS NMR spectrum. In addition, the appearance of an intense cross-peak pair between Ala and Alb (Red part in Fig. 6a) provided a direct evidence on the spatial proximity between the framework Brønsted and Lewis acid sites, where Ala and Alb had been confirmed from the Brønsted acid complex [(CH3)3PO••H-OAl≡] and Lewis acid complex [(CH3)3PO••Al≡] in the TMPO-loaded H-ZSM-5, respectively. Therefore, it implied that the superacidity characterized by 31 P NMR signals at 85-88 ppm (Fig. 3) on the specific BAS should be derived from the influence of the neighboring tri-coordinated FAL (LAS), generating a Brønsted/Lewis synergy effect. The assignment was also confirmed by the absence of the 31 P signals at 85-88 ppm on a quantitatively TMPO-adsorbed H-ZSM-5 zeolite with a low Si/Al of 140 (see Fig. S15), which has a much less probability of Al-Al pair in close proximity in the zeolite framework due to the much lower Al content.

In order to gain more insights into the Brønsted/Lewis synergy effect, we also performed DFT calculations (see Fig. S16). As shown in the local theoretically optimized structures (Fig. 6b and6c), for TMPO adsorbed on the BAS with a neighboring tri-coordinated FAL site, a 31 P chemical shift of 83 ppm was theoretically predicted, while for TMPO adsorbed on the isolated BAS, a 31 P chemical shift of 76 ppm was achieved, in consistence with our experimental observation. The calculated 31 P and 27 Al isotropic chemical shifts for different adsorption models were listed in Table S3. Except for the unreacted "NMRinvisible" framework LAS (tri-coordinated FAL), the chemical shift parameters extracted from DFT calculations were in consistence with our experimental results. These results indicate that the synergistic effect between framework BAS and framework LAS (tri-coordinated FAL) due to their spatial proximity (Fig. 5e) leads to an enhancement of Brønsted acidity and thus generates the superacidity in H-ZSM-5 zeolite.

Conclusions

Tri-coordinated framework aluminum species can also serve as Lewis acid sites in principle, however, one of fundamental challenges is to observe and distinguish these tri-coordinated framework aluminum species in zeolite. Their relatively low concentration and largely distorted local environment make them "NMR-invisible" in zeolites, thus the experimental study on their acidic property is rarely reported. Benefiting from the current state-of-the-art multi-nuclear and multi-dimensional solid-state NMR spectroscopy, unique insights into the structure and property of tri-coordinated FAL species acting as Lewis acid sites in H-ZSM-5 have been achieved in conjunction with TMPO probe molecule adsorptions. Since the unsaturated tricoordinated FAL species can act with the basic TMPO probe molecules to form distorted tetrahedral FAL sites, two tricoordinated FAL species with different Lewis acidities were discriminated by sensitivity-enhanced 2D 31 P- 27 Al D-HMQC MAS NMR experiments. In particular, the adsorption models of TMPO probe molecules on distinct acidic sites in H-ZSM-5 were clearly illustrated, and the strong 31 P- 27 Al dipolar interaction from the Lewis complex [(CH3)3PO••Al≡] was confirmed by measuring the distance between corresponding 31 P and 27 Al atoms. Furthermore, the remarkably improved resolution and sensitivity of 27 Al MAS NMR spectrum of TMPO-adsorbed H-ZSM-5 zeolite allowed us to obtain a well-resolved 2D 27 Al DQ-SQ homonuclear correlation spectrum, and thus to understand the structure of different framework acid sites in nonhydrated zeolite. Consequently, the superacidity of BAS (with acid strength stronger than 100% H2SO4) induced by the synergy effect between framework BAS and framework LAS (tri-coordinated FAL) was evidenced by their close spatial proximity in dehydrated H-ZSM-5 zeolites and the observation of 31 P chemical shift up to 85-88 ppm for adsorbed TMPO. Notably, apart from the extra-framework Al as LAS in zeolite, the existence of considerable amount of tri-coordinated framework Al serving as LAS in zeolite framework suggests that further studies are needed to understand its catalytic function, which might be essential to rationally design highly efficient zeolite catalysts. The detailed characterizations of tricoordinated framework Al in H-ZSM-5 in this work provide a useful strategy to explore the acidic nature of this type of "moisture-sensitive" species in zeolite catalysts, which may shed light on experimental exploration of its structure-function relationship in zeolite catalysis.

Experimental section

Sample preparations

Parent NH4-ZSM-5 with Si/Al = 25 and Si/Al =140 was purchased from Zeolyst International. The sample was calcined at 773 K under a dry air atmosphere for 10 h to remove the template and the H-form of ZSM-5 was obtained. Prior to adsorption of TMPO, the sample was dehydrated on a vacuum line. The temperature was gradually increased at a rate of 1 °C min -1 and the sample was kept at a final temperature of 673 K at a pressure below 10 -3 Pa overnight. Detailed procedures involved in introducing the TMPO probe molecule onto the sample can be found elsewhere. 45 In brief, a known amount of TMPO dissolved in anhydrous CH2Cl2 was first added into a vessel containing the dehydrated H-ZSM-5 sample in a N2 glovebox, followed by removal of the CH2Cl2 solvent by evacuation on a vacuum line at 323 K or 363 K. To ensure a uniform adsorption of TMPO molecules in the pores/channels of the zeolites, the sealed vessel was further subjected to a thermal treatment at 453 K for 2 h. Finally, the sample was transferred into a ZrO2 MAS rotor in the N2 glovebox prior to the solid-state NMR experiments.

ICP experiment

Atomic emission spectroscopy (AES-ICP) was used to determine the relative amounts of silicon, aluminum and phosphorus in TMPO-loaded H-ZSM-5 samples. Na and Al contents were measured by AES-ICP for a 10 mg sample dissolved in 50 mL HF solution.

FT-IR of pyridine adsorption

FT-IR of pyridine adsorption measurements were conducted on a Bruker Tensor 27 spectrometer. About 20 mg sample was pressed into a self-supported wafer with a diameter of 13 mm. Prior to the measurement, the catalysts were evacuated to 10 -2 Pa at 673 K for 2 h. After cooling down to room temperature, pyridine vapor was introduced into the sample cell at room temperature for 30 min to allow equilibrium, the residual pyridine was removed by vacuum. The FT-IR spectra of pyridineadsorbed samples were measured at 423, 523, and 623 after evacuation for 30 min, respectively.

Solid-state NMR experiments

For all 2D HETCOR experiments, the isotope within the braces represents that it is indirectly detected in F1 dimension. 1D 1 H MAS NMR and 31 P MAS NMR experiments were performed on an 9.4 T Bruker AVANCE-Ⅲ spectrometer using 4 mm rotors at a spinning frequency of 10 kHz (operating at a Larmor frequency of 399.3 MHz and 161.7 MHz, respectively). The 1 H single-pulse acquisition was employed with a pulse width of 3.9 µs (ca. π/2 pulse), a recycle delay of 3 s, and 32 scans. The 31 P single-pulse acquisition was employed with a pulse width of 3 µs (ca. π/2 pulse), a recycle delay of 30 s, and 1024 scans, in which the RF amplitude of 1 H TPPM decoupling was set to ca. 50 kHz.

The 29 Si MAS NMR experiments were performed on an 11.7 T Bruker AVANCE-Ⅲ spectrometer using 7 mm rotors at a spinning frequency of 6 kHz (operating at a Larmor frequency of 99.4 MHz). The 29 Si single-pulse acquisitions were employed with a pulse width of 8.4 µs (ca. π/2 pulse), a recycle delay of 30 s, and 320 scans, in which the RF amplitude of 1 H TPPM decoupling was set to ca. 38 kHz. The single-pulse 27 Al MAS NMR data were collected on an 11.7 T Bruker AVANCE-Ⅲ spectrometer with a 4 mm doubleresonance probe at a spinning rate of 10 kHz, a pulse width of 0.4 µs (ca. π/6 pulse), a recycle delay of 1 s, and 2048 scans. To improve resolution, the single-pulse 27 Al MAS NMR data were also collected on an 18.8 T Bruker AVANCE-Ⅲ spectrometer (corresponding to a 27 Al Larmor frequency of 208.6 MHz) with an 1.9 mm double-resonance probe at a spinning rate of 40 kHz, a pulse width of 0.3 µs (ca. π/9 pulse), a recycle delay of 1 s, and 4096 scans. The 2D 27 Al triple-quantum (3Q) MAS z-filtering NMR experiments with fast amplitude modulation (FAM) enhancement 50 were also carried out at 18.8 T. The pulse duration were set to 3.6 µs for the first strong pulse with νRF ≈ 125 kHz. A FAM pulse train consist of 4 loops of two short pulses (νRF ≈ 125 kHz) and two delays with the same duration of 0.7 µs to convert p = ± 3 to ± 1 coherences. The "soft" π/2 pulses were set to 10.4 µs.

The 1D 31 P { 27 Al} S-RESPDOR and 2D 31 P { 27 Al} PT-D-HMQC MAS NMR experiments were performed on an 11.7 T Bruker AVANCE-Ⅲ spectrometer using commercial 7 mm rotors with an o-ring cap at a spinning frequency of 6 kHz, tuning and matching to 1 H, 31 P and 27 Al Larmor frequencies (500.6 MHz, 202.6 MHz and 130.5 MHz, respectively). 1 H-31 P CP with a contact time of 8 ms was employed to prepare the initial 31 P signal.

For 31 P { 27 Al} S-RESPDOR experiments, a saturation pulse on the 27 Al channel with amplitude of ca. 25 kHz and a length of 166.67 µs (TR) was irradiated at ca. 30 ppm to transfer 31 P- 27 Al interactions. For 2D 31 P { 27 Al} PT-D-HMQC experiments, repetitive sideband-selective (SS) WURST-80 adiabatic pulses with a length 166.67 µs were employed on Al channel during SR4 recoupling to accelerate coherence transfers between 31 P and 27 Al. The values for peak RF amplitude and offset of WURST-80 were optimized with (ν1 max , νoffset ) = (7.5, 175) kHz on 27 Al channel. For 2D 27 Al { 31 P} D-HMQC and 27 Al DQ-SQ MAS NMR experiments, central transition (CT) enhancement of the initial magnetizations ( 27 Al) was obtained by using a SS-WURST-80 irradiation with a length of 1 ms. More details on the NMR parameters of 2D or double-resonance experiments can be found in Table S2.

All simulations of 1D MAS NMR spectra or slices were performed with the DMFIT software. 51 The chemical shift for 1 H, 31 P, 29 Si and 27 Al were referenced to adamantane, 85% H3PO4, kaolinite and 1 M aqueous Al(NO3)3, respectively.

DFT Calculations

All the structures of ZSM-5 zeolite are represented by a 64 T model (containing the complete double 10-membered ring (MR) intersection pores formed by 10-MR straight and 10-MR sinusoidal pore channels), which were extracted from their crystallographic structural data (http://www.iza-structure.org/databases/). The terminal Si-H was fixed at a bond length of 1.47 Å, oriented along the direction of the corresponding Si-O bond. The geometries of TMPO adsorbed at the isolated Brønsted acid site and the Brønsted/Lewis synergetic site of ZSM-5 were optimized over the 64 T cluster model by the combined theoretical ONIOM method. In the calculations, the TMPO and the high layer atoms of the zeolite framework were allowed to fully relax with the dispersion-corrected ωB97XD functional 52 with 6-31G(d, p) basis set, whereas the rest atoms were fixed at their crystallographic positions with AM1 method. Based on the optimized structures, the 31 P and 27 Al NMR chemical shifts were then calculated at the ωB97XD/TZVP level by GIAO method.

Scheme 1 .

 1 Scheme 1. Formation route for tri-coordinated framework Al (FAL) and extra-framework Al (EFAL) species in zeolite.

Figure 1 .

 1 Figure 1. 1D 27 Al MAS NMR spectra of parent NH4-ZSM-5 (a and d), calcined and dehydrated H-ZSM-5 (b and e), TMPO-adsorbed H-ZSM-5 with P/Al = 0.42 (c and f), obtained at B0 = 18.8 T (a-c) and B0 = 11.7 T (d-f). The asterisks denote spinning sidebands.

Figure 2 .

 2 Figure 2. Sheared 2D 27 Al 3QMAS NMR spectra of dehydrated H-ZSM-5 zeolites with different TMPO loadings (a) P/Al = 0.42 and (b) P/Al = 0.18, obtained at B0 = 18.8 T. Skyline projections are drawn along F1 and F2. The determined quadrupolar coupling constants are indicated in the brackets for different Al species.

Figure 3 .

 3 Figure 3. 1D one-pulse 31 P MAS spectra of dehydrated H-ZSM-5 zeolites with different TMPO loadings (a,c) and subsequent rehydration (b,d). And 2D 31 P { 1 H} CP HETCOR MAS NMR spectrum (e) of dehydrated H-ZSM-5 zeolite with a low TMPO loading (P/Al = 0.18).

Figure 4 .

 4 Figure 4. 1D one-pulse 31 P MAS (a) and 2D 31 P { 27 Al} PT-D-HMQC MAS (b) NMR spectra recorded at 11.7 T, and 2D 27 Al { 31 P} D-HMQC MAS NMR spectrum (c) recorded at 18.8 T of dehydrated H-ZSM-5 zeolite (Si/Al = 25) with adsorbed TMPO (P/Al = 0.42). Selected slices with simulations are shown on (d) and (e).

Figure 5 .

 5 Figure 5. 31 P { 27 Al} S-RESPDOR built-up curves of TMPO adsorbed on (a) LAS (tri-coordinated FAL) and (b) BAS of dehydrated H-ZSM-5 zeolite, fitted by analytical formula. DIS and rP-Al represent the dipolar interaction constant and internuclear distance of the 31 P-27 Al spin pair, respectively. And three adsorption structure models of (c) TMPO adsorbed on tri-coordinated FAL, forming a distorted tetrahedral FAL, (d) TMPO adsorbed on Brønsted acid site, and (e) TMPO adsorbed on Brønsted acid site with tri-coordinated FAL in close proximity. The red numbers in the parentheses were the distances obtained by DFT calculations.

Figure 6 .

 6 Figure 6. 2D 27 Al-27 Al DQ-SQ MAS NMR spectra (a) of dehydrated and TMPO-loaded H-ZSM-5 zeolite (P/Al = 0.42), recorded at 18.8 T with a spinning speed of 20 kHz. Local theoretically optimized structures of TMPO adsorbed at the isolated Brønsted site (b) and Brønsted/Lewis acid synergetic site (c) with tri-coordinated framework Al in close proximity. The predicted 31 P chemical shifts of TMPO (blue numbers) and 27 Al chemical shifts of corresponding acid sites (black numbers) were labeled in ppm.

Table 2 .

 2 Concentration of BAS and LAS on H-ZSM-5 (Si/Al = 25) samples.

	Sample	n(BAS+LAS) a μmol/g	n(BAS) b μmol/g	n(LAS) c μmol/g
	H-ZSM-5/25	647.2	511.4	86.1

a calculated by ICP-OES. b determined by 1 H MAS NMR. c calculated from pyridine-FTIR.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grants, 91745111, 21573278, 21733013, 21622311 and 21673283), Bureau of International Cooperation (112942KYSB20180009) and Youth Innovation Promotion Association, The Chinese Academy of Sciences. J. Trébosc and O. Lafon thank the Chevreul Institute (FR 2638) supported by the « Ministère de l'Enseignement Supérieur, de la Recherche et de l'innovation », the « Région Nord-Pas de calais » and the « Fonds Européen de Développement des Régions ». Financial support from the TGIR-RMN-THC Fr3050 CNRS for conducting the research is gratefully acknowledged. We gratefully thank Prof. Christophe Copéret for helpful discussions.

Conflicts of interest

The authors declare no conflict of interest.


Roland Lamberti

Yohan Petetin

A DOUBLE PROPOSAL NORMALIZED IMPORTANCE SAMPLING ESTIMATOR

Keywords: Monte Carlo integration, (self-normalized) importance sampling, variance minimization

Monte Carlo methods are widely used in signal processing for computing integrals of interest. Among Monte Carlo methods, Importance Sampling is a variance reduction technique which consists in sampling from an instrumental distribution and reweighting the samples in order to correct the discrepancy between the target and proposal distributions. When either the target or the proposal distribution is known only up to a constant, the moment of interest can be rewritten as a ratio of two expectations, which can be approximated via selfnormalized importance sampling. In this paper we show that it is possible to improve the self-normalized importance sampling estimate by approximating the two expectations in this ratio via two importance distributions. In order to tune them we optimize the variance of the final estimate under a reasonable constraint. Our results are validated via simulations.

INTRODUCTION 1.Importance Sampling (IS)

Let x be a random variable (r.v.) with density π(.). As far as notations are concerned we do not make a distinction between r.v. and their realizations. In many signal processing problems we are interested in computing the expectation of some function f (x) w.r.t. π:

µ = f (x)π(x)dx = E π (f (x)). (1)
In practice (1) can be very difficult to compute, so one often needs to resort to approximations. In this paper we focus on Monte Carlo (MC) integration methods (see e.g. [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]) which aim at approximating E π (f (x)) with a (possibly weighted) set of random samples. The variance of the estimate depends on the number N of samples, on function f and pdf π, and on the distribution from which the points are sampled from. As far as this last point is concerned, a natural distribution to sample from is the target distribution π, which results in the crude MC estimate μMC = 1 N N i=1 f (x i ) with x i i.i.d.

∼ π. However i) it is generally difficult to sample directly from π;

ii) the crude MC estimate can give poor results, particularly when the regions where π is large do not coincide with those where f is large. IS is a well established technique [START_REF] Kahn | Methods of reducing sample size in Monte Carlo computations[END_REF] [3] [4, §5.4] for bypassing these difficulties. Rewriting [START_REF] Robert | Monte Carlo Statistical Methods[END_REF] as

µ = f (x) π(x) q(x) q(x)dx = E q (f (x) π(x) q(x) ), (2)
where q is some importance distribution, leads to the IS estimator (the dependency in pdf q is stressed via the notation μIS (q)):

μIS (q) = 1 N N i=1 π(x i ) q(x i ) f (x i ), x i i.i.d. ∼ q. (3)
As far as variance reduction is concerned, one can easily show that N var(μ IS (q)) = π 2 f 2 q µ 2 and that the importance density q IS opt which minimizes var(μ IS (q)), and associated (scaled) minimal variance, are

q IS opt (x) ∝ |f (x)|π(x), (4)
N var(μ

IS (q IS opt )) = ( |f (x)|π(x)dx) 2 -µ 2 . ( 5
)
Even if in practice the optimal estimator μIS (q IS opt ) cannot be computed (we know q IS opt up to the unknown constant |f (x)|π(x)dx), this result tells us that the regions where it is important to sample from (whence the term "importance distribution") are not those where π is large, but rather those where |f |π is large (see point ii) above).

Normalized IS (NIS)

From (3) we see that μIS can be computed only if π is known exactly, or π and q are known up to a common constant. If this is not the case one can resort to self-normalized (NIS). More precisely, let π(x) = p(x) p(x)dx where p is known but the normalization constant is not. Then (2) can be rewritten as

µ = f (x) p(x) q(x) q(x)dx p(x) q(x) q(x)dx = E q (f (x) p(x) q(x) ) E q ( p(x) q(x) ) ; (6)
next each expectation can be approximated by a crude MC estimate from the same set of N i.i.d. samples {x i } N i=1 drawn from q. The NIS estimate μNIS (q) of µ reads

μNIS (q) = N i=1 p(x i ) q(x i ) f (x i ) N j=1 p(x j ) q(x j ) = 1 N N i=1 π(x i ) q(x i ) f (x i ) 1 N N j=1 π(x j ) q(x j )
. ( 7)

The first expression in [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] is used to compute the estimate in practice, while the second one is used to analyze its statistical properties.

The properties of μNIS (q) are well understood [5, Chapter 2] [6] [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. First, it is well-established [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF] that under mild assumptions, μNIS (q) a.s.

→ µ;

the bias of μNIS (q) is asymptotically null in 1/N ; its variance can be approximated at first-order by applying the deltamethod to the ratio (7) [8, p.35]:

N var(μ NIS (q)) = var( π(x)f (x) q(x) ) + µ 2 var( π(x) q(x) ) -2µcov( π(x)f (x) q(x) , π(x) q(x) ); (9)
finally the following Central Limit Theorem (CLT) [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF] (see also [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]) holds:

√ N (μ NIS (q) -µ) D → N (0, σ 2,NIS (q)), (10)
σ 2,NIS (q) = E q [ π 2 (x) q 2 (x) (f (x) -µ) 2 ]. (11)
Note that the first-order approximated variance in [START_REF] Veach | Optimally combining sampling techniques for monte carlo rendering[END_REF] and asymptotic variance in [START_REF] Douc | Convergence of adaptive mixtures of importance sampling schemes[END_REF] coincide:

N var(μ NIS (q)) = σ 2,NIS (q).

The choice of the importance distribution q(.) in [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] has been the subject of a vast literature, see e.g. [START_REF] Robert | Monte Carlo Statistical Methods[END_REF][START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF][START_REF] Liu | Monte Carlo Strategies in Scientific Computing[END_REF][START_REF] Veach | Optimally combining sampling techniques for monte carlo rendering[END_REF][START_REF] Cappé | Population Monte Carlo[END_REF][START_REF] Douc | Convergence of adaptive mixtures of importance sampling schemes[END_REF][START_REF] Owen | Monte Carlo theory, methods and examples[END_REF] (and [START_REF] Kroese | Chapter 2 -the cross-entropy method for estimation[END_REF] specifically in the non-normalized case) When the objective is to approximate µ for a large class of functions f , the importance distribution q is often chosen close to π. Now for a specific function f , the IS distribution which minimizes σ 2,NIS (q) in ( 11), and associated minimal asymptotic variance, are [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF]

q NIS opt (x) ∝ |f (x) -µ|π(x), (12)
σ 2,NIS (q NIS opt ) = ( |f (x) -µ|π(x)dx) 2 . ( 13
)

Scope of the paper

In this paper we start from the observation that the minimum variance var(μ IS (q IS opt )) in ( 5) is zero for all positive functions f , while in the normalized case σ 2,NIS (q NIS opt ) in ( 13) is strictly positive whatever function f (except of course in the degenerate situation where f (x) is constant). Heuristically, the reason why is that μNIS approximates the ratio [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF], so both integrals should be estimated as precisely as possible; but even though q NIS opt (x) in ( 12) is the best compromise between |f |p |f |p and p p , it is impossible for one single set of points {x i } N i=1 to be both in regions where |f |π is large, and in regions where π is large (-unless, again, if f is a constant);

In this paper we thus revisit the NIS mechanism by considering two importance distributions q 1 and q 2 , which leads to the Double Proposal Importance Sampling (DPIS) estimator μDPIS (q 1 , q 2 ). We next discuss on how to tune these importance distributions in order to improve the classical estimate μNIS in terms of asymptotic variance while keeping its computational attractiveness. The rest of this paper is organized as follows. Section 2.1 is devoted to the DPIS estimator and its properties. In section 2.2 we discuss the choice of importance distributions q 1 and q 2 and we propose a simple solution. Section 3 deals with simulations and we end the paper with a conclusion.

DOUBLE PROPOSAL IS (DPIS)

The DPIS estimator

In this section we start again from µ = f p p , and we introduce two importance distributions, q 1 for the numerator and q 2 for the denominator. So µ can be rewritten as

µ = f (x) p(x) q1(x) q 1 (x)dx p(x) q2(x) q 2 (x)dx = E q1 (f (x) p(x) q1(x) ) E q2 ( p(x) q2(x) ) . ( 14
)
An MC estimate of µ based on (14) can be obtained with N i.i.d. samples {x i 1 , x i 2 } N i=1 where x i 1 ∼ q 1 and x i 2 ∼ q 2 (note that for the moment, we do not make any assumption on the dependency for a given i between x i 1 and x i 2 ). The resulting estimate μDPIS reads

μDPIS (q 1 , q 2 ) = N i=1 p(x i 1 ) q1(x i 1 ) f (x i 1 ) N j=1 p(x j 2 ) q2(x j 2 ) = 1 N N i=1 π(x i 1 ) q1(x i 1 ) f (x i 1 ) 1 N N j=1 π(x j 2 ) q2(x j 2 )
.

(15) The following proposition describes the statistical properties of μDPIS (q 1 , q 2 ) and generalizes those recalled in Section 1.2 (the proofs are omitted but follow the same line as those of the results presented in Section 1.2).

Proposition 1 (First-order variance and CLT for DPIS) Let x 1 ∼ q 1 and x 2 ∼ q 2 . A first-order approximation of var(μ DPIS ) reads

N var(μ DPIS (q 1 , q 2 ) = var( π(x1)f (x1) q1(x1) ) + µ 2 var( π(x2) q2(x2) ) -2µcov( π(x1)f (x1) q1(x1) , π(x2) q2(x2) ). (16)
In addition, let us assume that

σ 2 (q 1 , q 2 ) = E[( π(x 1 ) q 1 (x 1 ) f (x 1 ) - π(x 2 ) q 2 (x 2 ) µ) 2 ] ( 17
)
is finite. Then the following CLT holds:

√ N (μ DPIS -µ) D → N (0, σ 2 (q 1 , q 2 )) (18)
and the (scaled) first-order approximated variance N × var( μDPIS )(q 1 , q 2 ) in ( 16) coincides with the asymptotic variance σ 2 (q 1 , q 2 ) in (17).

A linear sample adjustment solution

Proposition (1) holds whatever q 1 , q 2 , and the dependency between the samples {x i 1 } drawn from q 1 and {x i 2 } drawn from q 2 . In this section let us discuss how to tune the samples x i 1 and x i 2 . We begin with two simple remarks: i) μDPIS (q 1 , q 2 ) coincides with μNIS when q 1 = q 2 and x i 1 = x i 2 for all i, 1 ≤ i ≤ N ; ii) by contrast with the classical NIS mechanism, the variance of μDPIS (q 1 , q 2 ) reduces to 0 when q 1 = q DPIS 1,opt ∝ |f |p and q 2 = q DPIS 2,opt ∝ p. Of course, neither q DPIS 1,opt nor q DPIS 2,opt are computable in practice. Techniques to approach such distributions have been proposed in the context of unnormalized IS and could be adapted here, but in this paper we follow another direction; we observe that the presence of two (potentially different) importance densities q 1 and q 2 introduces more flexibility as compared to NIS, which leads us to the problem of minimizing the variance term ( 16), but at the price of a reasonable computational effort. To that end, remember that a computational advantage of the NIS mechanism is that one only needs to draw N samples from a single q. In this paper we do want to depart from this unique importance density q = q 1 = q 2 , but as a compromise between sampling diversity and sampling cost we will assume that x i 2 is determinist given x i 1 . So let us now consider an importance distribution q and N initial i.i.d. samples {x i } N i=1 drawn from q. Our basic idea is to transform the samples from q into samples {x i 1 } N i=1 from q 1 in order to approximate the numerator of (14), and next into samples {x i 2 } N i=1 from q 2 in order to approximate the denominator of (14). Here we will focus on linear transformations of the initial samples {x i } N i=1 relying on scaling parameters (α 1 , α 2 ) and translation parameters (β 1 , β 2 ):

x i 1 = α 1 x i + β 1 , x i 2 = α 2 x i + β 2 . ( 19
) It means that collections {x i 1 } N i=1 and {x i 2 } N i=1 are drawn i.i.d. from distributions q 1 (•) = | 1 α1 |q( •-β1 α1 ) and q 2 (•) = | 1 α2 |q( •-β2
α2 ), respectively but are dependent. The linear relationship between x i 1 and x i 2 for all i is given by x i 2 = α2 α1 (x i 1β 1 ) + β 2 . It remains to compute the optimal parameters (α 1 , α 2 , β 1 , β 2 ). Under this constraint (16) can be expressed as

N var(μ DPIS ) = 1 C 2 (E q [ 1 q 2 (x) (|α 1 |p(α 1 x + β 1 ) × f (α 1 x + β 1 ) -|α 2 |p(α 2 x + β 2 )µ) 2 ] ( 20
)
where C = p(x)dx. Note that if the minimization is performed correctly, the resulting estimate μDPIS necessary outperforms μNIS (provided N is large enough to satisfy the firstorder approximation ( 16)) which corresponds to the setting

α 1 = α 2 = 1, β 1 = β 2 = 0.
Equation ( 20) cannot be computed exactly because of the µ factor and of the integrals, but can be approximated via the available samples {x i } N i=1 initially drawn from q. An MC estimate of (20) reads

N var(μ DPIS ) = 1 C 2 ( 1 N N i=1 1 q 2 (x i ) × (|α 1 |p(α 1 x i + β 1 )f (α 1 x i + β 1 )-|α 2 |p(α 2 x i + β 2 )μ NIS ) 2 ). ( 21
)
To minimize this approximated variance, we use a subgradient descent method with backtracking line search for determining parameter move steps, using multiple random initializations over a given area in an effort to avoid local minima. Moreover, in order to improve the overall approximation we update our initial estimation of µ (needed in the first-order variance expression) at every gradient iteration using samples obtained with the current set of parameters.

SIMULATIONS

In this section we compare, given an importance distribution q, the classical NIS estimate μNIS (q) based on N samples {x i } N i=1 from q, to our DPIS estimate μDPIS (q 1 , q 2 ) in which the original samples x i are linearly adjusted by minimizing (20) or (21). The estimates are compared in terms of RMSE over P = 100 MC runs.

A toy example

Let us first compute an MC estimate of µ in (1) where f (x) = x 2 and π(x) = N (x, 0, 1) (N (x, m, σ 2 ) denotes the Gaussian distribution with mean m and variance σ 2 ). For the initial importance distribution we set q = π. Here var(μ DPIS ) can be computed in closed form as a function of (α 1 , β 1 , α 2 , β 2 ) because µ = 1 is known and the integrals are computable. More precisely, N var(μ DPIS ) reads

α 2 1 2α 2 1 -1 e β 2 1 2α 2 1 -1 [ β 4 1 (2α 2 1 -1) 4 + 6β 2 1 α 2 1 (2α 2 1 -1) 3 + 3α 4 1 (2α 2 -1) 2 ] + α 2 2 2α 2 2 -1 e β 2 2 2α 2 2 -1 - 2|α 1 ||α 2 | α 2 1 + α 2 2 -1 × exp( β 2 1 (1 -α 2 2 ) + β 2 2 (1 -α 2 1 ) + 2β 1 β 2 α 1 α 2 2(α 2 1 + α 2 2 -1)
)

× [ α 2 1 α 2 1 + α 2 2 -1 + β 2 α 1 α 2 -β 1 α 2 2 + β 1 α 2 1 + α 2 2 -1 2 ]. (22)
The optimal parameters are (α opt

1 ≈ 2.2, β opt 1 ≈ 0, α opt 2 ≈ 1.25, β opt 2 ≈ 0)
which means that only the variance of q is adjusted.

We compute two estimates based on a double proposal distribution: the first one uses the optimal parameters above while the second one uses the parameters (α 1 , α 2 , β 1 , β 2 ) estimated from (21) via a sub-gradient method (the parameters are estimated for each MC run). We compare these two DPIS estimates to an NIS estimate optimized using the same method (i.e., only one density q is optimized, which corresponds to setting α 2 = α 1 and β 2 = β 1 ), and finally to an NIS estimate with density q.

Fig. 1 displays the RMSE of the four estimates as a function of the number of samples N . We first observe that our first-order variance minimization produces translated and scaled samples {x i 1 , x i 2 } N i=1 that are better adapted to the estimation of µ than either the original samples {x i } N i=1 drawn from the importance distribution q, or the samples {x i } N i=1 adjusted from {x i } N i=1 in order to fit simultaneously the numerator and the denominator; in other words, adjusting two densities, one for the numerator and one for the denominator, is significantly better than adjusting one single density to both. Next, even if the parameters (α 1 , α 2 , β 1 , β 2 ) are deduced from an MC approximation of the first-order variance of μDPIS (see ( 21)), the corresponding estimate gives almost the same performance as that which uses the optimal parameters. This shows that our technique has the potential to remain relevant in more complex models where the integrals are not available in closed-form. In order to illustrate this claim, Section 3.2 presents a simple example of such a model. Fig. 1: Non-linear Gaussian model, f (x) = x 2 , π(x) = q(x) = N (x, 0, 1). Simulation with 100 MC runs comparing DPIS estimators, using either exact or approximated integrals in the first-order variance expressions for minimization, with classical NIS and optimized classical NIS.

Bayesian, non-linear Gaussian setting

We now consider a Bayesian problem where x is hidden, y is observed and the target distribution π(x) = p(x|y) ∝ p(x)p(y|x) with p(x) = N (x, 0, 1), p(y|x) = N (y, x 2 , 1) and we want to estimate the moment of f (x) = N (x, 3, 1).

In this setting, in order to compute RMSE values we take as a reference a classical NIS estimate using 10 7 samples. Fig. 2 displays the classical NIS estimate based on i.i.d. samples drawn from q, μDPIS where the parameters have been estimated by the optimization method described in §2. 4. CONCLUSIONS In this paper we observed that NIS approximates the ratio of two expectations with one single set of random samples, which by construction cannot fit simultaneously the numerator and the denominator. We thus proposed an NIS technique based on the use of two instrumental distributions. In order to tune these proposals in practice we adjusted an initial set of samples drawn from some initial distribution q via two linear transforms, one targeting the numerator and the other the denominator. These linear transforms are optimized by minimizing the first-order variance of our DPIS estimate. Our simulations show that the resulting DPIS estimate outperforms NIS estimates which take into account only one set of samples, even if this set is adjusted according to the DPIS procedure. Possible perspectives for future work include improving the variance approximation provided by the delta method, which is only of the first order and may thus fail in certain models if N is insufficiently high, and assessing the empirical relevance of our optimized proposal density adjustment in more complex, high-dimensional, models.

Fig. 2 :

 2 Fig.2: Bayesian non-linear Gaussian model, f (x) = N (x, 3, 1), π(x) = q(x) = N (x, 1, 1), p(y|x) = N (y, x 2 , 1). Simulation with 100 MC runs comparing the DPIS estimator using approximated variance integrals with classical NIS and optimized classical NIS.

  2, and an NIS estimate optimized over only two parameters (α 2 = α 1 and β 2 = β 1 ). The parameters obtained in this scenario are (α 1 ≈ 1.25, β 1 ≈ 1, α 2 ≈ 2.45, β 2 ≈ -0.1), on average.

		7	×10 -4			Classical NIS
						Optimized classical NIS
		6				Optimized DPIS with approximation
		5
	RMSE	4
		3
		2
		0.4 1	0.6	0.8	N (number of particles) 1.4 1.2 1	1.6	1.8	2 ×10 4


Ming Yin

Mohammed Naidjate
email: mohammed.naidjate@polymtl.ca

Nicolas Bracikowski
email: nicolas.bracikowski@univ-nantes.fr

Antoine Pierquin

Didier Trichet

Topology Optimization of Magnetic Actuator based on Reluctance Network Modeling and Adjoint Variable Method

Keywords: Topology optimization, magnetic actuator, reluctance network analysis, adjoint variable method

Topology Optimization (TO) has great significance in primarily the concept design of a product. In the electrical engineering community, conventional topology optimization is usually based on a spatial discretization that also is used as the mesh for Finite Element Analysis (FEA). In this paper, we proposed to use a local equivalent circuit modeling method -Reluctance Network Analysis (RNA). Afterward, a gradient-based optimization algorithm, line search method, was chosen as the optimizer where the sensitivity information was calculated by Adjoint Variable Method (AVM). Finally, the feasibility of this model for topology optimization was verified through a case study of a magnetic actuator.

I. INTRODUCTION

T OPOLOGY OPTIMIZATION has received wide atten- tion from electrical engineers since the seminal work by Bendsøe and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF]. Contrary to parametric optimization methods, the optimization has no initial constraint, so, it may result in an entirely novel topology that can't be imagined beforehand.

Many different approaches, such as Homogenization Method, Density Method, ON/OFF Method, Level-set Method and etc., were developed for TO [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF]. In the electromagnetic field, the Density Method was the most widely adopted for TO. A typical case is the so-called Solid Isotropic Material with Penalization (SIMP) which introduces an interpolation function that interprets the material properties as continuous variables. Design space is discretized into elements, then, different predefined materials are allocated into them.

To process TO, a numerical model is required for characterizing the performances. There are several approaches for modeling in electrical engineering, and the two popular used are the finite element analysis (FEA) and the reluctance network analysis (RNA). Their different characteristics were analyzed in [START_REF] Wang | A Review on Magnetic Gears: Topologies, Computational Models, and Design Aspects[END_REF]. This paper presents an effort to implement reluctance network modeling. In order to verify the feasibility of this model in TO, a classical problem of magnetic actuator, as shown in Fig. 2, is chosen as the case study due to its simple structure and magnetic path. It consists of three parts, the core (with 'C' shape), the coil (in the yellow regions), and the armature (the moving part). The chosen design region is marked out with a dashed line box as shown in Fig. 2 and the optimization objective is to maximize the magnetic force imposed on the armature. The density method was adopted and the material in the design space was updated toward the steepest descent direction. To speed up the solving, the sensitivity information usually is obtained through AVM. Many applications of this method are based on FEA, in this paper, the AVM will be developed for RNA cases.

II. RELUCTANCE NETWORK MODEL

The physical electromagnetic domain is decomposed into elements via a mesh as shown by the dashed line with four vertices in Fig. 1. In this case, each element can be represented by a magnetic equivalent circuit consisting of four lumped reluctances (R l , R r , R u , R d ) and two flux sources (φ l , φ r ) with a structure similar to Fig. 1. Every reluctance is expressed by ( 1) or (2), while the flux sources are expressed by [START_REF] Wang | A Review on Magnetic Gears: Topologies, Computational Models, and Design Aspects[END_REF]. The reluctance values depend on the discretization of the region to be studied and actuator dimensions. The coefficient N is chosen to respect Ampère's theorem in each loop of the circuit.

R l = R r = 1 µ 0 µ r 1 z ∆x/2 ∆y (1)
R u = R d = 1 µ 0 µ r 1 z ∆y/2 ∆x ( 2
)
φ l = φ r = N I R l = N I R r (3)
Where µ 0 is the vacuum permeability, µ r is the relative permeability of the physical material, z is the depth of the device, and ∆x, ∆y are the step size of the mesh.

Therefore, the whole system can be represented by an equivalent network. We apply nodal analysis to every node of the circuit. This equation for the whole system can be expressed by [START_REF] Korman | On the Use of Topology Optimization for Synchronous Reluctance Machines Design[END_REF], where the [Λ] is the matrix that contains the reluctances of the system, the {ε} is the magnetic potential vector and the {φ} is the flux sources vector due to current in the winding.

By solving the matrix system the values of magnetic potential in each node are obtained. Other information such as magnetic potential difference and magnetic flux can also be deduced.

[Λ]{ε} = {φ} (4)

III. SENSITIVITY COMPUTATION

The magnetic force (F) applied on the armature is computed through Maxwell stress tensor method. For the density method, the permeability of the predefined material is defined as [START_REF] Okamoto | Investigation of topology optimization of magnetic circuit using density method[END_REF] where ρ is the normalized density and the exponential term n of this mapping function is named penalization coefficient. By setting the different values of n, we would get different interpolation functions and its influence on the optimization result was investigated in the paper [START_REF] Korman | On the Use of Topology Optimization for Synchronous Reluctance Machines Design[END_REF]. In this paper, the line search was chosen as the optimization algorithm and the element's material was updated toward the steepest descent direction. Therefore, gradient information is needed. Usually, the sensitivity can be obtained through the finite difference method or AVM. The main drawback of the finite difference method is that it is computationally expensive for TO which gets large amounts of design variables [START_REF] Okamoto | Investigation of topology optimization of magnetic circuit using density method[END_REF]. In comparison, the AVM needs only an extra calculation to get the sensitivities of all the variables.

Firstly, differentiate a normalized density in the design region, ρ i , on both sides of the system equation. Then we would get an expression of (6), where the derivative of {φ} to ρ equals 0, cause the flux source is independent of the design variables. The objective function, F, can be seen as an implicit function, F (ε, ρ), about ε and ρ. A similar differential equation of the objective function to each variable i gets the form (7). Combine the two formulas, (6) and (7) ,finally we get a new expression (8). The introduced variable λ, is the adjoint variable. We could get the sensitivity of all the design variables by determining the value of adjoint variable once time.

µ = µ 0 (1 + (µ r -1)ρ n ) (5) ∂[Λ] ∂ρ i {ε} + [Λ] ∂{ε} ∂ρ i = 0 (6) dF dρ i = ∂ F (ε, ρ) ∂ρ i + ∂ F (ε, ρ) ∂ε T ∂ε ∂ρ i (7) dF dρ i = ∂ F (ε, ρ) ∂ρ i - ∂ F (ε, ρ) ∂ε T [Λ] -1 λ T ∂[Λ] ∂ρ i {ε} (8)
IV. RESULTS

The system in this case study is assumed to be linear. The optimized topology as shown in Fig. 2 seems logical because it reduces the flux leakage and links more flux to the armature. In order to remove the intermediate materials (0< ρ <1), a filter is required. Among all different filtering functions, the feasibility factor (FF) is adpoted in this paper. Fig. 3 depicts the Fig. 2. Optimal design and the used mesh Fig. 3. Force evolution evolution of magnetic force with iterations. The force reaches its relatively high value around the 8th iteration, but a clear boundary of the topology in the design space appears until around the 35th iteration. Regarding of time consumption, the optimization based on AVM gets considerable time saving compared with the finite difference method. As a whole, the topology optimization based on the RNA model is feasible in this case and it seems that it could be applied to more complex cases such as electrical machines, etc. Magnetic saturation effects and volume constraints would be discussed in the final paper.

Fig. 1 .

 1 Fig. 1. The mesh element and its associated magnetic equivalent circuit