首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 31 毫秒
1.
The tendency of today’s fuel injection systems to reach injection pressures up to 3000 bar in order to meet forthcoming emission regulations may significantly increase liquid temperatures due to friction heating; this paper identifies numerically the importance of fuel pressurization, phase-change due to cavitation, wall heat transfer and needle valve motion on the fluid heating induced in high pressure Diesel fuel injectors. These parameters affect the nozzle discharge coefficient (Cd), fuel exit temperature, cavitation volume fraction and temperature distribution within the nozzle. Variable fuel properties, being a function of the local pressure and temperature are found necessary in order to simulate accurately the effects of depressurization and heating induced by friction forces. Comparison of CFD predictions against a 0-D thermodynamic model, indicates that although the mean exit temperature increase relative to the initial fuel temperature is proportional to (1  Cd2) at fixed needle positions, it can significantly deviate from this value when the motion of the needle valve, controlling the opening and closing of the injection process, is taken into consideration. Increasing the inlet pressure from 2000 bar, which is the pressure utilized in today’s fuel systems to 3000 bar, results to significantly increased fluid temperatures above the boiling point of the Diesel fuel components and therefore regions of potential heterogeneous fuel boiling are identified.  相似文献   

2.
The thermodynamic driving force for domain growth in a rank-2 laminated ferroelectric crystal is derived in this article, and we used it, together with a homogenization theory, to study the issue of enhanced electrostrictive actuation recently reported by Burcsu et al. [2004. Large electrostrictive actuation of barium titanate single crystals. J. Mech. Phys. Solids 52, 823-846]. We derived this force from the reduction of Gibbs free energy with respect to the increase of domain concentration. It is shown that both the free energy and the thermodynamic force consist of three parts: the first arises from the difference in M0 and M1, the linear electromechanical compliances of the parent and product domains, respectively, at a given level of applied stress and electric field, the second stems from the electromechanical work associated with the change of spontaneous strain and spontaneous polarization during domain switch, and the third from the internal energy due to the distribution of polarizations strain and electric polarization inside the crystal. We prove that the first term is substantially lower than the second one, and the third one is identically zero with compatible domain pattern. The second one is, however, not exactly equal to the commonly written sum of the products of stress with strain, and electric field with polarization during switch, unless both domains have identical moduli in the common global axes. We also show that, with compatible domain patterns and when M1=M0, this driving force is identical to Eshelby's driving force acting on a flat interface due to the jump of energy-momentum tensor. Applications of the theory to a BaTiO3 crystal subjected to a fixed axial compression and decreasing electric field from the [0 0 1] state reveal that the crystal undergoes a three-stage switching process: (i) the 0→90° switch to form a rank-1 laminate, (ii) the 0→180° switch inside the 0° domain to form a laminate I with a concurrent 90°→−90° switch inside the 90° domain to form laminate II, creating a rank-2-laminated domain pattern, and (iii) finally the 90→180° switch. It is the exchange of stability between the 0, 90°, and 180° domains under compression and electric field that is the origin of the enhanced actuation. We illustrate these intrinsic features by showing the evolution of these domains, and demonstrate how the reported large actuation strain can be attained with a rank-2 laminate.  相似文献   

3.
Different measures have been suggested for quantifying the amount of fragmentation in randomly compacted crushable aggregates. A most effective and popular measure is to adopt variants of Hardin's [1985. Crushing of soil particles. J. Geotech. Eng. ASCE 111(10), 1177-1192] definition of relative breakage ‘Br’. In this paper we further develop the concept of breakage to formulate a new continuum mechanics theory for crushable granular materials based on statistical and thermomechanical principles. Analogous to the damage internal variable ‘D’ which is used in continuum damage mechanics (CDM), here the breakage internal variable ‘B’ is adopted. This internal variable represents a particular form of the relative breakage ‘Br’ and measures the relative distance of the current grain size distribution from the initial and ultimate distributions. Similar to ‘D’, ‘B’ varies from zero to one and describes processes of micro-fractures and the growth of surface area. However, unlike damage that is most suitable to tensioned solid-like materials, the breakage is aimed towards compressed granular matter. While damage effectively represents the opening of micro-cavities and cracks, breakage represents comminution of particles. We term the new theory continuum breakage mechanics (CBM), reflecting the analogy with CDM. A focus is given to developing fundamental concepts and postulates, and identifying the physical meaning of the various variables. In this part of the paper we limit the study to describe an ideal dissipative process that includes breakage without plasticity. Plastic strains are essential, however, in representing aspects that relate to frictional dissipation, and this is covered in Part II of this paper together with model examples.  相似文献   

4.
We develop a continuum-mechanical formulation and generalization of the Navier–Stokes-α equation based on a recently developed framework for fluid-dynamical theories involving higher-order gradient dependencies. Our flow equation involves two length scales α and β. The first of these enters the theory through the specific free-energy α 2|D|2, where D is the symmetric part of the gradient of the filtered velocity, and contributes a dispersive term to the flow equation. The remaining scale is associated with a dissipative hyperstress which depends linearly on the gradient of the filtered vorticity and which contributes a viscous term, with coefficient proportional to β 2, to the flow equation. In contrast to Lagrangian averaging, our formulation delivers boundary conditions and a complete structure based on thermodynamics applied to an isothermal system. For a fixed surface without slip, the standard no-slip condition is augmented by a wall-eddy condition involving another length scale characteristic of eddies shed at the boundary and referred to as the wall-eddy length. As an application, we consider the classical problem of turbulent flow in a plane, rectangular channel of gap 2h with fixed, impermeable, slip-free walls and make comparisons with results obtained from direct numerical simulations. We find that α/β ~ Re 0.470 and /h ~ Re −0.772, where Re is the Reynolds number. The first result, which arises as a consequence of identifying the specific free-energy with the specific turbulent kinetic energy, indicates that the choice β = α required to reduce our flow equation to the Navier–Stokes-α equation is likely to be problematic. The second result evinces the classical scaling relation η/L ~ Re −3/4 for the ratio of the Kolmogorov microscale η to the integral length scale L.   相似文献   

5.
Barenblatt and Botvina with elegant dimensional analysis arguments have elucidated that Paris’ power-law is a weak form of scaling, so that the Paris’ parameters C and m should not be taken as material constants. On the contrary, they are expected to depend on all the dimensionless parameters of the problem, and are really “constants” only within some specific ranges of all these. In the present paper, the dimensional analysis approach by Barenblatt and Botvina is generalized to explore the functional dependencies of m and C on more dimensionless parameters than the original Barenblatt and Botvina, and experimental results are interpreted for a wider range of materials including both metals and concrete. In particular, we find that the size-scale dependencies of m and C and the resulting correlation between C and m are quite different for metals and for quasi-brittle materials, as it is already suggested from the fact the fatigue crack propagation processes lead to m=2-5 in metals and m=10-50 in quasi-brittle materials. Therefore, according to the concepts of complete and incomplete self-similarities, the experimentally observed breakdowns of the classical Paris’ law are discussed and interpreted within a unified theoretical framework. Finally, we show that most attempts to address the deviations from the Paris’ law or the empirical correlations between the constants can be explained with this approach. We also suggest that “incomplete similarity” corresponds to the difficulties encountered so far by the “damage tolerant” approach which, after nearly 50 years since the introduction of Paris’ law, is still not a reliable calculation of damage, as Paris himself admits in a recent review.  相似文献   

6.
When dynamic force is applied to a saturated porous soil, drainage is common. In this paper, the saturated porous soil with a two-phase saturated medium is simulated, and Lamb's integral formulas with drainage and stress formulas for a two-phase saturated medium are given based on Biot's equation and Betti's theorem (the reciprocal theorem). According to the basic solution to Biot's equation, Green's function Gij and three terms of Green's function G4i, Gi4, and G44 of a two-phase saturated medium subject to a concentrated force on a spherical coordinate are presented. The displacement field with drainage, the magnitude of drainage, and the pore pressure of the center explosion source are obtained in computation. The results of the classical Sharpe's solutions and the solutions of the two-phase saturated medium that decays to a single-phase medium are compared. Good agreement is observed.  相似文献   

7.
Utilizing MXene to form the multifunctional derivative is a route to construct high-performance electrode materials. To address this issue, V2CTx MXene was employed to realize a derivative of VOx/V2CTx MXene via a partial oxidation process. Relying on the annealing in the air atmosphere, the controlled oxidation behavior transformed V2CTx MXene partially to V2O5 and formed a derivative hybrid of V2O5/V2CTx MXene. As a result, a package of capacity, rate performance, and cyclability can be enhanced. This work explores the derived behavior of MXenes and provides a route for constructing the hybrid with less interface contact. Furthermore, these findings can be extended to other MXene materials.  相似文献   

8.
An experimental study of evaporation heat transfer coefficients for single circular small tubes was conducted for the flow of C3H8, NH3, and CO2 under various flow conditions. The test matrix encompasses the entire quality range from 0.0 to 1.0, mass fluxes from 50 to 600 kg m−2 s−1, heat fluxes from 5 to 70 kW m−2, and saturation temperatures from 0 to 10 °C. The test section was made of circular stainless steel tubes with inner diameters of 1.5 mm and 3.0 mm, and a length of 2000 mm in a horizontal orientation. The test section was uniformly heated by applying electric power directly to the tubes. The effects of mass flux, heat flux, saturation temperature, and inner tube diameter on the heat transfer coefficient are reported. Among the working refrigerants considered in this study, CO2 has the highest heat transfer coefficient. Laminar flow was observed in the evaporative small tubes, and was considered in the modification of boiling heat transfer coefficients and pressure drop correlations.  相似文献   

9.
Among many presumed-shape pdf approaches for modeling non-premixed turbulent combustion, the presumed β-function pdf is widely used in the literature. However, numerical integration of the β-function pdf may encounter singularity difficulties at mixture fraction values of Z = 0 or 1. To date, this issue has been addressed by few publications. The present study proposes the Piecewise Integration Method (PIM), an efficient, robust and accurate algorithm to overcome these numerical difficulties with the added benefit of improving computational efficiency. Comparison of this method to the existing numerical integration methods shows that the PIM exhibits better accuracy and greatly increases computational efficiency. The PIM treatment of the β-function pdf integration is first applied to the Burke–Schumann solution in conjunction with the k − ε turbulence model to simulate a CH4/H2 bluff-body turbulent flame. The proposed new method is then applied to the same flow using a more complex combustion model, the laminar flamelet model. Numerical predictions obtained by using the proposed β-function pdf integration method are compared to experimental values of the velocity field, temperature and species mass fractions to illustrate the efficiency and accuracy of the present method.  相似文献   

10.
Multiscale material intends to enhance the strength and life of mechanical systems by matching the transmitted spatiotemporal energy distribution to the constituents at the different scale, say—macro, micro, nano, and pico,—, depending on the needs. Lower scale entities are, particularly, critical to small size systems. Large structures are less sensitive to microscopic effects. Scale shifting laws will be developed for relating test data from nano-, micro-, and macro-specimens. The benefit of reinforcement at the lower scale constituents needs to be justified at the macroscopic scale. Filling the void and space in regions of high energy density is considered.Material inhomogeneity interacts with specimen size. Their combined effect is non-equilibrium. Energy exchange between the environment and specimen becomes increasingly more significant as the specimen size is reduced. Perturbation of the operational conditions can further aggravate the situation. Scale transitional functions and/or fj/j+1 are introduced to quantify these characteristics. They are represented, respectively, by , and (fmi/ma,fna/mi,fpi/na). The abbreviations pi, na, mi, and ma refer to pico, nano, micro and macro.Local damage is assumed to initiate at a small scale, grows to a larger scale, and terminate at an even larger scale. The mechanism of energy absorption and dissipation will be introduced to develop a consistent book keeping system. Compaction of mass density for constituents of size 10−12, 10−9, 10−6, 10−3 m, will be considered. Energy dissipation at all scales must be accounted for. Dissipations at the smaller scale must not only be included but they must abide by the same physical and mathematical interpretation, in order to avoid inconsistencies when making connections with those at the larger scale where dissipations are eminent.Three fundamental Problems I, II, and III are stated. They correspond to the commonly used service conditions. Reference is made to a Representative Tip (RT), the location where energy absorption and dissipation takes place. The RT can be a crack tip or a particle. At the larger size scales, RT can refer to a region. Scale shifting of results from the very small to the very large is needed to identify the benefit of using multiscale materials.  相似文献   

11.
The mechanical model was established for the anti-plane fracture problem of a functionally graded coating–substrate system with a coating crack inclined to the weak/micro-discontinuous interface. The Cauchy singular integral equation for the crack was derived using Fourier integral transform, and the Lobatto–Chebyshev collocation method put up by Erdogan and Gupta was used to get its numerical solution. Finally, the effects of the weak/micro-discontinuity of the interface on SIFs were analyzed, the “affected regions” corresponding to the two crack tips have been obtained and their engineering significance was discussed. It was indicated that, for the crack tip in the corresponding “affected region”, to reduce the weak-discontinuity of the interface and to make the interface micro-discontinuous are the two effective ways to reduce the SIF, and the latter way always has more remarkable SIF-reduction effect. For the crack tip outside the “affected region”, its SIF is mainly influenced by material stiffness, and to prevent such a tip from growing toward the interface “softer coating and stiffer substrate” is a more advantageous combination than “stiffer coating and softer substrate”.  相似文献   

12.
We perform three-dimensional dislocation dynamics simulations of solid and annular pillars, having both free-surface boundary conditions, or strong barriers at the outer and/or inner surfaces. Both pillar geometries are observed to exhibit a size effect where smaller pillars are stronger. The scaling observed is consistent with the weakest-link activation mechanism and depends on the solid pillar diameter, or the annular pillar effective diameter, Deff = D − Di, where D and Di are the external and internal diameters of the pillar, respectively. An external strong barrier is observed to dramatically increase the dislocation density by an order of magnitude due to trapping dislocations at the surface. In addition, a considerable increase in the flow strength, by up to 60%, is observed compared to simulations having free-surface boundary conditions. As the applied load increases, weak spots form on the surface of the pillar by dislocations breaking through the surface when the RSS is greater than the barrier strength. The hardening rate is also observed to increase with increasing barrier strength. With cross-slip, we observe dislocations moving to other glide planes, and sometimes double-cross-slipping, producing a thickening of the slip traces at the surface. Finally the results are in qualitative agreement with recent compression experimental results of coated and centrally-filled micropillars.  相似文献   

13.
Finite element simulations are used to study strain localization during uniaxial tensile straining of a single crystal with properties representative of pure Al. The crystal is modeled using a constitutive equation incorporating self- and latent-hardening. The simulations are used to investigate the influence of the initial orientation of the loading axis relative to the crystal, as well as the hardening and strain rate sensitivity of the crystal on the strain to localization. We find that (i) the specimen fails by diffuse necking for strain rate exponents m < 100, and a sharp neck for m > 100. (ii) The strain to localization is a decreasing function of m for m < 100, and is relatively insensitive to m for m > 100. (iii) The strain to localization is a minimum when the tensile axis is close to (but not exactly parallel to) a high symmetry direction such as [1 0 0] or [1 1 1] and the variation of the strain to localization with orientation is highly sensitive to the strain rate exponent and latent-hardening behavior of the crystal. This behavior can be explained in terms of changes in the active slip systems as the initial orientation of the crystal is varied.  相似文献   

14.
The effects of exit-wedge angle on evolution, formation, pinch-off, propagation and diffusive mass entrainment of vortex rings in air were studied using digital particle image velocimetry. Vortex rings were generated by passing a solenoid-valve-controlled air jet through a cylindrical nozzle. Experiments were performed over a wide range of exit-wedge angles (10° ≤ α ≤ 90°) of the cylindrical nozzle, initial Reynolds numbers (450 ≤ Re ≤ 4,580) and length-to-diameter ratios (0.9 ≤ L/D ≤ 11) of the air jet. For sharp edges (α ≤ 10°), a secondary ring may emerge at high Reynolds numbers, which tended to distort the vortex ring if ingested into it. For blunt edges (α ≥ 45°), by contrast, stable vortex rings were produced. The formation phase of a vortex ring was found to be closely related to its evolution pattern. An exit-wedge angle of 45° was found to be optimal for rapid pinch-off and faster propagation and better stability of a vortex ring. Diffusive mass entrainment was found to be between 35% and 40% in the early stages of a vortex ring propagation and it gradually increased throughout the course of vortex ring propagation. Entrainment fraction was found to be sensitive to the L/D ratio of the initial jet and decreases when the L/D ratio is increased.  相似文献   

15.
We have measured time-integrated and time-gated electron temperature (Te) and density (Ne) spatial profiles from indirect-drive implosions. In our experiments, we used a multiple-pinhole two-dimensional imaging spectrometer to obtain multispectral X-ray images of the imploded core. Quantitative comparisons between quasi-monochromatic images in different energy bands allowed Te and Ne spatial profiles to be determined using two independent and validated techniques: a multi-objective search and reconstruction analysis, and an analytical analysis. We then compared the results to a simple one-dimensional (1D) mix-free hydrodynamics simulation in order to evaluate the ability of such a model to predict our experiments. Our data show spatial Te profiles that are qualitatively consistent with the predictions of our 1D simulations, but we observe central cores that are 10–25% cooler and emit X-rays as late as 200 ps after peak compression. We infer time-gated spatial Ne profiles that are consistent with our 1D simulations near the times of peak compression, but we find significant disagreement between time-integrated data and 1D simulation predictions at large radii. Careful analysis of the time-gated and time-integrated Te and Ne spatial profiles, together with streaked X-ray emission spectra from core and shell dopants, suggests mixing of shell material into the core is an important process that our 1D hydrodynamics simulations fail to capture, and comparison between image data and a simple analytical model suggests that 2–5 μm of the initial inner shell thickness mixes into the core during the time period of significant X-ray emission. This mix width is consistent with the predictions of a growth-factor analysis that treats instability growth seeded by capsule surface roughness, and points to the need to consider time-dependent mixing effects when interpreting Te and Ne spatial profiles derived from multispectral X-ray image data, particularly at large radii where mixing effects will be most significant.  相似文献   

16.
Injection of fluids into deep saline aquifers is practiced in several industrial activities, and is being considered as part of a possible mitigation strategy to reduce anthropogenic emissions of carbon dioxide into the atmosphere. Injection of CO2 into deep saline aquifers involves CO2 as a supercritical fluid that is less dense and less viscous than the resident formation water. These fluid properties lead to gravity override and possible viscous fingering. With relatively mild assumptions regarding fluid properties and displacement patterns, an analytical solution may be derived to describe the space–time evolution of the CO2 plume. The solution uses arguments of energy minimization, and reduces to a simple radial form of the Buckley–Leverett solution for conditions of viscous domination. In order to test the applicability of the analytical solution to the CO2 injection problem, we consider a wide range of subsurface conditions, characteristic of sedimentary basins around the world, that are expected to apply to possible CO2 injection scenarios. For comparison, we run numerical simulations with an industry standard simulator, and show that the new analytical solution matches a full numerical solution for the entire range of CO2 injection scenarios considered. The analytical solution provides a tool to estimate practical quantities associated with CO2 injection, including maximum spatial extent of a plume and the shape of the overriding less-dense CO2 front.  相似文献   

17.
Vortex shedding in the wake of two-dimensional bluff bodies is usually accompanied by three dimensional instabilities. These instabilities result in streamwise and vertical vorticity components which occur at a certain spanwise wavelength. The spanwise wavelength of the instabilities (λZ) depends on several parameters, including profile geometry and Reynolds number. The objective of the present work is to study the three dimensional wake instabilities for a blunt trailing edge profiled body, comprised of an elliptical leading edge and a rectangular trailing edge, and to manipulate these instabilities to control the aerodynamic forces. Results of numerical simulations of flow around the body at Re(d) = 400, 600, and 1000, as well as planar Laser Induced Fluorescence (LIF) flow visualizations at Re(d) = 600 and 1000 are analyzed to determine the wake vorticity structure and λZ. Based on the findings of these analyses, an active flow control mechanism for attenuation of the fluctuating aerodynamic forces on the body is proposed. The flow control mechanism is comprised of a series of trailing edge injection ports distributed across the span, with a spacing equal to λZ. Injection of a secondary flow leads to amplification of the three dimensional instabilities and disorganization of the von Kármán vortex street. Numerical simulations indicate that the flow control mechanism can attenuate the fluctuating aerodynamic forces at lower Reynolds numbers (Re(d) = 400 and 600) where λZ is constant in time. However, the control mechanism loses its effectiveness at Re(d) = 1000, due to the temporal variations of λZ.  相似文献   

18.
Detached-eddy simulation (DES) is well understood in thin boundary layers, with the turbulence model in its Reynolds-averaged Navier–Stokes (RANS) mode and flattened grid cells, and in regions of massive separation, with the turbulence model in its large-eddy simulation (LES) mode and grid cells close to isotropic. However its initial formulation, denoted DES97 from here on, can exhibit an incorrect behavior in thick boundary layers and shallow separation regions. This behavior begins when the grid spacing parallel to the wall Δ becomes less than the boundary-layer thickness δ, either through grid refinement or boundary-layer thickening. The grid spacing is then fine enough for the DES length scale to follow the LES branch (and therefore lower the eddy viscosity below the RANS level), but resolved Reynolds stresses deriving from velocity fluctuations (“LES content”) have not replaced the modeled Reynolds stresses. LES content may be lacking because the resolution is not fine enough to fully support it, and/or because of delays in its generation by instabilities. The depleted stresses reduce the skin friction, which can lead to premature separation.For some research studies in small domains, Δ is made much smaller than δ, and LES content is generated intentionally. However for natural DES applications in useful domains, it is preferable to over-ride the DES limiter and maintain RANS behavior in boundary layers, independent of Δ relative to δ. For this purpose, a new version of the technique – referred to as DDES, for Delayed DES – is presented which is based on a simple modification to DES97, similar to one proposed by Menter and Kuntz for the shear–stress transport (SST) model, but applicable to other models. Tests in boundary layers, on a single and a multi-element airfoil, a cylinder, and a backward-facing step demonstrate that RANS function is indeed maintained in thick boundary layers, without preventing LES function after massive separation. The new formulation better fulfills the intent of DES. Two other issues are discussed: the use of DES as a wall model in LES of attached flows, in which the known log-layer mismatch is not resolved by DDES; and a correction that is helpful at low cell Reynolds numbers.  相似文献   

19.
A modified second order viscoelastic constitutive equation is used to derive a kl type turbulence closure to qualitatively assess the effects of elastic stresses on fully-developed channel flow. Specifically, the second order correction to the Newtonian constitutive equation gives rise to a new term in the momentum equation involving the time-averaged elastic shear stress and in the turbulent kinetic energy transport equation quantifying the interaction between the fluctuating elastic stress and rate of strain tensors, denoted by P w , for which a closure is developed and tested. This closure is based on arguments of isotropic turbulence and equilibrium in boundary layer flows and a priori P w could be either positive or negative. When P w is positive, it acts to reduce the production of turbulent kinetic energy and the turbulence model predictions qualitatively agree with direct numerical simulation (DNS) results obtained for more realistic viscoelastic fluid models with memory which exhibit drag reduction. In contrast, P w  < 0 leads to a drag increase and numerical breakdown of the model occurs at very low values of the Deborah number, which signifies the ratio of elastic to viscous stresses. Limitations of the turbulence model primarily stem from the inadequacy of the kl formulation rather than from the closure for P w . An alternative closure for P w , mimicking the viscoelastic stress work predicted by DNS using the Finitely Extensible Nonlinear Elastic-Peterlin fluid model, which is mostly characterized by P w  > 0 but has also a small region of negative P w in the buffer layer, was also successfully tested. This second model for P w leads to predictions of drag reduction, in spite of the enhancement of turbulence production very close to the wall, but the equilibrium conditions in the inertial sub-layer were not strictly maintained.  相似文献   

20.
We present particle tracking velocimetry measurements and flow visualization of pulsatile flow fields in a stented cerebrovascular lateral aneurysm model with a wide ostium anchored on a curved parent vessel. Among the stent parameters, the blocking ratioC α ranging from 0% to 75% was selected to study its effect on the changes of intra-aneurysmal hemodynamics for the reference of minimally invasive endovascular aneurysm treatment. The Womersley number was 3.9 and the mean, peak, and minimal Reynolds numbers based on the bulk average velocity and diameter of the parent vessel were 600, 850, and 300, respectively. The results are characterized in terms of velocity vector field, coded streak images, region averaged velocity, vorticity, and wall shear stress. A critical range ofC α related to the inflow location as well as the shape and number of intra-aneurysmal vortices is identified. The intra-aneurysmal flow activity, vortex strength, and wall shear stress are found to decrease with increasingC α. Among theC α examined,C α=75% is the most favorable in attenuating the risk of aneurysmal rupture and promoting intra-aneurysmal thrombus.  相似文献   

设为首页 | 免责声明 | 关于勤云 | 加入收藏

Copyright©北京勤云科技发展有限公司  京ICP备09084417号