Literature DB >> 30441788

Influence of Oxygen Vacancy Behaviors in Cooling Process on Semiconductor Gas Sensors: A Numerical Analysis.

Jianqiao Liu1, Wanqiu Wang2, Zhaoxia Zhai3, Guohua Jin4, Yuzhen Chen5, Wusong Hong6, Liting Wu7, Fengjiao Gao8.   

Abstract

The influence of oxygen vacancy behaviors during a cooling process in semiconductor gas sensors is discussed by the numerical analysis method based on the gradient-distributed oxygen vacancy model. A diffusion equation is established to describe the behaviors of oxygen vacancies, which follows the effects of diffusion and exclusion in the cooling process. Numerical analysis is introduced to find the accurate solutions of the diffusion equation. The solutions illustrate the oxygen vacancy distribution profiles, which are dependent on the cooling rate as well as the temperature interval of the cooling process. The gas-sensing characteristics of reduced resistance and response are calculated. Both of them, together with oxygen vacancy distribution, show the grain size effects and the re-annealing effect. It is found that the properties of gas sensors can be controlled or adjusted by the designed cooling process. The proposed model provides a possibility for sensor characteristics simulations, which may be beneficial for the design of gas sensors. A quantitative interpretation on the gas-sensing mechanism of semiconductors has been contributed.

Entities:  

Keywords:  diffusion equation; gas sensor; numerical analysis; oxygen vacancy; semiconductor

Year:  2018        PMID: 30441788      PMCID: PMC6264112          DOI: 10.3390/s18113929

Source DB:  PubMed          Journal:  Sensors (Basel)        ISSN: 1424-8220            Impact factor:   3.576


1. Introduction

The successful invention of ZnO gas sensors by Seiyama in 1962 started the new era of the application of semiconductor in a gas detecting field [1]. From then on, several kinds of semiconductors were put into practice to develop advanced gas sensors, which contained a variety of sensing materials such as SnO2 [2,3,4], ZnO [5,6], and WO3 [7] as well as in different forms such as bulk [8], thick film [9,10], and thin film [11,12,13]. Recently, semiconductor quantum dots were also introduced into the investigation of gas sensors [14,15,16]. The gas-sensing mechanism of semiconductors was understood by Morrison who concluded that gas detection was completed by the change in resistance of a sensor [17]. A specific region near the surface of a semiconductor grain, which is usually called a depletion layer, would decrease its resistance when exposed to reducing gases while the resistance would increase if oxidizing gases were introduced. This transducing process was achieved by the reaction between reducing gas and adsorbed oxygen or the competitive adsorption between adsorbed oxygen and oxidizing gas [18]. In 2009, Yamazoe concluded the gas-sensing mechanism of semiconductor into three levels including the receptor function, the transducer function, and the utility factor [19]. Among these levels, the receptor function was the most complex one and described how a single semiconductor grain responded to a stimulant gas. It included several decisive factors such as the grain size and shape [20,21,22], surface composite [23,24], the type of adsorbed oxygen [25,26], and the characteristics of oxygen vacancies [27,28]. The oxygen vacancy (V) plays a vital role because it ionizes and provides semi-conductive nature to the metal-oxide system [29], which would be insulative without any defect inside. The behaviors of oxygen vacancies are essential questions about the gas-sensing mechanism of semiconductors. The formation of V usually takes place during sintering in which the oxygen atom escapes from the crystal lattice according to Equation (1). Then, some of V are involved into the grain, which grows up during the sintering. The V on the grain surface act as adsorption site of adsorbed oxygen while the ones inside the grain provide quasi-free electrons for conductance after their ionization. Furthermore, the defects would migrate in the grain if the temperature is above the absolute zero K especially when the operating temperature of a semiconductor gas sensor is usually above 200 °C. The importance of V was recognized and there were some specific studies on it [30,31,32]. However, in the theoretical investigations, the roles of V were simplified since its density had to be assumed as a uniform distribution throughout the grain [17,18]. This assumption was good for calculation, but it deviates from the actual condition. Therefore, a gradient-distributed oxygen vacancy model was proposed based on the experimental influences of the cooling rate on the gas-sensing characteristics of SnO2 thin films [33,34]. The migration of oxygen vacancies in a semiconductor grain is concluded as diffusion and exclusion. The former one results from the density gradient of defects in the grain, which follows Fick’s second law. The latter one is an excluding tendency, which is provided by the cooling process because the crystal lattice is stable if it is free of defects [35,36]. Thus, the diffusion equation of oxygen vacancies in the cooling process is established, as seen in Equation (2). In this scenario, N(r,t) represents the density of oxygen vacancies at a place of r and time of t in an ideal spherical grain in which the sphere coordinates are established as Figure 1. Considering the symmetry of sphere coordinates, only one-dimensional model is applied.
Figure 1

Schematic drawing of the migration of oxygen vacancies under the diffusion and exclusion effects in a sphere coordinates established in an ideal grain.

The first term of the right side of the diffusion equation indicates the diffusion effect, which is expressed by Fick’s second law. D(t) is the temperature-dependent diffusion coefficient and it is formulated as Equation (3). In this case, D0 is the pre-exponent constant and E is activation energy of diffusion. k is the Boltzmann constant. T0 is the initial temperature when the cooling process starts. β is the cooling rate and t is the time elapsed in the cooling process. The other term of the right side of the diffusion equation denotes the exclusion effect. P is the possibility of the exclusion for a defect moving outwards to a nearby position in unit time and it can be formulated as Equation (4). In this case, ν0 is the thermal vibration frequency of the oxygen atom. E represents the activated energy of defect migration and E0 is the unit energy decrease of the system, which resulted from the one-step exclusion of defects. Therefore, the diffusion equation of Equation (2) can be rewritten as Equation (5). The analytical solution of the diffusion equation is not available. Thus, several presumptions had to be made in the previous works [31,33,34] to simplify the calculation, which simulated the V distribution and gas-sensing characteristics of the semiconductor grain. However, some of the presumptions were far from the practical situations, which leads to an inaccuracy in the solutions that describe the distributions of oxygen vacancies. In the present work, an attempt on finding the accurate solutions is made by using numerical analysis methods. The numerical solutions that follow the diffusion equation of Equation (5) are found to describe the kinetics of oxygen vacancies and their time-dependent distributions. Furthermore, they are used to evaluate the electrical characteristics of semiconductor gas sensors such as resistance and response to reducing gases. The influences of the grain size effects, the re-annealing effect, and the controlled temperature interval on the sensor properties are also discussed.

2. Methods

The computational tool of MATLAB is used to find the numerical solutions of the diffusion equation of Equation (5). A discrete N(r,t) is established in the two-dimensional array of N(i,j), which derives from N(r,t) = N(iΔr,jΔt). In this scenario, Δr and Δt are the infinite small fragments in space and time. Therefore, Equations (3) and (4) can be transformed into the discrete expressions of Equations (6) and (7) after two one-dimensional arrays of D(j) and P(j) are established. The partial derivative of N(r,t) against r and t can be expressed as Equations (8) and (9) after discretization. Thus, the diffusion equation of Equation (5) can be transformed into Equation (10) from which Equation (11) can be easily obtained. Equation (11) tells how the oxygen vacancy density at a point interacts with nearby positions after a time of Δt elapses under the effects of diffusion and exclusion, which is shown in Figure 2. If N represents the V density on the grain surface, the discrete initial condition of Equation (12) and boundary conditions of Equation (13) are considered.
Figure 2

Schematic drawing of the calculation process of oxygen vacancy density in a semiconductor grain during the cooling process.

It is noted that N is also a time-dependent variable once the cooling process starts. The density of oxygen vacancy on the grain surface would change when V interacts with aerial oxygen, according to Equation (1). However, specific studies on this topic are still lacking. Therefore, a presumption is made that N is of linear dependence on time since N(t) = N × (1 + k). In this scenario, k is a constant indicating the increasing rate of V density on grain surface. Therefore, the numerical solution can be obtained by a MATLAB calculation, which provides possibilities to discuss the relationships among several important parameters that determine the gas-sensing characteristics of semiconductors. As described by the Poisson’s law, the potential at grain boundaries is correlated to the space charge density in the depletion layer with the abrupt model applied [37], which is shown in Equation (14), provided that the oxygen vacancy is assumed to be first-order ionized. A discrete one-dimensional array of V(k) is established by letting V(x) = V(kΔx) where Δx is the infinite small fragment in space. Thus, Equation (14) can be rewritten into the discrete expression of Equation (15). The boundary conditions can be expressed if w denotes the width of the depletion layer as Equation (16). By using k = w/Δx, Equation (16) transforms to Equation (17), which is used to calculate the potential barrier height together with Equation (15) when k falls into [0, k], as shown in Figure 3.
Figure 3

Schematic drawing of the calculation process of the potential barrier in the depletion layer of a semiconductor grain based on the density distribution of oxygen vacancies.

The reduced resistance (R/R0) of the sensor can be calculated after the potential barrier height at the grain surface (qV) is obtained by qV = qV(0) (see Equation (18)). In this case, R0 is the flat-band resistance. The sensor response (S) is defined as the ratio of the resistance in air (R) to the one in a reducing target gas (R), as S = R/R. The response is stimulated by the target gas, which may consume the adsorbed oxygen on the grain surface and release electrons back to the depletion layer. This leads to a decrease in the depletion layer width. If α is introduced to indicate the concentration of the target gas, the change of the depletion layer width can be expressed by Equation (19) in the one-dimensional model [38], which is shown in Figure 4. Furthermore, α denotes the percentage of the seized electrons that released back to the depletion layer from the adsorbed oxygen.
Figure 4

Potential barrier in the depletion layer of a semiconductor grain (a) in air and (b) in reducing gas.

Therefore, the response (S) can be calculated based on the grain resistance before and after exposure to target gas. In the calculations, the parameters are set as follows: D0 = 0.0431 m2/s [39], E = 2.7 eV [39], ν0 = 1014 s−1 [40], E − E + E0 = 0.1 eV [34], N = 2 × 1025 m−3 [18], k = 3 × 10−5 s−1, α = 0.5, and w = 4 nm [41,42].

3. Results

3.1. Oxygen Vacancy Distribution

The influence of the cooling rate on the gradient distribution profile of oxygen vacancies in a 25 nm semiconductor grain is shown in Figure 5. The oxygen vacancies have various distribution profiles, which are determined by the cooling rate. For the quenched sample with the largest cooling rate of 3600 °C/h, the initial oxygen vacancies do not get sufficient time to migrate. Therefore, they are frozen at the place where they are at the starting of the cooling process, which leaves an almost uniform distribution throughout the grain. For the slowly-cooled sample, the initial oxygen vacancies migrate under the effect of diffusion and exclusion, which forms a gradient distribution in the grain. The gradient is of negative dependence on the cooling rate.
Figure 5

Influence of the cooling rate on the gradient distribution profile of oxygen vacancies in the semiconductor grain with a radius of 25 nm.

Figure 6 shows the influence of the cooling rate on the steady state oxygen vacancy density at the center of the grain (N) where the defect density is the smallest throughout the grain. It controls the density difference between the surface and center. The influence of the cooling rate on the average oxygen vacancy density in the depletion layer (N) is also indicated in Figure 6. Within the cooling rate concerned, there are two regions where the vacancy densities have the linear correlation with the cooling rate. One is β = 1–200 °C/h and the other is β = 400–1000 °C/h. It infers that the V gradient and the amount in the depletion layer can be easily controlled by the cooling rate in the two regions above. Considering the determination of the oxygen vacancy amount in the depletion layer on the sensor performances, it is possible to control the sensor properties by adjusting the cooling rate in the fabrication process.
Figure 6

Influence of cooling rate on steady state oxygen vacancy density at the grain center (N) and the average oxygen vacancy density in the depletion layer (N) of a 25 nm grain.

The transient distribution of oxygen vacancies in a semiconductor grain during the cooling process is illustrated in Figure 7, according to the numerical solution of Equation (5) in which R = 25 nm and β = 100 °C/h. The initial uniform V distribution is driven to distribute in a gradient profile. It describes that the V distribution experiences three stages in the cooling process: (1) at the starting stage, the V density is uniform throughout the grain, (2) then V migrates under the effects of diffusion and exclusion at the transient stage, and, (3) lastly, the distribution reaches the steady state in a gradient profile.
Figure 7

Time-dependent oxygen vacancy distribution in a 25 nm semiconductor grain from the initial uniform distribution to a steady state gradient distribution.

3.2. Gas-Sensing Characteristics

The relationships between gas-sensing characteristics and cooling rate at various operating temperatures are shown in Figure 8 and Figure 9 in which the reduced resistance (R/R0) and response (S) show similar performances. They descend with a cooling rate at the beginning and then reach the lowest point at the cooling rate of 400 °C/h. After that, the slight increases are observed until the sensor properties remain constant when the cooling rate is above 1800 °C/h.
Figure 8

The relationship between reduced resistance (R/R0) of the gas sensor and the cooling rate at the operating temperatures of 25, 100, and 200 °C.

Figure 9

The relationship between the sensor response to reducing gas and the cooling rate at the operating temperatures of 25, 100, and 200 °C.

The transient states of the gas-sensing characteristics of the semiconductor gas sensor during the cooling process are simulated in Figure 10. Both of the reduced resistance and response to reducing gas are increasing during the cooling process with the V distribution driven from a uniform distribution to a gradient one, as described in Figure 5. However, sharp declines are observed at the start of the cooling process.
Figure 10

Transient states of gas-sensing characteristics of a semiconductor gas sensor during the cooling process.

The calculation results are compared with the experimental results in order to validate the applicability of the gradient-distributed oxygen vacancy model. The experimental sensor resistance and response are introduced from the previous study [34], which describes the preparation details of the SnO2 thin film gas sensors and its dependence of sensing performances on the cooling rate in the fabrication process. The actual results are plotted in Figure 11, which also illustrates the calculated influence of the cooling rate on sensor properties. In this case, R0 is assumed to be 0.5. It is observed that the experimental plots located around the calculated correlations show good agreement. Therefore, the calculation in the present work has the potential to interpret the gas-sensing mechanism of semiconductor gas sensors.
Figure 11

Calculation and experimental results of the sensor properties against the cooling rate of 25–1800 °C/h.

3.3. The Grain Size Effects

The steady-state distributions of oxygen vacancies in the semiconductor grains with various grain sizes from 1 to 60 nm are illustrated in Figure 12. The distribution profiles show a significant grain size effect. In a 60 nm grain, the density difference is rather large between the surface and the center. However, less than 1% difference in V density is found in a grain with R = 1 nm. This means that a large grain can maintain a large internal V density difference.
Figure 12

Steady-state distribution of oxygen vacancies in semiconductor grains with a radius of 1 to 60 nm.

The grain size also makes impacts on N and N. As shown in Figure 13, both N and N are of negative dependence with grain radius. The V density at the grain center indicates the gradient of distribution profile. The V amount in the depletion layer controls the supply of electrons, which results from the ionization of V. It is concluded that the supply of the electrons in the depletion layer determines the amount of adsorbed oxygen on the grain surface and, therefore, decides the gas-sensing characteristics [18,25]. Figure 13 shows that N keeps almost constant when the grain radius is larger than 30 nm. However, it increases with the reducing grain size when R < 30 nm. This will lead to the grain size effects of the sensing performances of semiconductor gas sensors.
Figure 13

The grain size effect on the steady-state oxygen vacancy density at the grain center (N) and average steady-state oxygen vacancy density in the depletion layer (N).

The grain size effects of calculated sensor properties at various operating temperatures are shown in Figure 14. When the grain radius is approaching the depletion layer width, both of the reduced resistance and response to reducing gas have significant grain size effects that are in agreement with the experimental observation [20,43]. The comparison between calculation results and experimental responses is shown in Figure 15. The experimental plots are extracted from C. Xu’s report [43] while the calculation is completed by using the following parameter settings: N = 5 × 1025 m−3, α = 0.27 (H2), and 0.37 (C4H10) [42] and T = 300 °C (H2) and 400 °C (C4H10). The experimental grain size effects are well described by the calculation results. Furthermore, if the grain radius is equal to or smaller than the depletion layer width, the grain falls into volume depletion. In this case, both the reduced resistance and response decrease dramatically with the descending grain size until they reach unit when R = 1 nm. For a volume-depleted grain, there are few free electrons inside it. The flat-band resistance R0 is rather large because of the descending of the Fermi level, which results in R/R0 to approach 1. On the other hand, the shortage of the electron supply prevents the adsorption of oxygen on the grain surface, which intercepts the gas-sensing mechanism. Therefore, a tiny calculated response is obtained. It infers that optimized gas-sensing properties could be obtained when the grain radius equals the depletion layer width by either controlling the grain size or adjusting the depletion layer width.
Figure 14

The grain size effects of reduced resistance and response at various operating temperatures of 25–200 °C.

Figure 15

The grain size effects of calculation results and experimental response of the SnO2 gas sensors to 800 ppm H2 at 300 °C and 800 ppm C4H10 at 400 °C, which are extracted from C. Xu’s work [43].

3.4. Re-Annealing Effect

The re-annealing effect is found in the previous study and it mentions that the quenched grain would have the same property as the normal sintered grain after a second annealing is conducted [33]. Figure 16 shows the calculated V distribution profiles in semiconductor grains with the radius of 25 nm. The quenched grain with a cooling rate of 1800 °C/h has an approximately uniform V distribution. After that, a re-annealing with a cooling rate of 25 °C/h is conducted and a gradient V distribution is obtained. The result is the same as the one in the normal sintered grain with a cooling rate of 25 °C/h. It infers that the V migration is frozen by quenching and restarted by re-annealing.
Figure 16

The steady-state oxygen vacancy distribution profiles in a 25-nm grain with normal sintering (β = 25 °C/h), quenching (β = 1800 °C/h), and re-annealing (β = 25 °C/h).

As shown in Figure 17 and Figure 18, after sintered with various cooling rates from 25 to 1800 °C/h, the grains are treated by a re-annealing process with a cooling rate of 25 °C/h. No matter what the cooling rates are used with the grain during the first time, the re-annealing process will lift the resistance and response to almost the same level as the one with a cooling rate of 25 °C/h. Although the experimental gas sensor properties appear as non-ideal plots when compared with the calculation, both of them show a similar tendency of the re-annealing effect.
Figure 17

The re-annealing effect of experimental and calculated resistance of semiconductor gas sensors.

Figure 18

The re-annealing effect of experimental and calculated response of semiconductor gas sensors.

3.5. Controlled Temperature Interval

The results above show that the temperature setting in the cooling process may make impacts on the V distribution as well as gas-sensing characteristics of semiconductor gas sensors. Thus, a series of designed cooling processes are simulated in order to investigate the influence of the temperature interval on the gas sensors. The temperature difference of the cooling process is set to be 300 °C while the starting and end temperatures are controlled, which is shown in Figure 19. The cooling rate is set to be 100 °C/h. In a semiconductor grain with the radius of 25 nm, the distribution profiles are dependent on the temperature interval. High temperature would drive V to appear in a profile of a large gradient. However, after treatment in a low temperature, the grain has a uniform V distribution in the center and gradient distribution near the surface. Figure 20 shows the influences of the temperature interval on N and N. Furthermore, the temperature interval can be used to control the gas-sensing characteristics of the semiconductor gas sensors. Figure 21 illustrates the dependence of reduced resistance and response on the temperature intervals. Both of the gas-sensing properties benefit from the low temperature treatment.
Figure 19

Oxygen vacancy distribution profiles in semiconductor grains with the radius of 25 nm are fabricated by using designed cooling processes, which have the same temperature difference of 300 °C but various starting and end temperatures.

Figure 20

Influences of the temperature interval on the steady state oxygen vacancy density at the grain center (N) and an average steady state oxygen vacancy density in the depletion layer (N).

Figure 21

Influences of the temperature interval on gas-sensing characteristics of semiconductor gas sensors.

4. Discussion

In the sections above, the numerical analysis was carried out to find the numerical solutions of the diffusion equation of V in the semiconductor grain. The solutions were used to illustrate the V distribution and also to simulate the gas-sensing properties. The experimental sensor resistance and response were employed to check the validity of the calculation results, which shows good agreement. However, some considerations need to be discussed. The gradient-distributed oxygen vacancy model was proposed to explain the influence of the cooling rate on the gas-sensing characteristics of semiconductor gas sensor. A diffusion equation of V migration was established in a one-dimensional model based on the diffusion and exclusion effects. In the previous studies [33,34,44], the parameter of T (end temperature of cooling process) had to be used to find the analytical solutions of the diffusion equation. This approximation is helpful to mathematical conduction but leads to an inaccuracy in the calculation results. The model cannot describe the V behaviors during the cooling process precisely. The imperfection drives us to implement the numerical analysis method in the research of semiconductor gas sensors. As described in Equation (5), the diffusion equation is able to simulate the actual cooling process in the gas sensor fabrication. The discrete expression of the diffusion equation is established to find the numerical solutions based on boundaries conditions from the initial condition to the steady state via a step-by-step process. The numerical solutions provide a precise description of the V behaviors in the cooling process and the gas-sensing characteristics. Furthermore, the influences of the cooling rate on the gas sensor properties are concluded in this work. Compared with previous research [33,34,44], the present results provide a more precise description of V behaviors during the cooling process. Therefore, it is possible to calculate the gas-sensing properties of the sensors prior to practical fabrication and the calculation results are beneficial for the design of sensor preparation. In the calculations, there are several presumptions that need to be considered. One of them is the value of N and its time dependence. N indicates the V density on the grain surface. It is a temperature-dependent parameter and is influenced by many other factors such as V formation and annihilation, according to Equation (1), as well as diffusion and exclusion effects based on Equation (5). However, the research on N is still lacking. Thus, a presumption has to be made that N is of linear dependence on time. A modification could be made based on the specific studies on N. The second one is the tunneling effect, which is not taken into consideration in the present study. Tunneling may take place in the semiconductor inside a grain or between grains. It provides possibilities for oxygen atoms and electrons to migrate through the potential barrier. Therefore, the diffusion equation of Equation (5) and the Poisson equation of Equation (14) should be amended. The third presumption is that the present calculation uses the same flat-band resistance (R0) when the grain is in volume depletion. In this situation, the reduced resistance (R/R0) and response (S) decrease rapidly to unit because of the lowered Fermi level, which would lead to a change in R0. However, it is not mathematically included in the calculation. A quantitative study for a volume-depleted grain is expected. Although there are some approximations in the calculation, the present work provides an opportunity to understand the gas-sensing mechanism of the semiconductor in a quantitative manner. The validity of the calculation is checked by experimental results. The proposed model has probable applications in the sensor design by simulating the fabrication process as well as calculating the sensor characteristics. Some modifications could be made to the model if the specific research studies are involved. Then, it can be used to simulate the actual sensor characteristics not only in the fabrication process but also in the working circumstance. Therefore, it is possible to interpret the gas-sensing mechanism of the semiconductor and also some experimental phenomena, which are not fully understood such as aging and the degeneration of sensors.

5. Conclusions

The behaviors of V in the cooling process of semiconductor gas sensors are described by the diffusion equation based on the gradient-distributed oxygen vacancy model. The numerical analysis method is used to find the accurate solutions, which illustrate the V migration, the time-dependent distributions, and their determination on the gas-sensing characteristics of sensors. Several conclusions have been drawn. The V distribution profile in semiconductor grain is determined by the cooling rate in the cooling process. Quenching or fast cooling can freeze the oxygen defects at the place where they form, which may results in an almost uniform distribution. A sintering or re-annealing with a low cooling rate may lead to a gradient V distribution. For the V density at the grain center and average density in the depletion layer, there are two regions where the defects densities have linear correlations with the cooling rate. The gas-sensing characteristics of the semiconductor correlate to the cooling rate. Negative relationships are observed for the reduced resistance and response to reducing gas against the cooling rate below 400 °C/h. The transient states of the gas-sensing properties during the cooling process show that there are sharp declines of reduced resistance and response at the starting of the cooling process. A designed re-annealing process can adjust the properties of semi-conductor gas sensors. The V distribution and gas sensor properties appear to have significant grain size effects. Large grains may maintain a great gradient of V distribution. When the grain radius approaches the depletion layer width, the gas-sensing properties increase until they reach the maximum. After the grain is in volume depletion, the sensor loses its gas-sensing properties rapidly. The optimized gas-sensing performances appear when R = w. The cooling process in the fabrication can be designed for controlling the sensor performances, which are dependent on the temperature interval. The sensor resistance and response benefit from the low temperature treatment, which reduces the V distribution gradient in the semiconductor grain.
  1 in total

Review 1.  Inhomogeneous Oxygen Vacancy Distribution in Semiconductor Gas Sensors: Formation, Migration and Determination on Gas Sensing Characteristics.

Authors:  Jianqiao Liu; Yinglin Gao; Xu Wu; Guohua Jin; Zhaoxia Zhai; Huan Liu
Journal:  Sensors (Basel)       Date:  2017-08-10       Impact factor: 3.576

  1 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.