首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 484 毫秒
1.
A double-porosity model is used to describe the oscillatory gas motion and associated contaminant transport induced by cyclical variations in the barometric pressure at the surface of a fractured porous medium. Flow along the fractures and within the permeable matrix blocks is locally one-dimensional. The interaction between fractures and blocks includes seepage of fluid as well as diffusion of contaminant. To guard against artificial numerical diffusion, the FRAM filtering remedy and methodology of Chapman is used in calculating the advective fluxes along fractures and within blocks. The entire system of equations, including the fracture-matrix interaction terms, is solved by a largely implicit non-iterative algorithm which remains stable and conservative even when the computational time step is large compared to the cross-block transit time of pressure waves. The numerical accuracy is tested by comparison with exact solutions for oscillatory and unidirectional flows, some of which include diffusion interaction between the fracture and the matrix. The method is used to estimate the rate of vertical transport of radioactive gases through the rubblized chimney produced by an underground nuclear explosion.  相似文献   

2.
In this work, we propose a novel hydraulic solver in order to simulate key mechanisms that control fluid-driven cracks in the framework of the combined finite-discrete element method (FDEM). The main innovative aspect of the present work is the independence of the fluid's critical time step size with the fracture opening. This advantage is extremely important because it means that very fine meshes can be used around areas of interest, such as boreholes, without penalizing the computational cost as fractures propagate (ie, open) and the fluid flows through them. This is a great advantage over other recently introduced approaches that exhibit a dependency of the time step in the form of Δtcrit ∝ (l/a)2 where l is the element size and a is the fracture aperture. This paper presents a series of benchmark cases for the proposed solver. The rationale adopted by the authors was to benchmark and validate the implementation of the hydraulic solver in an incremental fashion, starting from the simplest cases and building in complexity. The results shown in this work clearly demonstrate that the proposed approach is able to reproduce analytical results for fluid flow through a single crack. The results presented in this paper also demonstrate that the new approach is robust enough to deal with complex fracture patterns and complex geometries; the obtained fluid-driven fracture patterns in the vicinity of a borehole certainly stand to the scrutiny of human visual perception.  相似文献   

3.
Abstract Rock fracture enhances permeability and provides pathways through which fluids migrate. During contact metamorphism, fluids contained in isolated pores and fractures expand in response to temperature increases caused by the dissipation of heat from magmas. Heat transport calculations and thermomechanical properties of water-rich fluids demonstrate (1) that thermal energy is a viable mechanism to produce and maintain pore fluid pressure (Pf) in a contact metamorphic aureole; (2) that the magnitude of Pf generated is sufficient to propagate fractures during the prograde thermal history (cause hydrofracture) and enhance permeability; and (3) that Pf-driven fracture propagation is episodic with time-scales ranging from years to thousands of years. Because Pf dissipation is orders of magnitude faster than P, f buildup, Pf oscillations and cyclical behaviour are generated as thermal heating continues. The Pf cycle amplitude depends on the initial fracture length, geometry and the rock's resistance to failure whereas the frequency of fracture depends on the rate of heating. Consequently, oscillation frequency also varies spatially with distance from the heat source. Time series of fluid pressures caused by this process suggest that cyclical fracture events are restricted to an early time period of the prograde thermal event near the intrusive contact. In the far field, however, individual fracture events have a lower frequency but continue to occur over a longer time interval. Numerous fracture cycles are possible within a single thermal event. This provides a provisional explanation for multiple generations of veins observed in outcrop. P f cycling and oscillations may explain several petrological features. If pore fluids are trapped at various positions along a pressure cycle, the large amplitude of Pf variations for small fractures may account for different pressures recorded by fluid inclusions analysed from a single sample. Pf oscillations, during a single thermal episode, also drive chemical reactions which can produce complex mineral textures and assemblages for discontinuous reactions and/or zoning patterns for continuous reactions. These can mimic polymetamorphic or disequilibrium features. Temporal aspects of fracture propagation and permeability enhancement also constrain the likely timing of fluid flow and fluid-mineral interactions. These data suggest that fluid flow and fluid-mineral reactions are likely to be restricted to an early period in the prograde thermal history, characterized by high Pf coincident with relatively high temperatures, fracture propagation and consequent increases in permeability. This early prograde hydration event is followed by diffusional peak metamorphic reactions. This relationship is evident in the complex mineralogical textures common in some metamorphosed rocks.  相似文献   

4.
A numerical investigation was made of the relationships between fracture initiation, growth, stress field and boundary conditions. Two-dimensional plane strain continuum models were used in which fractures appeared as zones of strain localization developed through application of a strain softening Mohr–Coulomb constitutive model. R and R′ fractures developed first, followed by Y fractures at larger strains. The models showed that equal development of conjugate R and R′ fractures is easily changed to favor one or the other set by minor variations in model initial conditions. Strength loss in fractures caused stress field rotations in regions bounded by fractures, altering the orientation of subsequent fractures. The amount and sense of stress field rotation is dependent on the strength loss during displacement on the fractures, the orientation of fractures, and on the boundary conditions. Y oriented fractures could be explained on the basis of a Mohr–Coulomb failure criterion provided that stress field rotation is accounted for. Monitoring of fracture slip activity showed that, under conditions of constant boundary velocity, slip was discontinuous in time, alternating on fractures throughout the model.  相似文献   

5.
This study was based on the discrete fracture model to investigate the influence of fracture parameters on the solute transport in the fractured rocks of andesite in Lan-Yu island, Taiwan. In the simulation cases, the centers of fractures, fracture lengths and apertures were assumed to have Poisson’s distribution, negative exponential distribution and lognormal distribution, respectively. With the above assumptions, constructing the discrete fracture model became practicable. Using the mass-balance equation with specified boundary conditions, the flow field in the rock was solved. Then particles were released under the flow field. Monte Carlo method was used assuming that the amount of particles was proportional to the flow rates to get the particle accumulated percentage breakthrough curve and to estimate the dispersion coefficient. On the basis of the discrete fracture model, it was possible to evaluate the property of dispersion behavior of andesite in Lan-Yu Island with flow and transport mechanism. Properties of the dispersion behavior such as the relation between distance and traveling-time (ln〈r 2〉 and ln 〈t〉), anisotropic behavior, and the overall dispersion coefficient in a fracture network were characterized: the slope value of ln〈r 2〉 and ln〈t〉 was 1.64 an indication of non-Fickian dispersion, the particles dispersion along the flow (D11) was bigger than that perpendicular to the flow (D22), and the dispersion coefficient by this study was 0.91 m comparing the value 1 m from Sauty’s method.  相似文献   

6.
Analysis of contaminant transport through fractured crystalline rocks has received considerable attention, particularly with regard to subsurface nuclear waste repositories. Most of the studies have employed the dual continuum approach, with the fractures and the rock matrix as the two continuums, assuming that fractures control the overall conductivity of the rock and the porous matrix just provides storage. However, field observations of rock fractures have shown that the real situation can be very complex. Based on some recent investigations, it has been reported that the portion of the rock matrix adjacent to many open fractures is physically and chemically altered. These alterations, referred to as the fracture skin, can have different sorption and diffusion properties compared to those of the undisturbed rock matrix and this may influence the transport of solutes through such formations. In the present study, a numerical model is developed to simulate conservative solute transport in a fractured crystalline rock formation using the triple continuum approach ?? with the fracture, fracture skin and the rock matrix as the three continuums. The model is solved using a fully implicit finite difference scheme. Contaminant migration in the fractured formation with and without skin has been simulated. It is observed that contaminant penetration along the fracture is enhanced at large flow velocities. The effect of flow velocity on conservative solute transport is investigated for different fracture apertures and fracture skin thicknesses. The influence of flow velocity on contaminant transport is demonstrated to be more with change in fracture aperture than with change in skin thickness.  相似文献   

7.
Spreading of a non-aqueous phase liquid (NAPL) denser than water (DNAPL) lens (mound) in the unsaturated zone of double-porosity aquifer above an impervious plane boundary is investigated. The double-porosity aquifer is conceptualized as a fracture network surrounding pervious blocks. Vertical gravity equilibrium is assumed to prevail in each one of the two media, fractures and blocks. Through vertical integration, two coupled partial differential equations for the DNAPL content in each medium are obtained. The mass exchange rate between high- and low-permeability media is considered as a function of NAPL content. The dominant effect is gravity, whereas capillary forces are negligible. Analytical solutions for one-dimensional and axisymmetric problems of mound spreading are obtained.  相似文献   

8.
Analysis of fracture network connectivity using percolation theory   总被引:2,自引:0,他引:2  
Connectivity aspects of fracture networks are analyzed in terms of percolation theory. These aspects are of fundamental importance in characterization, exploitation, and management of fractured formations. General connectivity and power law relationships are determined that characterize the density of fractures and average number of intersections per fracture necessary to ensure network connectivity, the likelihood of a fractured formation being hydraulically connected, and the probability that any specific fracture is connected to the conducting portion of the network. Monte Carlo experiments with a two-dimensional fracture network model confirm the percolation theory predictions. These relationships may prove useful in formulating theoretically tractable approximations of fracture nerworks that capture the essential system properties.  相似文献   

9.
Hou  Michael Z.  Li  Mengting  Gou  Yang  Feng  Wentao 《Acta Geotechnica》2021,16(4):985-1000

Hydraulic fracturing is an essential technology for the development of unconventional resources such as tight gas. The evaluation of the fracture performance and productivity is important for the design of fracturing operations. However, the traditional dimensionless fracture conductivity is too simple to be applied in real fracturing operations. In this work, we proposed a new model of dimensionless fracture conductivity (FCD), which considers the irregular fracture geometry, proppant position and concentration. It was based on the numerical study of the multistage hydraulic fracturing and production in a tight gas horizontal well of the North German Basin. A self-developed full 3D hydraulic fracturing model, FLAC3Dplus, was combined with a sensitive/reliability analysis and robust design optimization tool optiSLang and reservoir simulator TMVOCMP to achieve an automatic history matching as well as simulation of the gas production. With this tool chain, the four fracturing stages were history matched. The simulation results show that all four fractures have different geometry and proppant distribution, which is mainly due to different stress states and injection schedule. The position and concentration of the proppant play important roles for the later production, which is not considered in the traditional dimensionless fracture conductivity FCD. In comparison, the newly proposed formulation of FCD could predict the productivity more accurately and is better for the posttreatment evaluation.

  相似文献   

10.
Summary This paper presents a theoretical study on normally incident elastic P-wave transmission across single dry fractures with a nonlinear normal deformational behavior. The effects of nonlinear fracture normal behavior on P-wave transmission are examined without the mixture of fracture shear behavior. The linear displacement discontinuity model for wave propagation across fractures is extended to a nonlinear model – the hyperbolic elastic model (BB model). Numeric solutions of magnitudes of transmission (|T non|) and reflection (|R non|) coefficients, for normally incident P-wave transmission across the nonlinear deformable fractures, are obtained and related to the closure behavior of fractures. Parametric studies are conducted to acquire an insight into the effects of the nonlinear fracture normal deformation on P-wave transmission, in terms of initial normal stiffness and the ratio of current maximum closure to maximum allowable closure of the fractures, as well as the incident wave amplitude and frequency. Comparisons between the linear and nonlinear models are presented. It is shown that, |T lin| and |R lin| for the linear model are special solutions of |T non| and |R non| for the nonlinear model, when the incident wave amplitude is so low that the current maximum closure of fracture incurred during the wave transmission is much smaller, relative to the maximum allowable closure. In addition, the nonlinear fracture behavior gives rise to a phenomenon of higher harmonics during the wave transmission across the fracture. The higher harmonics contribute to the increase of |T non| from |T lin|.  相似文献   

11.
Detailed mapping and biostratigraphic data provide new evidence of a major break below the base of the Caradoc succession along the northern margin of the Harlech Dome in south central Snowdonia. Within this outcrop the sequence is locally complicated by subsequent tectonic and volcanic events, but undisturbed sections indicate a break of at least 10 Ma between upper Arenig and middle, or upper, Llandeilo strata. The break is greatest between two N-S trending fracture systems, the Cwm Pennant Fault Zone in the west and Trawsfynydd Fault Zone in the east, which both have a persistent history of reactivation. Between these two fractures neither Llanvirn nor lower Llandeilo strata occur. This contrasts with the thick, and perhaps complete, sequence preserved in the Cadair Idris district on the southern margin of the dome and suggests that, during Llanvirn times, the Harlech Dome formed a major uplifted and tilted block, with a tectonically active northern margin. Subsequently, uplift and tectonism either ceased before, or was overwhelmed by, the ensuing sea-level rise associated with the gracilis (early Caradoc) transgression. The preservation of ooidal ironstones around the dome suggests that it may have became a large shoal or platformal area at this time. Renewed uplift and erosion along the northern margin of the Harlech Dome during the early Caradoc (gracilis to multidens) led to large-scale disruption of the stratigraphic succession by mass gravity flow and slumping, overprinting and locally accentuating the effects of the earlier hiatus.  相似文献   

12.
Empiricaldecline curve forecasts of future production rates and cumulative production to be expected from certain classes of petroleum reservoirs has been a methodology practiced and relied on by petroleum geologists for many years. New ways to consider the subject are described in this paper. Thereby it is seen that these forecasts sometimes can be made even when some of the record of the past performance history is sparse (say, because it has been lost, or never recorded).SI Units D, D i initial reference production rate,s –1 - n dimensionless constant [cf. Eq. (1)] - N cumulative production, m3 - Q, Q i production rate att and initially, m3/s - R dimensionless time variable [cf. Eq. (9)] - t time, s  相似文献   

13.
Summary.  The purpose of this paper is to establish a dynamic constitutive model of fracture normal behaviour, based on laboratory tests of artificial fractures cast by cement mortar. A series of tests are systematically carried out under quasi-static (10−1 MPa/s) up to highly dynamic (103 MPa/s) monotonic loading conditions. The normal stress-fracture closure response is measured at different loading rates. Based on the measured curves, a nonlinear (hyperbolic) dynamic model of fracture normal behaviour, termed as dynamic BB model, is proposed. The dynamic model is modified from the existing BB model of static normal behaviour of fractures by taking into account the loading-rate effect. Two important dynamic parameters of fractures, FSC d (dynamic fracture stiffness constant, which describes the incremental ratio of dynamic initial stiffness) and FCC d (dynamic fracture closure constant, which describes the decremental ratio of dynamic maximum allowable closure), are identified. They indicate the quantitative degree of loading-rate effect on fracture normal behaviour subjected to dynamic loads. For practical application, the new model is incorporated into the Universal Distinct Element Code (UDEC) and subsequently, UDEC modelling of normally incident P-wave transmission across single fractures with the dynamic BB model is conducted. Wave transmission coefficient is obtained for various combinations of fracture dynamic parameters, as well as different wave amplitudes and frequencies. The numerical results show that wave transmission coefficient for a fracture with the dynamic BB model is greater than that for a fracture with the static BB model. In addition, a fracture with higher values of FSC d and FCC d leads to higher transmission (lower attenuation). Author’s address: J. Zhao, Ecole Polytechnique Federale de Lausanne (EPFL), Rock Mechanics Laboratory, 1015 Lausanne, Switzerland  相似文献   

14.
The response of deformable fractures to changes in fluid pressure controls phenomena ranging from the flow of fluids near wells to the propagation of hydraulic fractures. We developed an analysis designed to simulate fluid flows in the vicinity of asperity‐supported fractures at rest, or fully open fractures that might be propagating. Transitions between at‐rest and propagating fractures can also be simulated. This is accomplished by defining contact aperture as the aperture when asperities on a closing fracture first make contact. Locations on a fracture where the aperture is less than the contact aperture are loaded by both fluid pressure and effective stress, whereas locations where the aperture exceeds the contact aperture are loaded only by fluid pressure. Fluid pressure and effective stress on the fracture are determined as functions of time by solving equations of continuity in the fracture and matrix, and by matching the global displacements of the fracture walls to the local deformation of asperities. The resulting analysis is implemented in a numerical code that can simulate well tests or hydraulic fracturing operations. Aperture changes during hydraulic well tests can be measured in the field, and the results predicted using this analysis are similar to field observations. The hydraulic fracturing process can be simulated from the inflation of a pre‐existing crack, to the propagation of a fracture, and the closure of the fracture to rest on asperities or proppant. Two‐dimensional, multi‐phase fluid flow in the matrix is included to provide details that are obscured by simplifications of the leakoff process (Carter‐type assumptions) used in many hydraulic fracture models. Execution times are relatively short, so it is practical to implement this code with parameter estimation algorithms to facilitate interpretation of field data. Copyright © 2006 John Wiley & Sons, Ltd.  相似文献   

15.
Natural fractures are characterized by rough surfaces and complex fluid flows. A large distribution of apertures (residual voids) within their walls and the presence of contact points (in situ normal loads) produce heterogeneous flows (channeling). The resulting permeabilities, porosities or fluid–rock exchange surfaces cannot be realistically modeled by parallel and smooth plate models. Four natural fractures are sampled at different depths and degrees of alteration in the Soultz sandstone and granite (EPS1 drillhole, Soultz-sous-Forêts, Bas-Rhin, France). The fracture surfaces are measured with mechanical profilometry and maps of asperity heights (XYZ). Resulting local apertures (XYe) are then calculated. A statistical study of the surface profiles (XZ) show that the fractures are more or less rough and tortuous according to the types of alteration. Altered samples are characterized by smoother surfaces of fractures. Such differences imply that (i) the average fracture aperture is not representative for the whole fracture and that (ii) the different local apertures should be integrated in hydraulic and mechanical models. A hydraulic model (finite difference calculations) of fluid flow, taking into account the elastic closure (Hertz contact theory) of fractures with depth, is used. Maps of contact points and relative local loads within the fracture planes are compared to flow maps. They show different channeling of fluid flows. Strongly altered fractures are characterized by homogeneous fluxes despite the presence of numerous contact zones during the closure of fracture. By contrast, fresh fractures develop, increasing fluid flow channels with depth.Fracture closure (increasing normal stress) does not systematically increase the channeling of fluid flow. There is evidence for a general smoothing out of the irregularities of the fracture walls due to precipitation of secondary minerals, indicating that the cubic law can be commonly valid, also at great crustal depth but this validity depends on the degree of fracture alteration. Mineralogical and geochemical observations, thus, should be taken into account to perform more accurate permeability calculations and models of fluid circulation in fracture networks.  相似文献   

16.
The geometry and thermal history of fractures have been determined at 59 stations from Reykjavik to Hvalfjördur in southwestern Iceland. The data provide information on crustal stress regimes in the vicinity of mid-ocean ridges.Two major, generalized fracture orientations are present
1. (1) a northeast system, trend 010°–030°, except on Akranes where the orientation is 040°–060°
2. (2) a broad east—west system containing one or more sets with strike between 070°–130°.
Thermal history of the host rock and fractures was determined from secondary minerals in vugs and fractures. The thermal history indicates that the northeast fracture set opened while the area was within the relatively hot axial zone of active volcanism and rifting. Some of the east—west trending fractures also opened at this time but many formed later, after the area had begun to cool and drift from the active zone.The northeast fracture set is essentially parallel to the trend of dikes and normal faults in southwestern Iceland. They have been interpreted as extension fractures (resulting in about 0.4% maximum extension) forming generally from the same stress field associated with normal faulting and dike injection in the active zone. Fracturing in an east-west direction (estimated 0.1% maximum extension), mainly near the edge and outside the active zone, indicates a reorientation of this stress field. The dominant mechanism related to the origin of the east—west fractures may be thermoelastic stresses arising from axial and basal accretion and cooling of lithospheric plates.Both fracture systems are inferred to have formed, in the Griffiths idealization, under nearly biaxial effective compressive loading on the order of 200 bar. The discrepancy between this value and the kilobar-order strengths of short-time laboratory tests reflects such factors as high temperature stress corrosion and fatigue. Fracture propagation is assumed to have been stable, but governed primarily by lateral load-diminishing mechanisms rather than by progressive loading. These relaxation mechanisms may have been episodic (northeast-system fissure swarm activity) or steady-state (thermoelastic contraction) in time.  相似文献   

17.
A two-dimensional boundary integral method to analyse the flow of contaminant in fractured media having a two- or three-dimensional orthogonal fracture network is presented. The method assumes that the fractures provide the paths of least resistance for transport of contaminants while the matrix, because of its low permeability, acts as ‘storage blocks’ into which the contaminant diffuses. Laplace transform is used to eliminate the time variable in the governing equation in order to facilitate the formulation of a boundary integral equation in the Laplace transform space. Conventional boundary element techniques are applied to solve for the contaminant concentrations at specified locations in the spatial domain. The concentration in the time domain is then obtained by using an efficient inversion technique developed by Talbot. The method is able to analyse the behaviour of waste repositories which have diminishing concentration due to the mass transport of the contaminant into the surrounding fractured media.  相似文献   

18.
Summary Thirteen natural rock profiles (Barton and Choubey, 1977) are analyzed for their fractal properties. Most of the profiles were found to approximate fractal curves but some also showed features of specific wavelengths and amplitudes superimposed on fractal characteristics. The profiles showed fractal dimensions from 1.1 to 1.5 covering a range of selfsimilar and self-affine curves. The analysis results suggest a negative correlation between fractal dimension,D, and amplitude,A. Joint roughness coefficients (JRC) show a positive correlation with amplitude,A, and a negative correlation with fractal dimension,D. A numerical model of fracture closure is used to investigate the effects of different profile characteristics (D, A and sample size) on the nature of dilation and contact area, using the natural profiles and synthetic fractional Brownian motion profiles. Smooth profiles (low JRC, highD, lowA) display many small contact regions whereas rough fractures (high JRC, lowD, highA) display few large contact areas. The agreement with published experimental data supports the suggested correlations between JRC and the fractal parameters,A andD. It is suggested that observed scale effects in JRC and joint dilation can be explained by small differential strain discontinuities across fractures, which originate at the time of fracture formation.  相似文献   

19.
Occurrence of springs in massifs of crystalline rocks,northern Portugal   总被引:1,自引:0,他引:1  
An inventory of artesian springs emerging from fractures (fracture springs) was conducted in the Pinh?o River Basin and Morais Massif, northern Portugal, comprising an area of approximately 650 km2. Over 1,500 springs were identified and associated with geological domains and fracture sets. Using cross-tabulation analysis, spring distributions by fracture sets were compared among geological environments, and the deviations related to differences in rock structure and, presumably, to differences in deformational histories. The relation between spring frequencies and rock structures was further investigated by spectral determination, the model introduced in this study. Input data are the spring frequencies and fracture lengths in each geological domain, in addition to the angles between fracture strikes and present-day stress-field orientation (θ). The model's output includes the so-called intrinsic densities, a parameter indexing spring occurrence to factors such as fracture type and associated deformational regime and age. The highest densities (12.2 springs/km of lineament) were associated with young shear fractures produced by brittle deformation, and the lowest (0.1) with old tensional and ductile fractures. Spectral determination also relates each orientation class to a dominant structural parameter: where spring occurrence is controlled by θ, the class is parallel to the present-day stress-field orientation; where the control is attributed to the length of fractures, the spring occurrence follows the strike of large-scale normal faults crossing the region. Electronic Publication  相似文献   

20.
In this paper, a fully coupled 2‐dimensional poroelastic displacement discontinuity method is used to investigate the refracturing process in horizontal wells. One of the objectives of refracturing is to access new reserves by adding new hydraulic fractures in zones that were bypassed in the initial fracturing attempt. Pore pressure depletion in the vicinity of old fractures directly affects the state of stress and eventually the propagation of newly created hydraulic fractures. Thus, a poroelastic analysis is required to identify guidelines for the refracturing process, in particular to understand the extension of the pore pressure depletion, and eventually, the orientation of new as well as old fractures. We propose a fully coupled approach to model the whole process of child fracture propagation in a depleted area between 2 parent fractures in the same wellbore. This approach omits the need of using multistep workflow that is regularly used to model the process. The maximum tensile stress criterion (σ criterion) is used for hydraulic fracture propagation. The proposed method is verified using available analytical solutions for total stress and pore pressure loading modes on a line fracture in drained and undrained conditions. Then, test cases of multifractured horizontal wells are studied to calculate the time evolution of the stress and pore pressure fields around old fractures and to understand the effect of these fields on the propagation path of newly created fractures. Finally, the effect of the pore pressure depletion on the propagation path of the newly created fractures in the bypassed area of the same wellbore is studied. The results show that the depleted areas around old fractures are highly affected by the extent and severity of the stress redistribution and pore pressure depletion. It is observed that a successful creation of new fractures may only happen in certain time frames. The results of this study provide new insights on the behavior of newly created fractures in depleted zones. They also clarify the relationship between stress change and pore pressure depletion in horizontal wells.  相似文献   

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

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