首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 31 毫秒
1.
CO2 injected into porous formations is accommodated by reduction in the volume of the formation fluid and enlargement of the pore space, through compression of the formation fluids and rock material, respectively. A critical issue is how the resulting pressure buildup will affect the mechanical integrity of the host formation and caprock. Building on an existing approximate solution for formations of infinite radial extent, this article presents an explicit approximate solution for estimating pressure buildup due to injection of CO2 into closed brine aquifers of finite radial extent. The analysis is also applicable for injection into a formation containing multiple wells, in which each well acts as if it were in a quasi-circular closed region. The approximate solution is validated by comparison with vertically averaged results obtained using TOUGH2 with ECO2N (where many of the simplifying assumptions are relaxed), and is shown to be very accurate over wide ranges of the relevant parameter space. The resulting equations for the pressure distribution are explicit, and can be easily implemented within spreadsheet software for estimating CO2 injection capacity.  相似文献   

2.
This article presents a numerical modeling application using the code TOUGHREACT of a leakage scenario occurring during a CO2 geological storage performed in the Jurassic Dogger formation in the Paris Basin. This geological formation has been intensively used for geothermal purposes and is now under consideration as a site for the French national program of reducing greenhouse gas emissions and CO2 geological storage. Albian sandstone, situated above the Dogger limestone is a major strategic potable water aquifer; the impacts of leaking CO2 due to potential integrity failure have, therefore, to be investigated. The present case–study illustrates both the capacity and the limitations of numerical tools to address such a critical issue. The physical and chemical processes simulated in this study have been restricted to: (i) supercritical CO2 injection and storage within the Dogger reservoir aquifer, (ii) CO2 upwards migration through the leakage zone represented as a 1D vertical porous medium to simulate the cement–rock formation interface in the abandoned well, and (iii) impacts on the Albian aquifer water quality in terms of chemical composition and the mineral phases representative of the porous rock by estimating fluid–rock interactions in both aquifers. Because of CPU time and memory constraints, approximation and simplification regarding the geometry of the geological structure, the mineralogical assemblages and the injection period (up to 5 years) have been applied to the system, resulting in limited analysis of the estimated impacts. The CO2 migration rate and the quantity of CO2 arriving as free gas and dissolving, firstly in the storage water and secondly in the water of the overlying aquifer, are calculated. CO2 dissolution into the Dogger aquifer induces a pH drop from about 7.3 to 4.9 limited by calcite dissolution buffering. Glauconite present in the Albian aquifer also dissolves, causing an increase of the silicon and aluminum in solution and triggering the precipitation of kaolinite and quartz around the intrusion point. A sensitivity analysis of the leakage rate according to the location of the leaky well and the variability of the petro-physical properties of the reservoir, the leaky well zone and the Albian aquifers is also provided.  相似文献   

3.
Dissolution of CO2 into brine is an important and favorable trapping mechanism for geologic storage of CO2. There are scenarios, however, where dissolved CO2 may migrate out of the storage reservoir. Under these conditions, CO2 will exsolve from solution during depressurization of the brine, leading to the formation of separate phase CO2. For example, a CO2 sequestration system with a brine-permeable caprock may be favored to allow for pressure relief in the sequestration reservoir. In this case, CO2-rich brine may be transported upwards along a pressure gradient caused by CO2 injection. Here we conduct an experimental study of CO2 exsolution to observe the behavior of exsolved gas under a wide range of depressurization. Exsolution experiments in highly permeable Berea sandstones and low permeability Mount Simon sandstones are presented. Using X-ray CT scanning, the evolution of gas phase CO2 and its spatial distribution is observed. In addition, we measure relative permeability for exsolved CO2 and water in sandstone rocks based on mass balances and continuous observation of the pressure drop across the core from 12.41 to 2.76 MPa. The results show that the minimum CO2 saturation at which the exsolved CO2 phase mobilization occurs is from 11.7 to 15.5%. Exsolved CO2 is distributed uniformly in homogeneous rock samples with no statistical correlation between porosity and CO2 saturation observed. No gravitational redistribution of exsolved CO2 was observed after depressurization, even in the high permeability core. Significant differences exist between the exsolved CO2 and water relative permeabilities, compared to relative permeabilities derived from steady-state drainage relative permeability measurements in the same cores. Specifically, very low CO2 and water relative permeabilities are measured in the exsolution experiments, even when the CO2 saturation is as high as 40%. The large relative permeability reduction in both the water and CO2 phases is hypothesized to result from the presence of disconnected gas bubbles in this two-phase flow system. This feature is also thought to be favorable for storage security after CO2 injection.  相似文献   

4.
We have used the TOUGH2-MP/ECO2N code to perform numerical simulation studies of the long-term behavior of CO2 stored in an aquifer with a sloping caprock. This problem is of great practical interest, and is very challenging due to the importance of multi-scale processes. We find that the mechanism of plume advance is different from what is seen in a forced immiscible displacement, such as gas injection into a water-saturated medium. Instead of pushing the water forward, the plume advances because the vertical pressure gradients within the plume are smaller than hydrostatic, causing the groundwater column to collapse ahead of the plume tip. Increased resistance to vertical flow of aqueous phase in anisotropic media leads to reduced speed of up-dip plume advancement. Vertical equilibrium models that ignore effects of vertical flow will overpredict the speed of plume advancement. The CO2 plume becomes thinner as it advances, but the speed of advancement remains constant over the entire simulation period of up to 400 years, with migration distances of more than 80 km. Our simulations include dissolution of CO2 into the aqueous phase and associated density increase, and molecular diffusion. However, no convection develops in the aqueous phase because it is suppressed by the relatively coarse (sub-) horizontal gridding required in a regional-scale model. A first crude sub-grid-scale model was developed to represent convective enhancement of CO2 dissolution. This process is found to greatly reduce the thickness of the CO2 plume, but, for the parameters used in our simulations, does not affect the speed of plume advancement.  相似文献   

5.
The reinjection of sour or acid gas mixtures is often required for the exploitation of hydrocarbon reservoirs containing remarkable amounts of acid gases (H2S and CO2) to reduce the environmental impact of field exploitation and provide pressure support for enhanced oil recovery (EOR) purposes. Sour and acid gas injection in geological structures can be modelled with TMGAS, a new Equation of State (EOS) module for the TOUGH2 reservoir simulator. TMGAS can simulate the two-phase behaviour of NaCl-dominated brines in equilibrium with a non-aqueous (NA) phase, made up of inorganic gases such as CO2 and H2S and hydrocarbons (pure as well as pseudo-components), up to the high pressures (~100 MPa) and temperatures (~200°C) found in deep sedimentary basins. This study is focused on the near-wellbore processes driven by the injection of an acid gas mixture in a hypothetical high-pressure, under-saturated sour oil reservoir at a well-sector scale and at conditions for which the injected gas is fully miscible with the oil. Relevant-coupled processes are simulated, including the displacement of oil originally in place, the evaporation of connate brine, the salt concentration and consequent halite precipitation, as well as non-isothermal effects generated by the injection of the acid gas mixture at temperatures lower than initial reservoir temperature. Non-isothermal effects are studied by modelling in a coupled way wellbore and reservoir flow with a modified version of the TOUGH2 reservoir simulator. The described approach is limited to single-phase wellbore flow conditions occurring when injecting sour, acid or greenhouse gas mixtures in high-pressure geological structures.  相似文献   

6.
Incompressible 3-D DNS is performed in non-decaying turbulence with single step chemistry to validate a new analytical expression for turbulent burning velocity. The proposed expression is given as a sum of laminar and turbulent contributions, the latter of which is given as a product of turbulent diffusivity in unburned gas and inverse scale of wrinkling at the leading edge. The bending behavior of U T at higher u′ was successfully reproduced by the proposed expression. It is due to decrease in the inverse scale of wrinkling at the leading edge, which is related with an asymmetric profile of FSD with increasing u′. Good agreement is achieved between the analytical expression and the turbulent burning velocities from DNS throughout the wrinkled, corrugated and thin reaction zone regimes. Results show consistent behavior with most experimental correlations in literature including those by Bradley et al. (Philos Trans R Soc Lond A 338:359–387, 1992), Peters (J Fluid Mech 384:107–132, 1999) and Lipatnikov et al. (Progr Energ Combust Sci 28:1–74, 2002).  相似文献   

7.
Kazemi et al. (SPE Reserv Eng 7(2):219–227, 1992) suggested an empirical matrix-fracture transfer function, verified based on experimental data of Mattax and Kyte (Trans AIME 225(15):177–184, 1962), to model fluid flow in naturally fractured dual porosity petroleum reservoirs using a dual-porosity numerical simulator. Their generalized shape factor should be valid for all possible irregular matrix blocks. The factor is calculated based on the volume of the matrix block, the surface open to flow in all directions and the distances of these surfaces to the centre of the matrix block. The summation is done over all open surfaces of a matrix block. Kazemi et al. (1992) showed that for rectangles and cylinders the formula reduces to the well-known forms of the shape factor. By the time, many authors indicated the validity of the formula, but no theoretical proof was offered for that so far. This study derives the Kazemi et al. (1992) shape factor using control volume finite difference discretization on the fracture-matrix dual continuum. The matrix blocks are handled as Voronoi polyhedra. The derivation is given for both isotropic and tensorial matrix permeability. Based on this derivation the authors conclude that the Kazemi et al. (SPE Reserv Eng 7(2):219–227, 1992) formula is exact under pseudo-steady-state conditions within the dual continuum mathematical concept of natural fractured dual porosity systems.  相似文献   

8.
Recently, the tube diameter relaxation time in the evolution equation of the molecular stress function (MSF) model (Wagner et al., J Rheol 49: 1317–1327, 2005) with the interchain pressure effect (Marrucci and Ianniruberto, Macromolecules 37:3934–3942, 2004) included was shown to be equal to three times the Rouse time in the limit of small chain stretch. From this result, an advanced version of the MSF model was proposed, allowing modeling of the transient and steady-state elongational viscosity data of monodisperse polystyrene melts without using any nonlinear parameter, i.e., solely based on the linear viscoelastic characterization of the melts (Wagner and Rolón-Garrido 2009a, b). In this work, the same approach is extended to model experimental data in shear flow. The shear viscosity of two polybutadiene solutions (Ravindranath and Wang, J Rheol 52(3):681–695, 2008), of four styrene-butadiene random copolymer melts (Boukany et al., J Rheol 53(3):617–629, 2009), and of four polyisoprene melts (Auhl et al., J Rheol 52(3):801–835, 2008) as well as the shear viscosity and the first and second normal stress differences of a polystyrene melt (Schweizer et al., J Rheol 48(6):1345–1363, 2004), are analyzed. The capability of the MSF model with the interchain pressure effect included in the evolution equation of the chain stretch to model shear rheology on the basis of linear viscoelastic data alone is confirmed.  相似文献   

9.
The turbulence structure near a wall is a very active subject of research and a key to the understanding and modeling of this flow. Many researchers have worked on this subject since the fifties Hama et al. (J Appl Phys 28:388–394, 1957). One way to study this organization consists of computing the spatial two-point correlations. Stanislas et al. (C R Acad Sci Paris 327(2b):55–61, 1999) and Kahler (Exp Fluids 36:114–130, 2004) showed that double spatial correlations can be computed from stereoscopic particle image velocimetry (SPIV) fields and can lead to a better understanding of the turbulent flow organization. The limitation is that the correlation is only computed in the PIV plane. The idea of the present paper is to propose a new method based on a specific stereoscopic PIV experiment that allows the computation of the full 3D spatial correlation tensor. The results obtained are validated by comparison with 2D computation from SPIV. They are in very good agreement with the results of Ganapthisubramani et al. (J Fluid Mech 524:57–80, 2005a).  相似文献   

10.
This paper is concerned with the decay structure for linear symmetric hyperbolic systems with relaxation. When the relaxation matrix is symmetric, the dissipative structure of the systems is completely characterized by the Kawashima–Shizuta stability condition formulated in Umeda et al. (Jpn J Appl Math 1:435–457, 1984) and Shizuta and Kawashima (Hokkaido Math J 14:249–275, 1985) and we obtain the asymptotic stability result together with the explicit time-decay rate under that stability condition. However, some physical models which satisfy the stability condition have non-symmetric relaxation term (for example, the Timoshenko system and the Euler–Maxwell system). Moreover, it had been already known that the dissipative structure of such systems is weaker than the standard type and is of the regularity-loss type (see Duan in J Hyperbolic Differ Equ 8:375–413, 2011; Ide et al. in Math Models Meth Appl Sci 18:647–667, 2008; Ide and Kawashima in Math Models Meth Appl Sci 18:1001–1025, 2008; Ueda et al. in SIAM J Math Anal 2012; Ueda and Kawashima in Methods Appl Anal 2012). Therefore our purpose in this paper is to formulate a new structural condition which includes the Kawashima–Shizuta condition, and to analyze the weak dissipative structure for general systems with non-symmetric relaxation.  相似文献   

11.
The theory of thin wires developed in Dret and Meunier (Comptes Rendus de l’Académie des Sciences. Série I. Mathématique 337:143–147, 2003) is adapted to phase-transforming materials with large elastic moduli in the sense discussed in James and Rizzoni (J Elast 59:399–436, 2000). The result is a one-dimensional constitutive model for shape memory wires, characterized by a small number of material constants. The model is used to analyze self-accommodated and detwinned microstructures and to study superelasticity. It also turns out that the model successfully reproduces the behavior of shape memory wires in experiments of restrained recovery (Tsoi et al. in Mater Sci Eng A 368:299–310, 2004; Tsoi in 50:3535–3544, 2002; S̆ittner et al. in Mater Sci Eng A 286:298–311, 2000; vokoun in Smart Mater Struct 12:680–685, 2003; Zheng and Cui in Intermetallics 12:1305–1309, 2004; Zheng et al. in J Mater Sci Technol 20(4):390–394, 2004). In particular, the model is able to predict the shift to higher transformation temperatures on heating. The model also captures the effect of prestraining on the evolution of the recovery stress and of the martensite volume fraction.  相似文献   

12.
Carbon storage in saline formations is considered as a promising option to ensure the necessary decrease of CO2 anthropogenic emissions. Its industrial development in those formations is above all conditioned by its safety demonstration. Assessing the evolution of trapped and mobile CO2 across time is essential in the perspective of reducing leakage risks. In this work, we focus on residual trapping phenomenon occurring during the wetting of the injected CO2 plume. History dependent effects are of first importance when dealing with capillary trapping. We then apply the classical fractional flow theory (Buckley–Leverett type model) and include trapping and hysteresis models; we derive an analytical solution for the temporal evolution of saturation profile and of CO2 trapped quantity when injecting water after the gas injection (“artificial imbibition”). The comparison to numerical simulations for different configurations shows satisfactory match and justifies, in the case of industrial CO2 storage, the assumptions of incompressible flow with no consideration of capillary pressure. The obtained analytical solution allows the quick assessment of both the quantity and the location of mobile gas left during imbibition.  相似文献   

13.
This paper is dedicated to the study of viscous compressible barotropic fluids in dimension N ≧ 2. We address the question of the global existence of strong solutions for initial data close to a constant state having critical Besov regularity. First, this article shows the recent results of Charve and Danchin (Arch Ration Mech Anal 198(1):233–271, 2010) and Chen et al. (Commun Pure Appl Math 63:1173–1224, 2010) with a new proof. Our result relies on a new a priori estimate for the velocity that we derive via the intermediary of the effective velocity, which allows us to cancel out the coupling between the density and the velocity as in Haspot (Well-posedness in critical spaces for barotropic viscous fluids, 2009). Second, we improve the results of Charve and Danchin (2010) and Chen et al. (2010) by adding as in Charve and Danchin (2010) some regularity on the initial data in low frequencies. In this case we obtain global strong solutions for a class of large initial data which rely on the results of Hoff (Arch Rational Mech Anal 139:303–354, 1997), Hoff (Commun Pure Appl Math 55(11):1365–1407, 2002), and Hoff (J Math Fluid Mech 7(3):315–338, 2005) and those of Charve and Danchin (2010) and Chen et al. (2010). We conclude by generalizing these results for general viscosity coefficients.  相似文献   

14.
We used the multiphase and multicomponent TOUGH2/EOS7CA model to carry out predictive simulations of CO2 injection into the shallow subsurface of an agricultural field in Bozeman, Montana. The purpose of the simulations was to inform the choice of CO2 injection rate and design of monitoring and detection activities for a CO2 release experiment. The release experiment configuration consists of a long horizontal well (70 m) installed at a depth of approximately 2.5 m into which CO2 is injected to mimic leakage from a geologic carbon sequestration site through a linear feature such as a fault. We estimated the permeability of the soil and cobble layers present at the site by manual inversion of measurements of soil CO2 flux from a vertical-well CO2 release. Based on these estimated permeability values, predictive simulations for the horizontal well showed that CO2 injection just below the water table creates an effective gas-flow pathway through the saturated zone up to the unsaturated zone. Once in the unsaturated zone, CO2 spreads out laterally within the cobble layer, where liquid saturation is relatively low. CO2 also migrates upward into the soil layer through the capillary barrier and seeps out at the ground surface. The simulations predicted a breakthrough time of approximately two days for the 100kg d−1 injection rate, which also produced a flux within the range desired for testing detection and monitoring approaches. The seepage area produced by the model was approximately five meters wide above the horizontal well, compatible with the detection and monitoring methods tested. For a given flow rate, gas-phase diffusion of CO2 tends to dominate over advection near the ground surface, where the CO2 concentration gradient is large, while advection dominates deeper in the system.  相似文献   

15.
A model for the rheological properties of a concentrated suspension in weakly viscoelastic fluid matrices is proposed. The model is derived according to the Roscoe differential procedure described in 1952. The analytical results produced recently by Greco et al. (J Non-Newton Fluid Mech 147:1–10, 2007) and Housiadas and Tanner (J Non-Newton Fluid Mech 162:88–92, 2009) for dilute suspensions of neutrally buoyant, non-Brownian rigid spheres in weakly viscoelastic matrix fluids are the key results which are used as a base to predict the properties of concentrated suspensions. The results are compared with the few available experimental data from the literature, showing promising trends for the viscometric properties of the suspensions. In particular, one sees the rapidly increasing value of −N2/N1 as concentration increases.  相似文献   

16.
Predicting fluid replacement by two-phase flow in heterogeneous porous media is of importance for issues such as supercritical CO2 sequestration, the integrity of caprocks and the operation of oil water/brine systems. When considering coupled process modelling, the location of the interface is of importance as most of the significant interaction between processes will be happening there. Modelling two-phase flow using grid based techniques presents a problem as the fluid–fluid interface location is approximated across the scale of the discretisation. Adaptive grid methods allow the discretisation to follow the interface through the model, but are computationally expensive and make coupling to other processes (thermal, mechanical and chemical) complicated due to the constant alteration in grid size and effects thereof. Interface tracking methods have been developed that apply sophisticated reconstruction algorithms based on either the ratio of volumes of a fluid in an element (Volume of Fluid Methods) or the advective velocity of the interface throughout the modelling regime (Level set method). In this article, we present an “Analytical Front Tracking” method where a generic analytical solution for two-phase flow is used to “add information” to a finite element model. The location of the front within individual geometrical elements is predicted using the saturation values in the elements and the velocity field of the element. This removes the necessity for grid adaptation, and reduces the need for assumptions as to the shape of the interface as this is predicted by the analytical solution. The method is verified against a standard benchmark solution and then applied to the case of CO2 pooling and forcing its way into a heterogeneous caprock, replacing hot brine and eventually breaking through. Finally the method is applied to simulate supercritical CO2 injected into a brine saturated heterogeneous reservoir rock leading to significant viscous fingering and developement of preferential flow paths. The results are compared with to a finite volume simulation.  相似文献   

17.
On the basis of observations at four enhanced coalbed methane (ECBM)/CO2 sequestration pilots, a laboratory-scale study was conducted to understand the flow behavior of coal in a methane/CO2 environment. Sorption-induced volumetric strain was first measured by flooding fresh coal samples with adsorptive gases (methane and CO2). In order to replicate the CO2–ECBM process, CO2 was then injected into a methane-saturated core to measure the incremental “swelling.” As a separate effort, the permeability of a coal core, held under triaxial stress, was measured using methane. This was followed by CO2 flooding to replace the methane. In order to best replicate the conditions in situ, the core was held under uniaxial strain, that is, no horizontal strain was permitted during CO2 flooding. Instead, the horizontal stress was adjusted to ensure zero strain. The results showed that the relative strain ratio for CO2/methane was between 2 and 3.5. The measured volumetric strains were also fitted using a Langmuir-type model, thus enabling calculation of the strain at any gas pressure and using the analytical permeability models. For permeability work, effort was made to increase the horizontal stress to achieve the desired zero horizontal strain condition expected under in situ condition, but this became impossible because the “excess” stress required to maintain this condition was very large, resulting in sample failure. Finally, when CO2 was introduced and horizontal strain was permitted, permeability reduction was an order of magnitude greater, suggesting that the “excess” stress would have reduced it significantly further. The positive finding of the work was that the “excess” stresses associated with injection of CO2 are large. The excess stresses generated might be sufficient to cause microfracturing and increased permeability, and improved injectivity. Also, there might be a weakening effect resulting from repeated CO2 injection, as has been found to be the case with thermal cycling of rocks.  相似文献   

18.
According to the research theory of improved black oil simulator, a practical mathematical model for C02 miscible flooding was presented. In the model, the miscible process simulation was realized by adjusting oil/gas relative permeability and effective viscosity under the condition of miscible flow. In order to predict the production performance fast, streamline method is employed to solve this model as an alternative to traditional finite difference methods. Based on streamline distribution of steady-state flow through porous media with complex boundary confirmed with the boundary element method (BEM), an explicit total variation diminishing (TVD) method is used to solve the one-dimensional flow problem. At the same time, influences of development scheme, solvent slug size, and injection periods on CO2 drive recovery are discussed. The model has the advantages of less information need, fast calculation, and adaptation to calculate CO2 drive performance of all kinds of patterns in a random shaped porous media with assembly boundary. It can be an effective tool for early stage screening andmiscible oil field.reservoir dynamic management of the CO2 miscible oil field.  相似文献   

19.
In this paper, a non-isobaric Marangoni boundary layer flow that can be formed along the interface of immiscible nanofluids in surface driven flows due to an imposed temperature gradient, is considered. The solution is determined using a similarity solution for both the momentum and energy equations and assuming developing boundary layer flow along the interface of the immiscible nanofluids. The resulting system of nonlinear ordinary differential equations is solved numerically using the shooting method along with the Runge-Kutta-Fehlberg method. Numerical results are obtained for the interface velocity, the surface temperature gradient as well as the velocity and temperature profiles for some values of the governing parameters, namely the nanoparticle volume fraction φ (0≤φ≤0.2) and the constant exponent β. Three different types of nanoparticles, namely Cu, Al2O3 and TiO2 are considered by using water-based fluid with Prandtl number Pr =6.2. It was found that nanoparticles with low thermal conductivity, TiO2, have better enhancement on heat transfer compared to Al2O3 and Cu. The results also indicate that dual solutions exist when β<0.5. The paper complements also the work by Golia and Viviani (Meccanica 21:200–204, 1986) concerning the dual solutions in the case of adverse pressure gradient.  相似文献   

20.
Competitive Methane Desorption by Supercritical CO2 Injection in Coal   总被引:1,自引:0,他引:1  
A large diameter (∼70 mm) dry coal sample was used to study the competitive displacement of CH4 by injection of supercritical CO2, and CO2–CH4 counter-diffusion in coal matrix. During the test, a staged loading procedure, which allows the calibration of the key reservoir modelling parameters in a sequential and progressive manner, was employed. The core-flooding test was history matched using an Enhanced Coalbed Methane (ECBM) simulator, in which Fick’s Law for mixed gas diffusion and the extended Langmuir equations are implemented. The system pressure rise during the two loading stages and the CO2 breakthrough time in the final production stage were matched by using the pair of constant sorption times (9 and 3.2 days) for CH4 and CO2, respectively. The corresponding diffusion coefficients for CH4 and CO2 were estimated to be 1.6 ×  10−12 and 4.6 ×  10−12 m2/s, respectively. Comparison was made with published gas diffusion coefficients for dry ground samples (ranging from < 0.063 to ∼3 mm) of the same coal at relatively low pressures (< 4 MPa). The CO2/CH4 gas diffusion coefficient ratio was well within the reported range (2–3), whereas the CH4 diffusion coefficient obtained from history matching of the core-flooding test is approximately 15 times smaller than that arrived by curve-fitting the measured sorption uptake rate using a unipore diffusion model. The calibrated model prediction of the effluent gas composition was in good agreement with the test data for CO2 mole fraction of up to 20%.  相似文献   

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

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