Literature DB >> 31316062

Acoustic meta-atom with experimentally verified maximum Willis coupling.

Anton Melnikov1,2,3,4, Yan Kei Chiang5, Li Quan6, Sebastian Oberst7, Andrea Alù6,8, Steffen Marburg9, David Powell10.   

Abstract

Acoustic metamaterials are structures with exotic acoustic properties, with promising applications in acoustic beam steering, focusing, impedance matching, absorption and isolation. Recent work has shown that the efficiency of many acoustic metamaterials can be enhanced by controlling an additional parameter known as Willis coupling, which is analogous to bianisotropy in electromagnetic metamaterials. The magnitude of Willis coupling in a passive acoustic meta-atom has been shown theoretically to have an upper limit, however the feasibility of reaching this limit has not been experimentally investigated. Here we introduce a meta-atom with Willis coupling which closely approaches this theoretical limit, that is much simpler and less prone to thermo-viscous losses than previously reported structures. We perform two-dimensional experiments to measure the strong Willis coupling, supported by numerical calculations. Our meta-atom geometry is readily modeled analytically, enabling the strength of Willis coupling and its peak frequency to be easily controlled.

Entities:  

Year:  2019        PMID: 31316062      PMCID: PMC6637156          DOI: 10.1038/s41467-019-10915-5

Source DB:  PubMed          Journal:  Nat Commun        ISSN: 2041-1723            Impact factor:   14.919


Introduction

Acoustic metamaterials[1,2] have demonstrated unique material properties which do not exist naturally, such as negative bulk modulus[3] and negative dynamic mass density[4]. These properties have enabled the development of acoustic superlenses[5-7], barriers[8], cloaking devices[9], and the enhancement of non-linear effects[10]. Metamaterials are typically arrays of sub-wavelength structures, known as meta-atoms, with geometry engineered to control their dynamic mass and stiffness. It has recently been shown that more efficient metamaterial designs can be created by incorporating an additional degree of freedom, represented by the Willis coupling parameter. Willis coupling is a term in the acoustic and elastic constitutive relations that couples potential and kinetic energy[11-13], analogous to the bianisotropy parameter in electromagnetism[14,15]. The Willis coupling and bianisotropy parameters can only be non-zero in structures with low symmetry[16], for example, in one dimension a structure lacking mirror symmetry in the propagation direction is required[17]. The inclusion of these terms into the constitutive relations has been shown to resolve violations of causality and passivity in metamaterial homogenization[18,19]. Recent work has demonstrated that the incorporation of Willis coupling or bianisotropy into metamaterial structures of sub-wavelength thickness, known as metasurfaces, can improve their efficiency when refracting at large angles[20-22]. While bianisotropy has been demonstrated and engineered in a wide range of electromagnetic meta-atom designs[23], approaches for controlling the degree of Willis coupling in acoustic meta-atoms are not well-established. Recently, a bound on the maximum value of the Willis coupling parameter was derived, based on the conservation of energy[21]. It was shown how meta-atoms can be designed to reach this theoretical bound, using space-coiling structures with long meander-line channels[21,24-26]. While such structures are advantageous for achieving resonance in a sub-wavelength volume, they are difficult to manufacture reproducibly, typically requiring additive manufacturing techniques. Moreover, their channel widths are often comparable to the viscous and thermal boundary layer thickness, and their channel lengths are of the order of the wavelength, leading to high thermo-viscous losses, and a significant reduction in scattering efficiency[27-29]. It was shown numerically in ref. [21] that the thermo-viscous losses in space-coiling meta-atoms can reduce their Willis coupling magnitude to be significantly lower than the theoretical bound. Experimental evidence of Willis coupling has been reported in both one-dimensional[2,17] and two-dimensional metamaterial structures[20]. It has also been shown how the effective medium properties of a bulk metamaterial incorporating Willis coupling can be derived from the polarizability of individual acoustic meta-atoms[30]. However, the theoretical limit on the strength of Willis coupling has not been tested experimentally, and it remains unknown how closely this limit may be approached in practice. To resolve this question we propose an acoustic meta-atom, which is designed to realize maximum Willis coupling and to minimize thermo-viscous losses. Experimental results obtained from a fabricated sample are compared to numerical calculations, showing good agreement of the resonant frequency and lineshape, with a reduction in magnitude that we attribute primarily to the energy leakage through the top and bottom waveguide plates, which are not perfectly hard. The simplicity of the structure enables us to present an analytical model for its polarizability, showing how the Willis coupling can be tailored to have any value between zero and the theoretical bound.

Results

Meta-atom design

Acoustic wave interaction with a meta-atom is conveniently described by its polarizability tensor. Due to the sub-wavelength size of meta-atoms, their scattering is dominated by monopole and dipole moments, and their polarizability can be written as where M is the scalar monopole moment, D is the vector dipole moment, is the polarizability tensor, and are the incident pressure and the velocity at the center of the meta-atom[21]. The off-diagonal terms and represent the Willis coupling between the dipolar and monopolar responses. We show how the reported structures exhibiting Willis coupling[2,17,21] can be replaced by a simpler, and more reproducible structure, by avoiding thin channels and large areas of fluid-structure interfaces. The structure can be tuned to achieve any value of Willis coupling up to the theoretical bounds and can be readily analyzed using a closed form analytical solution which assists in understanding the physical mechanisms behind its operation. Achieving a strong acoustic polarizability in a small volume requires a resonant structure. Inspired by resonant sonic crystals[8] and Helmholtz resonators with multiple apertures[31-33], our novel meta-atom design is presented in Fig. 1. This meta-atom exhibits Willis coupling due to the asymmetrical neck openings. The air within each neck is treated as an incompressible mass, while the air in the cavity acts as a spring, together creating an oscillating system excited by an incident acoustic wave. As the oscillation occurs in the fluid domain only, this avoids wave coupling through fluid-structure interfaces. Since the structure avoids long and thin channels, thermo-viscous losses are expected to be greatly reduced.
Fig. 1

Meta-atom geometry and dimensions. Here a is the cylinder radius and ri is the cavity radius. The neck widths w are in general different for each aperture, the neck length l is common to all apertures and the cavity volume

Meta-atom geometry and dimensions. Here a is the cylinder radius and ri is the cavity radius. The neck widths w are in general different for each aperture, the neck length l is common to all apertures and the cavity volume Peak Willis coupling is expected to occur close to the Helmholtz resonator’s eigenfrequency, which is dependent on the air mass moving within the apertures and the inner cavity volume. In two dimensions (2D), the moving mass is determined by the aperture cross-sections A = w and the neck length l, which is equal for all apertures due to the inner and outer cylindrical boundaries being concentric. The 2D-volume is determined by the inner radius ri = a – l, where a is the outer radius. By neglecting radiation damping, the eigenfrequency for N apertures can be approximated as (see Supplementary Note 3) A more accurate model including radiation damping is developed subsequently in this work. Generally, an arbitrary number of apertures may be included. However, for controlling Willis coupling along one axis, two oppositely arranged apertures are sufficient, as illustrated in Fig. 1. The apertures can be used to adjust the resonant frequency and the level of Willis coupling. The maximum Willis coupling magnitude is exhibited when the structure is shown in Fig. 1 becomes maximally asymmetric, consistent with the single aperture meta-atoms shown in ref. [21].

Experimental verification

To experimentally demonstrate Willis coupling in the presented meta-atom, its polarizability is determined in a 2D experiment. A single aperture meta-atom having maximum asymmetry is manufactured and investigated with the dimensions a = 20 mm, ri = 10 mm, and w = 12 mm (see Fig. 2a). The incident and scattered pressure fields are measured in a 2D parallel-plate waveguide[29] (see Fig. 2b) with subsequent extraction of the polarizability tensor as detailed in Supplementary Note 5.
Fig. 2

Experimental system and results. a Stainless steel sample: single aperture meta-atom with a = 20 mm, ri = 10 mm, w = 12 mm, and h = 66 mm. b, c Schematic and photograph of the experimental set-up. d–g Theoretically predicted and experimentally determined components of the normalized polarizability tensor. The error bars show the standard deviation resulting from the least squares fit over 12 incident angles

Experimental system and results. a Stainless steel sample: single aperture meta-atom with a = 20 mm, ri = 10 mm, w = 12 mm, and h = 66 mm. b, c Schematic and photograph of the experimental set-up. d–g Theoretically predicted and experimentally determined components of the normalized polarizability tensor. The error bars show the standard deviation resulting from the least squares fit over 12 incident angles The polarizability tensor defined in Eq. (1) has elements with different units and values differing by many orders of magnitude. Therefore it is convenient to introduce the normalized polarizability tensor ′, the elements of which are shown for our meta-atom structure in Fig. 2d–g. The normalized values are defined as , , , and  [21]. As required by reciprocity[21], the tensor satisfies ′ = ′, since the off-diagonal terms are equal to each other with a sign reversal (see Fig. 2e, f). The error bars in Fig. 2d–g show the standard deviation of the experimentally extracted polarizability terms, obtained from the least squares fit (see Methods section). To verify that the magnitude of Willis coupling is close to the theoretical maximum, Fig. 3 shows the numerically and analytically calculated and the experimentally determined and . For the numerical analysis, the 2D Boundary Element Method (BEM) has been applied. Here, the peak Willis coupling at k ⋅ a = 0.75 achieves 90% of the theoretical bound 4ω−2.
Fig. 3

Willis coupling magnitude. Result for a single aperture meta-atom obtained from BEM simulation (cyan solid line), polarizability theory (blue dashed line) and experimentally (red solid and black dashed lines), normalized to the theoretical bound 4ω−2

Willis coupling magnitude. Result for a single aperture meta-atom obtained from BEM simulation (cyan solid line), polarizability theory (blue dashed line) and experimentally (red solid and black dashed lines), normalized to the theoretical bound 4ω−2

Polarizability theory

To show how the Willis coupling and other polarizability components can be tailored by adjusting the meta-atom geometry, we develop a polarizability theory. The predictions of this theory are shown in Figs. 2 and 3. It is based on multiple aperture Helmholtz resonator dynamics for the outward directed particle displacement ξ within each aperture (see Supplementary Note 3) Here is the external pressure at aperture n, ρ0 is the medium density, c is the speed of sound, A = w is the cross-section of the aperture, K = c2ρ0 is the bulk modulus, V is the inner cavity volume, is the effective neck length (including a radiative correction) and is the radiative loss coefficient (see Fig. 1). The ω3-term corresponds to radiative losses, the ω2-term represents Newton’s second law applied to the air in the neck, and the summation over m accounts for coupling between apertures via compression of the center cavity using Hooke’s law. Considering all apertures results in a matrix equation Keqξ = pext, where Keq is the dynamic stiffness matrix given by Supplementary Equation (38). This matrix allows the displacement to be solved for an arbitrary incident pressure field and hence, the contribution of the oscillating air masses (see Fig. 1) to the scattering of the meta-atom. See Supplementary Note 3 for full details of the derivation. We consider now the shape illustrated in Fig. 1. This restricts the Willis coupling to a single axis and simplifies the problem to a 2 × 2 matrix. In this case, the resonator polarizability tensor res is obtained as The ratio of aperture widths determines the strength of Willis coupling and can be used to optimize the structure. Expanding Eq. (4) for a single aperture meta-atom with w2 = 0 leads to the following expression for the polarizability The dynamic stiffness matrix becomes a scalar equation with the remaining projection matrix being singular. The singularity arises because both the monopole and dipole moments are determined from the air movement within a single aperture, but it presents no computational difficulties, since the inverse of the polarizability tensor is not required. As the Helmholtz resonator is embedded within a cylinder (see Fig. 1), there is an additional influence on the meta-atom polarizability due to the background scattering from the cylinder. Considering only the dipole and monopole terms, the polarizability tensor of a cylinder of radius a is (see Supplementary Note 4)where J is a Bessel function, is a Hankel function of the first kind and is the wavenumber. As expected for a symmetrical geometry, the off-diagonal terms corresponding to Willis coupling are zero. Since the resonator and cylinder are superimposed, they influence the effective incident fields of each other through an additional scattered term: This results in a coupled formulation for the monopole and dipole moments of the cylinder and the Helmholtz resonator aswhere is the incident field vector and represents the acoustic propagation from the cylinder to each of the apertures. Finally, adding the monopole and dipole moments from Eq. (8) results inwhich gives the total polarizability tensor tot. The components of tot for a single aperture meta-atom with a = 20 mm, ri = 10 mm, w = 12 mm, c = 343 m ⋅ s−1, and ρ0 = 1.2 kg ⋅ m−3 are shown together with the experimental results in Fig. 2.

Discussion

It has been shown based on passivity and reciprocity conditions[21] that maximum Willis coupling can be achieved only if the polarizability components share the same magnitude and the meta-atom has no losses. Figure 4 shows the magnitude of the experimentally extracted polarizability components, where a crossing point can be observed at k⋅a = 0.75. At this point the magnitudes are very close to each other, giving a precise indication of the frequency of maximum Willis coupling.
Fig. 4

Experimental polarizability. Amplitudes of experimentally determined polarizability components showing shared magnitude closely to k · a = 0.75

Experimental polarizability. Amplitudes of experimentally determined polarizability components showing shared magnitude closely to k · a = 0.75 The lossless condition cannot be perfectly satisfied, since any fluid medium exhibits thermal and viscous losses, which are greatly magnified close to the boundaries. Furthermore, our two-dimensional parallel-plate waveguide shown in Fig. 2c will leak some energy through the top and bottom plates, since they can only approximate the perfectly hard boundary condition. Although these losses are not treated in our theory, Fig. 2d–g shows that the experimental line-shape and resonant frequency are well described by our theoretical model. To illustrate how closely our meta-atom approaches the theoretical bound, the magnitude of the Willis coupling is plotted in Fig. 3. The slight frequency shift observable in Fig. 2d–g between the experimental and theoretical values can be attributed to the thermo-viscous losses of air. To investigate these mechanisms within the meta-atom, we model the thermo-viscous losses using the Finite Element Method (FEM), with the resulting Willis coupling shown in Fig. 5. The results reveal a downshifting of the frequency by 2.2% and show that thermo-viscous losses lead to a reduction of 0.32% in the magnitude of Willis coupling. We note that this reduction in magnitude is much less than that previously reported for space-coiling meta-atoms with thin channels, which was ~21% at k ⋅ a ≈ 0.75 (see Supplementary Material of ref. [21]). In Fig. 5 we also plot the Willis coupling magnitude for the space-coiling meta-atom from ref. [21] where the much higher radiative quality factor of the space-coiling structure leads to higher internal dissipation, hence a much greater impact of thermo-viscous losses on the Willis coupling.
Fig. 5

Sensitivity to thermo-viscous losses. Numerical comparison of Willis coupling showing the influence of thermo-viscous losses in air, which cause a reduction in magnitude of only 0.32% for the c-shape meta-atom, whereas for the space-coiling structure reported in ref. [21] the magnitude drops by ∼21%

Sensitivity to thermo-viscous losses. Numerical comparison of Willis coupling showing the influence of thermo-viscous losses in air, which cause a reduction in magnitude of only 0.32% for the c-shape meta-atom, whereas for the space-coiling structure reported in ref. [21] the magnitude drops by ∼21% In addition to the experimentally demonstrated maximum Willis coupling, the structure presented in Fig. 1 can be tailored to have Willis coupling from zero up to the theoretical bound. To demonstrate this property, normalized Willis coupling for four different parameter sets is shown in Fig. 6. The red meta-atom on top illustrates the single aperture configuration similar to the experimentally investigated structure from Fig. 2a. Once a second aperture with neck width w2 is introduced, Willis coupling is significantly reduced (Fig. 6 magenta and violet lines). This would result in a frequency shift, as expected from Eq. (2). To avoid this, w1 is tuned to match the peak frequency of the single aperture case. When w2 is further increased and the shape starts to converge to the symmetrical case (Fig. 6 blue, ), the Willis coupling becomes very weak and disappears completely when . This behavior is of practical importance, since it allows tailoring of the Willis coupling to any desired values. A full parametric analysis of the influence of w1 and w2 on the resonant frequency and peak Willis coupling is presented in Supplementary Note 6.
Fig. 6

Control of Willis coupling. Willis coupling of four different meta-atom geometries, where w2 is varied to control Willis coupling and w1 is tuned to keep the peak frequency fixed

Control of Willis coupling. Willis coupling of four different meta-atom geometries, where w2 is varied to control Willis coupling and w1 is tuned to keep the peak frequency fixed In conclusion, we introduced and experimentally validated a novel meta-atom exhibiting strong Willis coupling and providing a low radiative Q factor. The experiment revealed a Willis coupling magnitude reaching ~90% of the theoretical bound. In this structure, thermo-viscous losses in air are quite small, whereas they are much stronger in previously reported space-coiling meta-atoms. Additionally, the simple shape of our structure facilitates manufacturing and enables accurate analytical modeling. Combining the models of a Helmholtz resonator and a cylindrical scatterer, a theory was developed and shown to agree well with numerical simulations. Since our structure enables Willis coupling to be tailored, this theory can be used to engineer Willis coupling for specific applications.

Methods

Extraction of polarizability tensor

The polarizability relationship given in Eq. (1) is used to illustrate Willis coupling, where it appears as the off-diagonal terms and . An extraction method is necessary to obtain the polarizability of a scatterer from an experiment or numerical model. For simplicity, only the 2D case is considered. We build on the method for extracting polarizability from highly symmetric 2D structures in ref. [29] which does not account for Willis coupling. This method makes use of the incident and scattered pressure fields around the object and fits Bessel and Hankel functions to them as and . Since only monopole and dipole components are of importance, the problem can be reduced to considering terms. Following this, the incident pressure at the meta-atom center is and the particle velocity . The monopole moment and the dipole moments can be retrieved from the scattered field. For further information see Supplementary Note 1 and 2. The expansion coefficients β and γ can be obtained from measured or numerically determined pressure on circles with radii Rinc and Rscat respectively. From the orthogonality of exponential functions, the coefficients can be found asand To avoid the singularity of Eq. (10) due to zeros of the Bessel function[29], we ensure . However, Rscat should be significantly larger than the meta-atom outer radius to reduce near field contributions. These conflicting requirements mean that Rinc and Rscat are different in general. To fully determine the polarizability in 2D, , , M and D must be determined for at least 3 incident angles. The incident field quantities for all available angles θ1.. are arranged in a matrix as Knowing M and D for each θ1.., allows the polarizability tensor α to be determined by inversion of . For increased robustness, we take additional angles. The polarizability tensor is then determined via least squares as

Numerical model

To obtain a numerical solution for the polarizability a custom 2D BEM code is used, as described in ref. [29]. It uses continuous elements with quadratic interpolation functions[34] and discretization by collocation method[35] with an adaptive integration scheme[36]. Initially the solids are treated as acoustic hard boundaries. Thermo-viscous losses in air are not included in this formulation. The incident field is a unit intensity plane wave, therefore and are known explicitly. The medium density and the speed of sound are set to ρ0 = 1.2 kg ⋅ m−3 and c = 343 m ⋅ s−1. To calculate the influence of thermo-viscous losses in air on the polarizability of the meta-atom (Fig. 5) 2D FEM calculations are performed with the COMSOL Multiphysics Thermoviscous Acoustics Module. The acoustic boundary of the meta-atom is treated as rigid. For the lossless case, the mechanical boundary condition of the structure is set to be slip and the thermal boundary is set to be adiabatic. Both the viscosity coefficient and the thermal conductivity coefficient are set as zero in the simulation. For the lossy case, the mechanical boundary is set to be no slip and the thermal boundary is set to be isothermal. The thermal and viscous coefficients used for air are ρ0 = 1.2043 kg · m−3, c = 343.14 m · s−1, μ = 1.814 × 10−5 Pa · s, μ = 1.0884 × 10−5 Pa · s, k = 0.025768 W · m−1 · K−1, C = 1005.4 J · kg−1 · K−1, α = 0.0034112 K−1, γ = 1.4 and β = 9.8692 × 10−6 Pa−1. The finite element type is a triangular element with quadratic interpolation functions. The maximum element size is set to approximately 28 elements per wavelength at 2500 Hz.

Waveguide scattering experiment

To experimentally determine the incident and scattered pressure fields we use the 2D anechoic waveguide chamber presented in ref. [29]. The propagation medium is air with an acoustic velocity of 343 m · s−1. The height of the chamber (66 mm) supports only a single propagation mode at frequencies up to 2600 Hz, making it a 2D wave propagation system. The excitation source is a speaker, excited by a continuous wave with frequency varying between 1500 Hz and 2500 Hz in 50 Hz steps. The response is measured by a microphone, which is moved in two axes by belts driven by stepper motors. The sample is a single aperture meta-atom shown in Fig. 1a. This meta-atom was manufactured from stainless steel with precisely machined surface with maximum surface roughness of 4 μm (Rz4). Additionally, rubber seals (black) are glued on top and bottom to prevent air leakage from the resonator cavity and to achieve more homogeneous clamping. The incident field is measured at a radius Rinc = 40 mm. The scattering of the sample is measured at 12 different incident angles (0°, 30°, 60°, 90°, 120°, 150°, 180°, 210°, 240°, 270°, 300°, 330°) at a radius of Rscat = 200 mm
  14 in total

1.  Matryoshka locally resonant sonic crystal.

Authors:  Daniel P Elford; Luke Chalmers; Feodor V Kusmartsev; Gerry M Swallowe
Journal:  J Acoust Soc Am       Date:  2011-11       Impact factor: 1.840

2.  Ultra-sparse metasurface for high reflection of low-frequency sound based on artificial Mie resonances.

Authors:  Y Cheng; C Zhou; B G Yuan; D J Wu; Q Wei; X J Liu
Journal:  Nat Mater       Date:  2015-08-31       Impact factor: 43.841

3.  Focusing ultrasound with an acoustic metamaterial network.

Authors:  Shu Zhang; Leilei Yin; Nicholas Fang
Journal:  Phys Rev Lett       Date:  2009-05-15       Impact factor: 9.161

4.  Broadband acoustic cloak for ultrasound waves.

Authors:  Shu Zhang; Chunguang Xia; Nicholas Fang
Journal:  Phys Rev Lett       Date:  2011-01-10       Impact factor: 9.161

5.  Radiative feedback in Helmholtz resonators with more than one opening.

Authors:  Hans Günter Dosch
Journal:  J Acoust Soc Am       Date:  2016-11       Impact factor: 1.840

6.  Maximum Willis Coupling in Acoustic Scatterers.

Authors:  Li Quan; Younes Ra'di; Dimitrios L Sounas; Andrea Alù
Journal:  Phys Rev Lett       Date:  2018-06-22       Impact factor: 9.161

7.  Quasi-phase-matched backward second-harmonic generation by complementary media in nonlinear metamaterials.

Authors:  Li Quan; Xiaozhou Liu; Xiufen Gong
Journal:  J Acoust Soc Am       Date:  2012-10       Impact factor: 1.840

8.  Metagratings: Beyond the Limits of Graded Metasurfaces for Wave Front Control.

Authors:  Younes Ra'di; Dimitrios L Sounas; Andrea Alù
Journal:  Phys Rev Lett       Date:  2017-08-10       Impact factor: 9.161

9.  Experimental evidence of Willis coupling in a one-dimensional effective material element.

Authors:  Michael B Muhlestein; Caleb F Sieck; Preston S Wilson; Michael R Haberman
Journal:  Nat Commun       Date:  2017-06-13       Impact factor: 14.919

10.  Systematic design and experimental demonstration of bianisotropic metasurfaces for scattering-free manipulation of acoustic wavefronts.

Authors:  Junfei Li; Chen Shen; Ana Díaz-Rubio; Sergei A Tretyakov; Steven A Cummer
Journal:  Nat Commun       Date:  2018-04-09       Impact factor: 14.919

View more
  4 in total

1.  Odd Willis coupling induced by broken time-reversal symmetry.

Authors:  Li Quan; Simon Yves; Yugui Peng; Hussein Esfahlani; Andrea Alù
Journal:  Nat Commun       Date:  2021-05-10       Impact factor: 14.919

2.  Microacoustic Metagratings at Ultra-High Frequencies Fabricated by Two-Photon Lithography.

Authors:  Anton Melnikov; Sören Köble; Severin Schweiger; Yan Kei Chiang; Steffen Marburg; David A Powell
Journal:  Adv Sci (Weinh)       Date:  2022-04-24       Impact factor: 17.521

3.  Asymptotics of the meta-atom: plane wave scattering by a single Helmholtz resonator.

Authors:  M J A Smith; P A Cotterill; D Nigro; W J Parnell; I D Abrahams
Journal:  Philos Trans A Math Phys Eng Sci       Date:  2022-10-10       Impact factor: 4.019

4.  Cloaking In-Plane Elastic Waves with Swiss Rolls.

Authors:  Younes Achaoui; André Diatta; Muamer Kadic; Sébastien Guenneau
Journal:  Materials (Basel)       Date:  2020-01-17       Impact factor: 3.623

  4 in total

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