Self-demagnetization due to strongly magnetic bodies can seriously affect the interpretation of magnetic anomalies. Conventional numerical methods often neglect the self-demagnetization effects and limit their use to low susceptibilities (χ<1 SI). We have developed a novel iterative method based on the integral equation and the Gauss-fast Fourier transform (FFT) technique for calculating the magnetic field from finite bodies of high magnetic susceptibility and arbitrary shapes. The method uses a segmented model consisting of prismatic voxels to approximate a complex target region. In each voxel, the magnetization is assumed to be constant, so that the integral equation in the spatial domain can reduce to a simple form with lots of merit in numerical calculation after the 2D Fourier transform. Moreover, a contraction operator is derived to ensure the convergence of the iterative calculation, and the Gauss-FFT technique is applied to reduce numerical errors due to edge effects. Our modeling results indicate that this new iterative scheme performs well in a wide range of magnetic susceptibilities (1–1000 SI). For lower susceptibilities (χ≤10 SI), the iterative algorithm converges rapidly and indicates very good correlation with the analytical solutions. At higher susceptibilities (10<χ≤100 SI), our method still performs well, but the numerical accuracy improves with a relatively slow speed. In the extreme case (χ=1000 SI), an acceptable result is also obtained after sufficient iterative computation. A further improvement in the numerical precision can be achieved by increasing the number of prismatic voxels, but at the same time, the computational time increases linearly with the size of the voxels.

To take the sphericity of the Earth into account, tesseroids are often utilized as grid elements in large-scale gravitational forward modeling. However, such elements in a latitude–longitude mesh suffer from degenerating into poorly shaped triangles near poles. Moreover, tesseroids have limited flexibility in describing laterally variable density distributions with irregular boundaries and also face difficulties in achieving completely equivalent division over a spherical surface that may be desired in a gravity inversion. We develop a new method based on triangular spherical prisms (TSPs) for 3D gravitational modeling in spherical coordinates. A TSP is defined by two spherical surfaces of triangular shape, with one of which being the radial projection of the other. Due to the spherical triangular shapes of the upper and lower surfaces, TSPs enjoy more advantages over tesseroids in describing mass density with different lateral resolutions. In addition, such an element also allows subdivisions with nearly equal weights in spherical coordinates. To calculate the gravitational effects of a TSP, we assume the density in each element to be polynomial along radial direction so as to accommodate a complex density environment. Then, we solve the Newton’s volume integral using a mixed Gaussian quadrature method, in which the surface integral over the spherical triangle is calculated using a triangle-based Gaussian quadrature rule via a radial projection that transforms the spherical triangles into linear ones. A 2D adaptive discretization strategy and an extension technique are also combined to improve the accuracy at observation points near the mass sources. The numerical experiments based on spherical shell models show that the proposed method achieves good accuracy from near surface to a satellite height in the case of TSPs with various dimensions and density variations. In comparison with the classical tesseroid-based method, the proposed algorithm enjoys better accuracy and much higher flexibility for density models with laterally irregular shapes. It shows that to achieve the same accuracy, the number of elements required by the proposed method is much less than that of the tesseroid-based method, which substantially speeds up the calculation by more than 2 orders. The application to the tessellated LITHO1.0 model further demonstrates its capability and practicability in realistic situations. The new method offers an attractive tool for gravity forward and inverse problems where the irregular grids are involved.

Fast and accurate calculation of gravitational effects on a regional or global scale with complex density environment is a critical issue in gravitational forward modelling. Most existing significant developments with tessroid-based modelling are limited to homogeneous density models or polynomial ones of a limited order. Moreover, the total gravitational effects of tesseroids are often calculated by pure summation in these methods, which makes the calculation extremely time-consuming. A new efficient and accurate method based on tesseroids with a polynomial density up to an arbitrary order in depth is developed for 3D large-scale gravitational forward modelling. The method divides the source region into a number of tesseroids, and the density in each tesseroid is assumed to be a polynomial function of arbitrary degree. To guarantee the computational accuracy and efficiency, two key points are involved: (1) the volume Newton’s integral is decomposed into a one-dimensional integral with a polynomial density in the radial direction, for which a simple analytical recursive formula is derived for efficient calculation, and a surface integral over the horizontal directions evaluated by the Gauss–Legendre quadrature (GLQ) combined with a 2D adaptive discretization strategy; (2) a fast and flexible discrete convolution algorithm based on 1D fast Fourier transform (FFT) and a general Toepritz form of weight coefficient matrices is adopted in the longitudinal dimension to speed up the computation of the cumulative contributions from all tesseroids. Numerical examples show that the gravitational fields predicted by the new method have a good agreement with the corresponding analytical solutions for spherical shell models with both polynomial and non-polynomial density variations in depth. Compared with the 3D GLQ methods, the new algorithm is computationally more accurate and efficient. The calculation time is significantly reduced by 3 orders of magnitude as compared with the traditional 3D GLQ methods. Application of the new algorithm in the global crustal CRUST1.0 model further verifies its reliability and practicability in real cases. The proposed method will provide a powerful numerical tool for large-scale gravity modelling and also an efficient forward engine for inversion and continuation problems.

Wall thinning, as one of the key defects determined directly the usability of the bent tube, should be strictly controlled. Elastic modulus changes with the plastic deformation and its influence the plastic forming quality of tube bending. To precisely predict the wall thinning of tube bending, a finite element (FE) model for numerical control (NC) bending of high strength 21-6-9 stainless steel tube (21-6-9-HS tube) was established considering the variation of elastic modulus. Using the model, the effects of process parameters on wall thinning of the 21-6-9-HS tube in NC bending were investigated. The results show that the variation of elastic modulus has no obvious effect on the change tendency of the wall thinning degree, but only increases the value of the wall thinning degree. The wall thinning degree enhances sharply at first, then occurs a platform deforming characteristics, and finally reduces abruptly from the bending section to the initial bending section, and it increases with the increase of the clearance of tube-bending die, the friction coefficient of tube-bending die and the friction coefficient of tube-wiper die or with the decrease of the clearance of tube-wiper die, but the increased degree is not significant in all conditions.

The magnetic gradient tensor provides a powerful tool for detecting magnetic bodies because of its ability to emphasize detailed features of the magnetic anomalies. To interpret field measurements obtained by magnetic gradiometry, the forward calculation of magnetic gradient fields is always necessary. In this paper, we present a contraction integral equation method to simulate the gradient fields produced by 3-D magnetic bodies of arbitrary shapes and high susceptibilities. The method employs rectangular prisms to approximate the source region with the assumption that the magnetization in each element is homogeneous. The gradient fields are first solved in the Fourier domain and then transformed into the spatial domain by 2-D Gauss-FFT. This calculation is performed iteratively until the required accuracy is reached. The convergence of the iterative procedure is ensured by a contraction operator. To facilitate application, we introduce a FORTRAN program to implement the algorithm. This program is intended for users who show interests in 3D magnetic modeling at high susceptibility. The performance of the program, including its computational accuracy, efficiency and convergence behavior, is tested by several models. Numerical results show that the code is computationally accurate and efficient, and performs well at a wide range of magnetic susceptibilities from 0 SI to 1000 SI. This work, therefore, provides a significant tool for 3D forward modeling of magnetic gradient fields at high susceptibility.

Multidimensional Fourier transform on an irregular grid is a useful tool for various seismic forward problems caused by complex media and wavefield distributions. Using a shape-function-based strategy, we have developed four different algorithms for 1D and 2D nonuniform Fourier transforms, such as two high-accuracy Fourier transforms (the linear shape-function-based Fourier transform [LSF-FT] and the quadratic shape-function-based Fourier transform [QSF-FT]) and two nonuniform fast Fourier transforms (NUFFTs) (the linear shape-function-based NUFFT [LSF-NUFFT] and the quadratic shape-function-based NUFFT [QSF-NUFFT]), respectively, based on their linear and quadratic shape functions. The main advantage of incorporating shape functions into the Fourier transform is that triangular elements can be used to mesh any complex wavefield distribution in the 2D case. Therefore, these algorithms can be used in conjunction with any irregular sampling strategies. The accuracy and efficiency of the four nonuniform Fourier transforms are investigated and compared by applying them in frequency-domain seismic wave modeling. All of the algorithms are compared with the exact solutions. Numerical tests indicate that the quadratic shape-function-based algorithms are more accurate than those based on the linear shape function. Moreover, LSF-FT/QSF-FT exhibits higher accuracy but much slower calculation speed, whereas LSF-NUFFT/QSF-NUFFT is highly efficient but has lower accuracy at near-source points. In contrast, a combination of these algorithms, by using QSF-FT at near-source points and LSF-NUFFT/QSF-NUFFT at others, achieves satisfactory efficiency and high accuracy at all points. Although our tests are restricted to seismic models, these improved NUFFT algorithms may also have potential applications in other geophysical problems, such as forward modeling in complex gravity and magnetic models.