Impact of difference accuracy on computational properties of vertical grids for a nonhydrostatic model

Starting from nonhydrostatic Boussinesq approximation equations, a general method is introduced to deduce the dispersion relationships. A comparative investigation is performed on inertia-gravity wave with horizontal lengths of 100, 10 and 1 km. These are examined using the second-order central difference scheme and the fourth-order compact difference scheme on vertical grids that are currently available from the perspectives of frequency, horizontal and vertical component of group velocity. These findings are compared to analytical solutions. The obtained results suggest that whether for the second-order central difference scheme or for the fourth-order compact difference scheme, Charny-Phillips and Lorenz ( L) grids are suitable for studying waves at the above-mentioned horizontal scales; the Lorenz time-staggered and Charny-Phillips time staggered (CPTS) grids are applicable only to the horizontal scales of less than 10 km, and N grid ( unstaggered grid) is unsuitable for simulating waves at any horizontal scale. Furthermore, by using fourth-order compact difference scheme with higher difference precision, the errors of frequency and group velocity in horizontal and vertical directions produced on all vertical grids in describing the waves with horizontal lengths of 1, 10 and 100 km cannot inevitably be decreased. So in developing a numerical model, the higher-order finite difference scheme, like fourth-order compact difference scheme, should be avoided as much as possible, typically on L and CPTS grids, since it will not only take many efforts to design program but also make the calculated group velocity in horizontal and vertical directions even worse in accuracy.

Impact of spatial resolution and spatial difference accuracy on the performance of Arakawa A-D grids

This paper alms at illustrating the impact of spatial difference scheme and spatial resolution on the performance of Arakawa A-D grids in physical space. Linear shallow water equations are discretized and forecasted on Arakawa A-D grids for 120-minute using the ordinary second-order (M and fourth-order (C4) finite difference schemes with the grid spacing being 100 km, 10 km and I km, respectively. Then the forecasted results are compared with the exact solution, the result indicates that when the grid spacing is I kin, the inertial gravity wave can be simulated on any grid with the same results from C2 scheme or C4 scheme, namely the impact of variable configuration is neglectable; while the inertial gravity wave is simulated with lengthened grid spacing, the effects of different variable configurations are different. However, whether for C2 scheme or for C4 scheme, the RMS is minimal (maximal) on C (D) grid. At the same time it is also shown that when the difference accuracy increases from C2 scheme to C4 scheme, the resulted forecasts do not uniformly decrease, which is validated by the change of the group A velocity relative error from C2 scheme to C4 scheme. Therefore, the impact of the grid spacing is more important than that of the difference accuracy on the performance of Arakawa A-D grid.

裂隙/孔隙型弹性介质中地震波场模拟与叠前逆时偏移

With the development of oil/gas seismic exploration, seismic survey for fracture/porosity type reservoir is becoming more and more important. As for China, since it has over 60% store of low porosity and low permeability oil/gas reservoir, it’s more urgent to validly describe fracture/porosity type oil/gas trap and proposing the related, developed seismic technique. To achieve mapping fracture/porosity region and its development status, it demands profound understanding of seismic wave propagation discipline in complex fractured/pored media. Meanwhile, it has profound scientific significance and applied worth to study forward modeling of fracture/porosity type media and pre-stacked reverse time migration. Especially, pre-stacked reverse-time migration is the lead edge technique in the field of seismology and seismic exploration. In this paper, the author has summarized the meaning, history and the present state of numerical simulation of seismic propagation in fractured/pored media and seismic exploration of fractured/pored reservoirs. Extensive Dilatancy Anisotropy (EDA) model is selected as media object in this work. As to forward modeling, due to local limitation of solving spatial partial derivative when using finite-difference and finite-element method, the author turns to pseudo-spectral method (PSM), which is based on the global characteristic of Fourier transform to simulate three-component elastic wave-field. Artifact boundary effect reduction and simulation algorithm stability are also discussed in the work. The author has completed successfully forward modeling coding of elastic wave-field and numerical simulation of two-dimensional and three-dimensional EDA models with different symmetric axis. Seismic dynamic and kinematical properties of EDA media are analyzed from time slices and seismic records of wave propagation. As to pre-stacked reverse-time migration for elastic wave-field in fractured/pored media, based on the successful experience in forward modeling results with PSM, the author has studied pre-stacked reverse-time depth-domain migration technique using PSM of elastic wave-field in two dimensional EDA media induced by preferred fracture/pore distribution. At the same time, different image conditions will bring up what kind of migration result is detailed in this paper. The author has worded out software for pre-stacked reverse-time depth-domain migration of elastic wave-field in EDA media. After migration processing of a series of seismic shot gathers, influences to migration from different isotropic and anisotropy models are described in the paper. In summary, following creative research achievements are obtained:  Realizing two-dimensional and three-dimensional elastic wave-field modeling for fractured/pored media and related software has been completed.  Proposed pre-stacked reverse-time depth-domain migration technique using PSM of elastic wave-field.  Through analysis of the seismic dynamic and kinematical properties of EDA media, the author made a conclusion that collection of multi-component seismic data can provide important data basis for locating and describing the fracture/pore regions and their magnitudes and the preferred directions.  Pre-stacked reverse-time depth-domain migration technique has the ability to reconstruct complex geological object with steep formations and tilt fracture distribution. Neglecting seismic anisotropy induced by the preferred fracture/pore distribution, will lead to the disastrous imaging results.

裂缝性油气藏储层预测方法研究

In this paper, based on the E & P situation in the oilfield and the theory of geophysical exploration, a series researches are conducted on fracture reservoir prediction technology in general，and it especially focus on some difficult points. The technological series which integrated amplitude preserved data processing、interpretation and its comprehensive application research as a whole were developed and this new method can be applied to the other similar oilfield exploration and development. The contents and results in this paper are listed as follows: 1. An overview was given on the status and development of fracture reservoir estimation technique, compare and analyze those geophysical prediction methods. This will be very helpful to the similar reservoir researches. 2. Analyze and conclude the characters of geologies and well logging response of burial hills fracture reservoir, those conclusions are used to steer the geophysical research and get satisfying results. 3. Forward modeling anisotropy seismic response of fracture reservoir. Quantitatively describe the azimuthal amplitude variation. Amplitude ellipse at each incidence angle is used to identify the fracture orientation. 4. Numerical simulation of structure stress based on finite difference method is carried out. Quantitatively describe and analyze the direction and intensity of fracture. 5. Conventional attributes extraction of amplitude preserved seismic data、attributes with different azimuthal angle and different offset are used to determine the relationship between the results and fracture distribution. 6. With spectrum decomposition method based on wavelet transform， the author disclose the reservoir distribution in space. It is a powerful tool to display its anisotropy. 7. Integrated seismic wave impendence、elastic impendence、spectrum decomposition、attribute extraction、fracture analysis result as a whole to identify and evaluate the fracture reservoir. An optimum workflow is constructed. It is used to practical oil&gas production and good results are obtained. This can indicate the wide foreground of this technique series.

二维自组织介质中地震波的传播特性

The topic of this study is simulation in the two dimensional self-organized media. The study in complexity of the earth plays an important role in structures, sources and energy seismic detection. And it mainly focuses on vertical or horizontal heterogeneous, anisotropic and linear media. Based on 2D self-organized velocity model and four-order finite-difference method, we simulate different types self-organized media and the same type mode with various parameters such as horizontal relative length, vertical relative length, variations, and velocity background gradient. Also we analyze the seismograms with complexity methods with instant information including amplitude, energy and frequency. The results can be summarized as the fallows: (1) The waveforms fluctuate with the velocity variations; (2) Different type self-organized media bring different effects on the amplitudes, energy and waveforms; (3) Different parameters also produce various influences to seismograms. (4) The layer contains their self-organized features, from which we can investigate the quality of the earth.

基于广义Lippmann-Schwinger方程的地震波场模拟算法与应用

In the last several decades, due to the fast development of computer, numerical simulation has been an indispensable tool in scientific research. Numerical simulation methods which based on partial difference operators such as Finite Difference Method (FDM) and Finite Element Method (FEM) have been widely used. However, in the realm of seismology and seismic prospecting, one usually meets with geological models which have piece-wise heterogeneous structures as well as volume heterogeneities between layers, the continuity of displacement and stress across the irregular layers and seismic wave scattering induced by the perturbation of the volume usually bring in error when using conventional methods based on difference operators. The method discussed in this paper is based on elastic theory and integral theory. Seismic wave equation in the frequency domain is transformed into a generalized Lippmann-Schwinger equation, in which the seismic wavefield contributed by the background is expressed by the boundary integral equation and the scattering by the volume heterogeneities is considered. Boundary element-volume integral method based on this equation has advantages of Boundary Element Method (BEM), such as reducing one dimension of the model, explicit use the displacement and stress continuity across irregular interfaces, high precision, satisfying the boundary at infinite, etc. Also, this method could accurately simulate the seismic scattering by the volume heterogeneities. In this paper, the concrete Lippmann-Schwinger equation is specifically given according to the real geological models. Also, the complete coefficients of the non-smooth point for the integral equation are introduced. Because Boundary Element-Volume integral equation method uses fundamental solutions which are singular when the source point and the field are very close，both in the two dimensional and the three dimensional case, the treatment of the singular kernel affects the precision of this method. The method based on integral transform and integration by parts could treat the points on the boundary and inside the domain. It could transform the singular integral into an analytical one both in two dimensional and in three dimensional cases and thus it could eliminate the singularity. In order to analyze the elastic seismic wave scattering due to regional irregular topographies, the analytical solution for problems of this type is discussed and the analytical solution of P waves by multiple canyons is given. For the boundary reflection, the method used here is infinite boundary element absorbing boundary developed by a pervious researcher. The comparison between the analytical solutions and concrete numerical examples validate the efficiency of this method. We thoroughly discussed the sampling frequency in elastic wave simulation and find that, for a general case, three elements per wavelength is sufficient, however, when the problem is too complex, more elements per wavelength are necessary. Also, the seismic response in the frequency domain of the canyons with different types of random heterogeneities is illustrated. We analyzed the model of the random media, the horizontal and vertical correlation length, the standard deviation, and the dimensionless frequency how to affect the seismic wave amplification on the ground, and thus provide a basis for the choice of the parameter of random media during numerical simulation.

离散裂缝诱导各向异性粘弹性地震波模拟与特征分析

The real media always attenuate and distort seismic waves as they propagate in the earth. This behavior can be modeled with a viscoelastic and anisotropic wave equation. The real media can be described as fractured media. In this thesis, we present a high-order staggered grid finite-difference scheme for 2-D viscoelastic wave propagation in a medium containing a large number of small finite length fractures. We use the effective medium approach to compute the anisotropic parameters in each grid cell. By comparing our synthetic seismogram by staggered-grid finite-difference with that by complex-ray parameter ray tracing method, we conclude that the high-order staggered-grid finite-difference technique can effectively used to simulate seismic propagation in viscoelastic-anisotropic media. Synthetic seismograms demonstrate that strong attenuation and significant frequency dispersion due to viscosity are important factors of reducing amplitude and delaying arrival time varying with incidence angle or offset. On the other hand, the amount of scattered energy not only provides an indicator of orientation of fracture sets, but can also provide information about the fracture spacing. Analysis of synthetic seismograms from dry- and fluid-filled fractures indicates that dry-filled fractures show more significant scattering on seismic wavefields than fluid-filled ones, and offset-variations in P-wave amplitude are observable. We also analyze seismic response of an anticlinal trap model that includes a gas-filled fractured reservoir with high attenuation, which attenuates and distorts the so-called bright spot.

地震P-S波时差耦合作用的斜坡崩滑效应研究

China locates between the circum-Pacific and the Mediterranean-Himalayan seismic belt. The seismic activities in our country are very frequent and so are the collapses and slides of slope triggered by earthquakes. Many collapses and slides of slope take place mainly in the west of China with many earthquakes and mountains, especially in Sichuan and Yunnan Provinces. When a strong earthquake happening, the damage especially in mountains area caused by geological hazards it triggered such as rock collapses, landslides and debris flows is heavier than that it caused directly. A conclusion which the number of lives lost caused by geological hazards triggered by a strong earthquake in mountains area often accounts for a half even more of the total one induced by the strong earthquake can be made by consulting the statistical loss of several representative earthquakes. As a result, geological hazards such as collapses and slides of slope triggered by strong earthquakes attract wide attention for their great costs. Based on field geological investigation, engineering geological exploration and material data analysis, chief conclusions have been drawn after systematic research on formation mechanism, key inducing factors, dynamic characteristics of geological hazards such as collapses and slides of slope triggered by strong earthquakes by means of engineering geomechanics comprehensive analysis, finite difference numerical simulation test, in-lab dynamic triaxial shear test of rock, discrete element numerical simulation. Based on research on a great number of collapses and landslides triggered by Wenchuan and Xiaonanhai Earthquake, two-set methods, i.e. the method for original topography recovering based on factors such as lithology and elevation comparing and the method for reconstructing collapsing and sliding process of slope based on characteristics of seism tectonic zone, structural fissure, diameter spatial distribution of slope debris mass, propagation direction and mechanical property of seismic wave, have been gotten. What is more, types, formation mechanism and dynamic characteristics of collapses and slides of slope induced by strong earthquakes are discussed comprehensively. Firstly, collapsed and slided accumulative mass is in a state of heavily even more broken. Secondly, dynamic process of slope collapsing and sliding consists of almost four stages, i.e. broken, thrown, crushed and river blocked. Thirdly, classified according to failure forms, there are usually four types which are made up of collapsing, land sliding, land sliding-debris flowing and vibrating liquefaction. Finally, as for key inducing factors in slope collapsing and sliding, they often include characteristics of seism tectonic belts, structure and construction of rock mass, terrain and physiognomy, weathering degree of rock mass and mechanical functions of seismic waves. Based on microscopic study on initial fracturing of slope caused by seismic effect, combined with two change trends which include ratio of vertical vs. horizontal peak ground acceleration corresponding to epicentral distance and enlarging effect of peak ground acceleration along slope, key inducing factor of initial slope fracturing in various area with different epicentral distance is obtained. In near-field area, i.e. epicentral distance being less than 30 km, tensile strength of rock mass is a key intrinsic factor inducing initial fracturing of slope undergoing seismic effect whereas shear strength of rock mass is the one when epicentral distance is more than 30 km. In the latter circumstance, research by means of finite difference numerical simulation test and in-lab dynamic triaxial shear test of rock shows that initial fracture begins always in the place of slope shoulder. The fact that fracture strain and shear strength which are proportional to buried depth of rock mass in the place of slope shoulder are less than other place and peak ground acceleration is enlarged in the place causes prior failure at slope shoulder. Key extrinsic factors inducing dynamic fracture of slope at different distances to epicenter have been obtained through discrete element numerical simulation on the total process of collapsing and sliding of slope triggered by Wenchuan Earthquake. Research shows that combined action of P and S seismic waves is the key factor inducing collapsing and sliding of slope at a distance less than 64 km to initial epicenter along earthquake-triggering structure. What is more, vertical tensile action of P seismic wave plays a leading role near epicenter, whereas vertical shear action of S seismic wave plays a leading role gradually with epicentral distance increasing in this range. On the other hand, single action of P seismic wave becomes the key factor inducing collapsing and sliding of slope at a distance between 64 km and 216 km to initial epicenter. Horizontal tensile action of P seismic wave becomes the key factor gradually from combined action between vertical and horizontal tensile action of P seismic wave with epicentral distance increasing in this distance range. In addition, initial failure triggered by strong earthquakes begins almost in the place of slope shoulder. However, initial failure beginning from toe of slope relates probably with gradient and rock occurrence. Finally, starting time of initial failure in slope increases usually with epicentral distance. It is perhaps that the starting time increasing is a result of attenuating of seismic wave from epicenter along earthquake-triggering structure. It is of great theoretical and practical significance for us to construct towns and infrastructure in fragile geological environment along seism tectonic belts and conduct risk management on earthquake-triggered geological hazards by referring to above conclusions.

三维油藏数值模拟正反演算法研究

The dynamic prediction of complex reservoir development is one of the important research contents of dynamic analysis of oil and gas development. With the increase development of time, the permeabilities and porosities of reservoirs and the permeability of block reservoir at its boundaries are dynamically changing. How to track the dynamic change of permeability and porosity and make certain the permeability of block reservoir at its boundary is an important practical problem. To study developing dynamic prediction of complex reservoir, the key problem of research of dynamic prediction of complex reservoir development is realizing inversion of permeability and porosity. To realize the inversion, first of all, the fast forward and inverse method of 3-dimension reservoir simulation must be studied. Although the inversion has been widely applied to exploration and logging, it has not been applied to3-dimension reservoir simulation. Therefore, the study of fast forward and inverse method of 3-dimension reservoir simulation is a cutting-edge problem, takes on important realistic signification and application value. In this dissertation, 2-dimension and 3-dimension fluid equations in porous media are discretized by finite difference, obtaining finite difference equations to meet the inner boundary conditions by Peaceman's equations, giving successive over relaxation iteration of 3-dimension fluid equations in porous media and the dimensional analysis. Several equation-solving methods are compared in common use, analyzing its convergence and convergence rate. The alternating direction implicit procedure of 2-dimension has been turned into successive over relaxation iteration of alternating direction implicit procedure of 3-dimension fluid equations in porous media, which possesses the virtues of fast computing speed, needing small memory of computer, good adaptability for heterogeneous media and fast convergence rate. The geological model of channel-sandy reservoir has been generated with the help of stochastic simulation technique, whose cross sections of channel-sandy reservoir are parabolic shapes. This method makes the hard data commendably meet, very suit for geological modeling of containing complex boundary surface reservoir. To verify reliability of the method, theoretical solution and numerical solution are compared by simplifying model of 3-dimension fluid equations in porous media, whose results show that the only difference of the two pressure curves is that the numerical solution is lower than theoretical at the wellbore in the same space. It proves that using finite difference to solve fluid equations in porous media is reliable. As numerical examples of 3-dimension heterogeneous reservoir of the single-well and multi-well, the pressure distributions have been computed respectively, which show the pressure distributions there are clearly difference as difference of the permeabilities is greater than one order of magnitude, otherwise there are no clearly difference. As application, the pressure distribution of the channel-sandy reservoir have been computed, which indicates that the space distribution of pressure strongly relies on the direction of permeability, and is sensitive for space distributions of permeability. In this dissertation, the Peaceman's equations have been modified into solving vertical well problem and horizontal well problem simultaneously. In porous media, a 3D layer reservoir in which contain vertical wells and horizontal wells has been calculated with iteration. For channel-sandy reservoir in which there are also vertical wells and horizontal wells, a 3D transient heterogeneous fluid equation has been discretized. As an example, the space distribution of pressure has been calculated with iteration. The results of examples are accord with the fact, which shows the modification of Peaceman's equation is correct. The problem has been solved in the space where there are vertical and horizontal wells. In the dissertation, the nonuniform grid permeability integration equation upscaling method, the nonuniform grid 2D flow rate upscaling method and the nonuniform grid 3D flow rate upscaling method have been studied respectively. In those methods, they enhance computing speed greatly, but the computing speed of 3D flow rate upscaling method is faster than that of 2D flow rate upscaling method, and the precision of 3D flow rate upscaling method is better than that of 2D flow rate upscaling method. The results also show that the solutions of upscaling method are very approximating to that of fine grid blocks. In this paper, 4 methods of fast adaptive nonuniform grid upscaling method of 3D fluid equations in porous media have been put forward, and applied to calculate 3D heterogeneous reservoir and channel-sandy reservoir, whose computing results show that the solutions of nonuniform adaptive upscaling method of 3D heterogeneous fluid equations in porous media are very approximating to that of fine grid blocks in the regions the permeability or porosity being abnormity and very approximating to that of coarsen grid blocks in the other region, however, the computing speed of adaptive upscaling method is 100 times faster than that of fine grid block method. The formula of sensitivity coefficients are derived from initial boundary value problems of fluid equations in porous media by Green's reciprocity principle. The sensitivity coefficients of wellbore pressure to permeability parameters are given by Peaceman's equation and calculated by means of numerical calculation method of 3D transient anisotropic fluid equation in porous media and verified by direct method. The computing results are in excellent agreement with those obtained by the direct method, which shows feasibility of the method. In the dissertation, the calculating examples are also given for 3D reservoir, channel-sandy reservoir and 3D multi-well reservoir, whose numerical results indicate: around the well hole, the value of the sensitivity coefficients of permeability is very large, the value of the sensitivity coefficients of porosity is very large too, but the sensitivity coefficients of porosity is much less than the sensitivity coefficients of permeability, so that the effect of the sensitivity coefficients of permeability for inversion of reservoir parameters is much greater than that of the sensitivity coefficients of porosity. Because computing the sensitivity coefficients needs to call twice the program of reservoir simulation in one iteration, realizing inversion of reservoir parameters must be sustained by the fast forward method. Using the sensitivity coefficients of permeability and porosity, conditioned on observed valley erosion thickness in wells (hard data), the inversion of the permeabilities and porosities in the homogeneous reservoir, homogeneous reservoir only along the certain direction and block reservoir are implemented by Gauss-Newton method or conjugate gradient method respectively. The results of our examples are very approximating to the real data of permeability and porosity, but the convergence rate of conjugate gradient method is much faster than that of Gauss-Newton method.

单（多）井抽灌对浅部地温场的影响研究

The technique of energy extraction using groundwater source heat pumps, as a sustainable way of low-grade thermal energy utilization, has widely been used since mid-1990's. Based on the basic theories of groundwater flow and heat transfer and by employing two analytic models, the relationship of the thermal breakthrough time for a production well with the effect factors involved is analyzed and the impact of heat transfer by means of conduction and convection, under different groundwater velocity conditions, on geo-temperature field is discussed.A mathematical model, coupling the equations for groundwater flow with those for heat transfer, was developed. The impact of energy mining using a single well system of supplying and returning water on geo-temperature field under different hydrogeological conditions, well structures, withdraw-and-reinjection rates, and natural groundwater flow velocities was quantitatively simulated using the finite difference simulator HST3D. Theoretical analyses of the simulated results were also made. The simulated results of the single well system indicate that neither the permeability nor the porosity of a homogeneous aquifer has significant effect on the temperature of the production segment provided that the production and injection capability of each well in the aquifers involved can meet the designed value. If there exists a lower permeable interlayer, compared with the main aquifer, between the production and injection segments, the temperature changes of the production segment will decrease. The thicker the interlayer and the lower the interlayer permeability, the longer the thermal breakthrough time of the production segment and the smaller the temperature changes of the production segment. According to the above modeling, it can also be found that with the increase of the aquifer thickness, the distance between the production and injection screens, and/or the regional groundwater flow velocity, and/or the decrease of the production-and-reinjection rate, the temperature changes of the production segment decline. For an aquifer of a constant thickness, continuously increase the screen lengths of production and injection segments may lead to the decrease of the distance between the production and injection screens, and the temperature changes of the production segment will increase, consequently.According to the simulation results of the single well system, the parameters, that can cause significant influence on heat transfer as well as geo-temperature field, were chosen for doublet system simulation. It is indicated that the temperature changes of the pumping well will decrease as the aquifer thickness, the distance between the well pair and/or the screen lengths of the doublet increase. In the case of a low permeable interlayer embedding in the main aquifer, if the screens of the pumping and the injection wells are installed respectively below and above the interlayer, the temperature changes of the pumping well will be smaller than that without the interlay. The lower the permeability of the interlayer, the smaller the temperature changes. The simulation results also indicate that the lower the pumping-and-reinjection rate, the greater the temperature changes of the pumping well. It can also be found that if the producer and the injector are chosen reasonably, the temperature changes of the pumping well will decline as the regional groundwater flow velocity increases. Compared with the case that the groundwater flow direction is perpendicular to the well pair, if the regional flow is directed from the pumping well to the injection well, the temperature changes of the pumping well is relatively smaller.Based on the above simulation study, a case history was conducted using the data from an operating system in Beijing. By means of the conceptual model and the mathematical model, a 3-D simulation model was developed and the hydrogeological parameters and the thermal properties were calibrated. The calibrated model was used to predict the evolution of the geo-temperature field for the next five years. The simulation results indicate that the calibrated model can represent the hydrogeological conditions and the nature of the aquifers. It can also be found that the temperature fronts in high permeable aquifers move very fast and the radiuses of temperature influence are large. Comparatively, the temperature changes in clay layers are smaller and there is an obvious lag of the temperature changes. According to the current energy mining load, the temperature of the pumping wells will increase by 0.7°C at the end of the next five years. The above case study may provide reliable base for the scientific management of the operating system studied.

Influence of Local Magnetic Fields on P-Doped Si Floating Zone Melting Crystal Growth in Microgravity

A two-dimensional axisymmetric numerical model is presented to study the influence of local magnetic fields on P-doped Si floating zone melting crystal growth in microgravity. The model is developed based on the finite difference method in a boundary-fitted curvilinear coordinate system. Extensive numerical simulations are carried out, and parameters studied include the curved growth interface shape and the magnetic field configurations. Computed results show that the local magnetic field is more effective in reducing the impurity concentration nonuniformity at the growth interface in comparison with the longitudinal magnetic field. Moreover, the curved growth interface causes more serious impurity concentration nonuniformity at the growth interface than the case with a planar growth interface.

对流扩散方程的变步长摄动有限差分格式

摄动有限差分(PFD)方法是构造高精度差分格式的一种新方法。变步长摄动有限差分方法是等步长摄动有限差分方法的发展和推广。对需要局部加密网格的计算问题，变步长PFD格式不需要对自变量进行数学变换，且和等步长PFD格式一样，具有如下的共同特点：从变步长一阶迎风格式出发，通过把非微商项(对流系数和源项)作变步长摄动展开，展开幂级数系数通过消去摄动格式修正微分方程的截断误差项求出，由此获得高精度变步长PFD格式。该格式在一、二和三维情况下分别仅使用三、五和七个基点，且具有迎风性。文中利用变步长PFD格式对对流扩散反应模型方程，变系数方程及Burgers方程等进行了数值模拟，并与一阶迎风和二阶中心格式及其问题的精确解作了比较。数值试验表明，与一阶迎风和二阶中心格式相比，变步长PFD格式具有精度高，稳定性与收敛性好的特点。变步长PFD格式与等步长PFD格式相比，变步长PFD解在薄边界层型区域的分辨率得到了明显的提高。

A Numerical Study of Indentation Using Indenters of Different Geometry

Finite element simulation of the Berkovich, Vickers, Knoop, and cone indenters was carried out for the indentation of elastic-plastic material. To fix the semiapex angle of the cone, several rules of equivalence were used and examined. Despite the asymmetry and differences in the stress and strain fields, it was established that for the Berkovich and Vickers indenters, the load-displacement relation can closely be simulated by a single cone indenter having a semiapex angle equal to 70.3degrees in accordance with the rule of the volume equivalence. On the other hand, none of the rules is applicable to the Knoop indenter owing to its great asymmetry. The finite element method developed here is also applicable to layered or gradient materials with slight modifications.

双曲守恒型方程的二阶摄动有限差分格式

对双曲守恒型方程，将其一阶迎风格式空间差商的常系数摄动展开为时间步长和空间步长的幂级数，通过确定幂级数系数而获得二阶精度的摄动有限差分（PFD）格式。进而从双曲守恒型方程的通量分裂型一阶迎风格式出发，通过娄似的摄动展开方法，获得空间精度为二阶的通量分裂形式的摄动有限差分（FPFD）格式。这两类格式保留了一阶守恒迎风格式的简洁结构形式，使用三节点即可达到二阶精度，又避免了三点二阶格式的非物理数值振荡。并将这两类格式推广应用到双曲守恒型方程组，最后通过模型方程和一维激波管流动的数值算例验证了格式的高精度、高分辨率性质。

NS方程激波计算的摄动有限差分方法

摄动有限差分（PFD）方法从一阶迎风差分格式出发，将差分系数展开为网格步长的幂级数，通过提高修正微分方程的逼近精度来获得更高精度的差分格式。由于格式基于一阶迎风格式，因此具有迎风效应、网格节点少等特点。本文首先通过对Burgers方程的摄动差分格式的推导，将摄动有限差分格式引入时间相关法的计算，并构造了守恒形式的摄动有限差分格式，然后推广到一维Navier-Stokes方程组的计算。数值比较研究表明：本文构造的NS方程摄动有限差分格式具有比一阶迎风较高的精度和分辨率，而且保持了一阶迎风格式的无振荡性质。