# Solution of Axisymmetric Potential Problem in Oblate Spheroid Using the Exodus Method

Authors: O. D. Momoh; M. N. O. Sadiku; S. M. Musa
Journal: Journal of Computational Engineering
 (2014)

---

## Abstract

This paper presents the use of Exodus method for computing potential distribution within a conducting oblate spheroidal system. An explicit finite difference method for solving Laplace’s equation in oblate spheroidal coordinate systems for an axially symmetric geometry was developed. This was used to determine the transition probabilities for the Exodus method. A strategy was developed to overcome the singularity problems encountered in the oblate spheroid pole regions. The potential computation results obtained correlate with those obtained by exact solution and explicit finite difference methods.


---

## Body

## 1. Introduction

An oblate spheroid is the surface generated by the rotation of an ellipse about its minor axis, and depending upon the ellipse’s eccentricity, the spheroid will be flattened about the minor axis [1]. An oblate spheroidal shell, for instance, is considered as a continuous system constructed from two spherical shell caps by matching the continuous boundary conditions [2].Oblate and prolate spheroidal coordinates are widely used in many fields of science and engineering, such as potential theory, fluid mechanics, heat and mass transfer, thermal stress, and elastic inclusions. For example, oblate and prolate spheroids being surfaces of revolution can be more easily conformed to most districts of human body (e.g., extremities) which is of interest for dedicated MRI systems [3]. Oblate spheroidal coordinates are the natural choice for the translation of any ellipsoid parallel to a principal axis [4]. There is a more recent improvement in the lightning ground tracking systems based on the time-of-arrival (TOA) technique because of the refinement in the mathematics to more accurately accommodate the oblate shape of the earth spheroid. Approximating the earth as a perfect sphere affects not only the accuracy of time clock offset calculations, but also the accuracy of stroke coordinate computation given receiver time differences. Oblate solution mathematics can provide a substantial systematic error reduction of up to 50% percent [5].In this paper, Exodus method is used to compute potential distribution inside conducting oblate spheroidal shells maintained at two potentials. This work is a continuation of our previous work in which a fixed random walk Monte Carlo method (MCM) was used for numerical computation of potential distribution with two conducting oblate spheroidal shells maintained at two potentials [6]. An explicit Neumann boundary condition was imposed at the pole regions (η=-π/2,π/2) of the oblate spheroid to treat the presence of singularities in those regions.The Exodus method is one of Monte Carlo methods which are nondeterministic (probabilistic or stochastic) numerical methods employed in solving mathematical and physical problems. The fixed random walk and Exodus methods are the most frequently used Monte Carlo methods for solving heat conduction and potential problems. The Exodus method is however preferred to fixed random walk method because of its computational efficiency. It yields more accurate results with less computing time as compared to the original Monte Carlo method [7, 8].


## 2. Oblate Spheroidal Coordinate Systems

The geometry of an oblate spheroid showing surfaces of constant oblate spheroidal coordinates is illustrated in Figure1. The oblate spheroidal coordinates are related to the rectangular coordinates as follows:
(1)x=c′cosh⁡ξsinηcos⁡φ,y=c′cosh⁡ξsinηsinφ,x=c′sinh⁡ξcos⁡η,
where c′ is the focal length of the oblate spheroid.Figure 1
The oblate spheroidal coordinate system [9].An arbitrary grid pointr on the oblate spheroid is given by
(2)r-i,j,k=ξia-ξ+ηja-η+ϕka-ϕ,
where a-ξ, a-η, and a-ϕ are the unit vectors in the direction of ξ, η, and ϕ coordinates, respectively.The respective scale factor for each of the three coordinates(ξ,η,ϕ) is
(3)hξ=c′sinh⁡2ξ+sin2η,hη=c′sinh⁡2ξ+sin2η,hϕ=c′cosh⁡ξcos⁡η.


## 3. Potential Distribution Computation

### 3.1. Finite Difference Transformation of Oblate Spheroid Laplace’s Equation

The Laplacian equation in oblate spheroidal coordinate systems is(4)0=∇2V=(1sinh⁡2ξ+sin2η)*[(sech2ξtan2η+sec2ηtanh2ξ)∂2V∂ϕ2+∂2V∂ξ2+tanhξ∂V∂ξ+∂2V∂η2-tanη∂V∂η].
The term outside the square bracket may be ignored. Also the first term inside square bracket is ignored due to the rotational symmetry about the vertical (z) axis. Therefore, (4) reduces to
(5)∂2V∂ξ2+tanhξ∂V∂ξ+∂2V∂η2-tanη∂V∂η=0.
Equation (5) governs potential distribution in an axisymmetric oblate spheroid potential problem.Two oblate spheroidal shells made up of two constant conducting surfacesξ1 and ξ2 are shown in Figure 2. The two equipotential surfaces are maintained at 50 V and 100 V (Dirichlet boundary conditions), respectively. The choice of these Dirichlet boundary conditions is arbitrary. Any potential value can be assigned. Also, the value of the constant oblate spheroidal surfaces equipotential that constitutes the two conducting shells are arbitrarily chosen as ξ1=1.0 and ξ1=2.0, respectively.Figure 2
Two oblate spheroidal shells.The explicit finite difference transformation of (5) is
(6)V(i,j)=0.5[1[(1/(Δξ)2)+(1/(Δη)2)]]*[[1(Δη)2+tan(ηj)2Δη]V(i+1,j)[1(Δξ)2+tanh(iΔξ)2Δξ]+V(i-1,j)[1(Δξ)2-tanh(iΔξ)2Δξ]+V(i,j+1)[1(Δη)2-tan(ηj)2Δη]+V(i,j-1)[1(Δη)2+tan(ηj)2Δη]].
Figure 3 shows one-quarter of the constant oblate spheroidal surfaces. The figure exhibits symmetry with respect to the η coordinate. Therefore, two lines of symmetries will be encountered in this range of η. They are η=90° and η=0°. On these lines of symmetries, the condition ∂V/∂η=0 is imposed. This strategy eliminates the singularity causing term (tan(ηj)/2Δη) at the oblate spheroid poles as seen in (6).Figure 3
Oblate spheroidal surface path used in computing the charge enclosed.Consequently, the finite difference equations along the two lines of symmetries become as follows.Alongη=90°, j=jmax⁡,
(7)V(i,jmax⁡)=0.5[1[(1/(Δξ)2)+(1/(Δη)2)]]*[[1(Δη)2]V(i+1,jmax⁡)[1(Δξ)2+tanh⁡(iΔξ)2Δξ]+V(i-1,jmax⁡)[1(Δξ)2-tanh(iΔξ)2Δξ]+2V(i,jmax⁡-1)[1(Δη)2]].Alongη=0°, j=0,
(8)V(i,0)=0.5[1[(1/(Δξ)2)+(1/(Δη)2)]]*[[1(Δη)2]V(i+1,0)[1(Δξ)2+tanh⁡(iΔξ)2Δξ]+V(i-1,0)[1(Δξ)2-tanh⁡(iΔξ)2Δξ]+2V(i,1)[1(Δη)2]].


### 3.2. Transition Probabilities Determination

In a more compact form, (6) can be rewritten as
(9)V(i,j)=Pξ+V(i+1,j)+Pξ-V(i-1,j)+Pη+V(i,j+1)+Pη-V(i,j-1),
where(10a)Pξ-=0.5([(1/(Δξ)2)+(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(10b)Pξ-=0.5([(1/(Δξ)2)-(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(10c)Pη-=0.5([(1/(Δη)2)+(tan⁡(jΔη)/2Δη)][(1/(Δξ)2)+(1/(Δη)2)]),(10d)Pη-=0.5([(1/(Δη)2)-(tan⁡(jΔη)/2Δη)][(1/(Δξ)2)+(1/(Δη)2)]).Note that Pξ++Pξ-+Pη++Pη-=1. Equations (10a), (10b), (10c), and (10d) serve as the transition probabilities for the Exodus method used in this work. If N particles are dispersed at node (ξ0,η0), they have probabilities Pξ+,Pξ-,Pη+, andPη- of moving to points (ξ+Δξ,η),(ξ-Δξ,η),(ξ,η+Δη), and (ξ,η-Δη), respectively. The direction of movement of the particles is determined by generating a random number U,0<U<1. The particles are instructed to walk as follows:
(11)(ξ0,η0)⟹(ξ0+Δξ,η0),if(0<U<Pξ+),(ξ0,η0)⟹(ξ0-Δξ,η0),if(Pξ+<U<PT),(ξ0,η0)⟹(ξ0,η0+Δη),if(PT<U<PTT),(ξ0,η0)⟹(ξ0,η0-Δη),if(PTT<U<1),
where
(12)PT=Pξ++Pξ-,PTT=Pξ++Pξ-+Pη+.For the random walk computation along the line of symmetry (η=90°), (7) becomes
(13)V(i,jmax⁡)=Pξ+V(i+1,jmax⁡)+Pξ-V(i-1,jmax⁡)+Pη-V(i,jmax⁡-1),
where(14a)Pξ-=0.5([(1/(Δξ)2)+(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(14b)Pξ-=0.5([(1/(Δξ)2)-(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(14c)Pη-=([(1/(Δη)2)][(1/(Δξ)2)+(1/(Δη)2)]).Note that Pξ++Pξ-+Pη-=1 along this line of symmetry. If a random-walking particle is instantaneously at the point (ξ0,90°), it has probabilities Pξ+,Pξ-, and Pη- of moving to points (ξ+Δξ,90°),(ξ-Δξ,90°), and (ξ,90°-Δη), respectively. The particles are instructed to walk as follows:
(15)(ξ0,90°)⟹(ξ0+Δξ,90°),if(0<U<Pξ+),(ξ0,90°)⟹(ξ0-Δξ,90°),if(Pξ+<U<Pmax⁡),(ξ0,90°)⟹(ξ0,90°-Δη),if(Pmax⁡<U<1),
where
(16)Pmax⁡=Pξ++Pξ-.For the random walk computation along the line of symmetry(ξ0,0°), (8) becomes
(17)V(i,0)=Pξ+V(i+1,0)+Pξ-V(i-1,0)+Pη+V(i,1).
The radial parts of the transition probabilities (Pξ+ and Pξ-) in (17) are the same as those in (14a) and (14b). The transition probability for the η coordinate component of (17) is
(18)Pη+=([1/(Δη)2][(1/(Δξ)2)+(1/(Δη)2)]).
Note also that Pξ++Pξ-+Pη+=1. For dispersed particles at point (ξ0,0°), they have probabilities Pξ+,Pξ-, and Pη+ of moving to points (ξ+Δξ,0°),(ξ-Δξ,0°), and (ξ,Δη), respectively. The particle walking pattern of the particles is as follows:
(19)(ξ0,0)⟹(ξ0+Δξ,0),if(0<U<Pξ+),(ξ0,0)⟹(ξ0-Δξ,0),if(Pξ+<U<PT0),(ξ0,0)⟹(ξ0,Δθ),if(PT0<U<1),
where
(20)PT0=Pξ++Pξ-.
The potential at a specific point (ξ0,η0) is to be determined. We define the transition probability Pn as the probability that a random walk starting at the point of interest (ξ0,η0) ends at the boundary node (ξn,ηn) with prescribed potential (Vb(n)).Large particles are dispatched at the free node (ξ0,η0). The application of Exodus method begins by setting particles P(ξi,ηj)=0 at all other nodes (both fixed and free), except at the free node (ξ0,η0), where we assume a large value N. By scanning the mesh as in finite difference analysis (FDM), the particles are dispatched at each free node to its neighboring nodes according to the random walk transition probabilities described above. Detailed description of this process is given in [10]. If N and Nn designate the total number of particles dispersed and the number of particles that have reached the boundary n, respectively, the probability that a random walk terminates on the boundary is
(21)Pn=NnN.
If there are M boundaries or fixed nodes (excluding the corner points), the potential at the specified node (ξ0,η0) is
(22)V(ξ0,η0)=∑n=1MPnVb(n)=∑n=1MNnVb(n)N.
Since there are just two boundaries (the inner and outer oblate shell surfaces) in this potential problem, (22) simplifies to
(23)V(ξ0,η0)=(N1Vb(1)+N2Vb(2))N,
where Vb(1) and Vb(2) are the prescribed potentials at the inner and outer oblate spheroidal shells, respectively.The equation for the exact solution of the potential computation in oblate spheroidal shells is [9]
(24)V=Vb(1)+(Vb(2)-Vb(1))*([gd(ξ)-gd(ξ1)]gd(ξ2)-gd(ξ1)),
where Vb(1) and Vb(2) represent the nonzero potential (Dirichlet boundary conditions) at the inner and outer conducting oblate spheroidal shells, respectively, as shown in Figure 3. Also, gd(·) denotes the Gudermannian function [6] and is represented as
(25)gd(ξ)=sin-1(tanh(ξ)).


## 3.1. Finite Difference Transformation of Oblate Spheroid Laplace’s Equation

The Laplacian equation in oblate spheroidal coordinate systems is(4)0=∇2V=(1sinh⁡2ξ+sin2η)*[(sech2ξtan2η+sec2ηtanh2ξ)∂2V∂ϕ2+∂2V∂ξ2+tanhξ∂V∂ξ+∂2V∂η2-tanη∂V∂η].
The term outside the square bracket may be ignored. Also the first term inside square bracket is ignored due to the rotational symmetry about the vertical (z) axis. Therefore, (4) reduces to
(5)∂2V∂ξ2+tanhξ∂V∂ξ+∂2V∂η2-tanη∂V∂η=0.
Equation (5) governs potential distribution in an axisymmetric oblate spheroid potential problem.Two oblate spheroidal shells made up of two constant conducting surfacesξ1 and ξ2 are shown in Figure 2. The two equipotential surfaces are maintained at 50 V and 100 V (Dirichlet boundary conditions), respectively. The choice of these Dirichlet boundary conditions is arbitrary. Any potential value can be assigned. Also, the value of the constant oblate spheroidal surfaces equipotential that constitutes the two conducting shells are arbitrarily chosen as ξ1=1.0 and ξ1=2.0, respectively.Figure 2
Two oblate spheroidal shells.The explicit finite difference transformation of (5) is
(6)V(i,j)=0.5[1[(1/(Δξ)2)+(1/(Δη)2)]]*[[1(Δη)2+tan(ηj)2Δη]V(i+1,j)[1(Δξ)2+tanh(iΔξ)2Δξ]+V(i-1,j)[1(Δξ)2-tanh(iΔξ)2Δξ]+V(i,j+1)[1(Δη)2-tan(ηj)2Δη]+V(i,j-1)[1(Δη)2+tan(ηj)2Δη]].
Figure 3 shows one-quarter of the constant oblate spheroidal surfaces. The figure exhibits symmetry with respect to the η coordinate. Therefore, two lines of symmetries will be encountered in this range of η. They are η=90° and η=0°. On these lines of symmetries, the condition ∂V/∂η=0 is imposed. This strategy eliminates the singularity causing term (tan(ηj)/2Δη) at the oblate spheroid poles as seen in (6).Figure 3
Oblate spheroidal surface path used in computing the charge enclosed.Consequently, the finite difference equations along the two lines of symmetries become as follows.Alongη=90°, j=jmax⁡,
(7)V(i,jmax⁡)=0.5[1[(1/(Δξ)2)+(1/(Δη)2)]]*[[1(Δη)2]V(i+1,jmax⁡)[1(Δξ)2+tanh⁡(iΔξ)2Δξ]+V(i-1,jmax⁡)[1(Δξ)2-tanh(iΔξ)2Δξ]+2V(i,jmax⁡-1)[1(Δη)2]].Alongη=0°, j=0,
(8)V(i,0)=0.5[1[(1/(Δξ)2)+(1/(Δη)2)]]*[[1(Δη)2]V(i+1,0)[1(Δξ)2+tanh⁡(iΔξ)2Δξ]+V(i-1,0)[1(Δξ)2-tanh⁡(iΔξ)2Δξ]+2V(i,1)[1(Δη)2]].


## 3.2. Transition Probabilities Determination

In a more compact form, (6) can be rewritten as
(9)V(i,j)=Pξ+V(i+1,j)+Pξ-V(i-1,j)+Pη+V(i,j+1)+Pη-V(i,j-1),
where(10a)Pξ-=0.5([(1/(Δξ)2)+(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(10b)Pξ-=0.5([(1/(Δξ)2)-(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(10c)Pη-=0.5([(1/(Δη)2)+(tan⁡(jΔη)/2Δη)][(1/(Δξ)2)+(1/(Δη)2)]),(10d)Pη-=0.5([(1/(Δη)2)-(tan⁡(jΔη)/2Δη)][(1/(Δξ)2)+(1/(Δη)2)]).Note that Pξ++Pξ-+Pη++Pη-=1. Equations (10a), (10b), (10c), and (10d) serve as the transition probabilities for the Exodus method used in this work. If N particles are dispersed at node (ξ0,η0), they have probabilities Pξ+,Pξ-,Pη+, andPη- of moving to points (ξ+Δξ,η),(ξ-Δξ,η),(ξ,η+Δη), and (ξ,η-Δη), respectively. The direction of movement of the particles is determined by generating a random number U,0<U<1. The particles are instructed to walk as follows:
(11)(ξ0,η0)⟹(ξ0+Δξ,η0),if(0<U<Pξ+),(ξ0,η0)⟹(ξ0-Δξ,η0),if(Pξ+<U<PT),(ξ0,η0)⟹(ξ0,η0+Δη),if(PT<U<PTT),(ξ0,η0)⟹(ξ0,η0-Δη),if(PTT<U<1),
where
(12)PT=Pξ++Pξ-,PTT=Pξ++Pξ-+Pη+.For the random walk computation along the line of symmetry (η=90°), (7) becomes
(13)V(i,jmax⁡)=Pξ+V(i+1,jmax⁡)+Pξ-V(i-1,jmax⁡)+Pη-V(i,jmax⁡-1),
where(14a)Pξ-=0.5([(1/(Δξ)2)+(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(14b)Pξ-=0.5([(1/(Δξ)2)-(tanh⁡(iΔξ)/2Δξ)][(1/(Δξ)2)+(1/(Δη)2)]),(14c)Pη-=([(1/(Δη)2)][(1/(Δξ)2)+(1/(Δη)2)]).Note that Pξ++Pξ-+Pη-=1 along this line of symmetry. If a random-walking particle is instantaneously at the point (ξ0,90°), it has probabilities Pξ+,Pξ-, and Pη- of moving to points (ξ+Δξ,90°),(ξ-Δξ,90°), and (ξ,90°-Δη), respectively. The particles are instructed to walk as follows:
(15)(ξ0,90°)⟹(ξ0+Δξ,90°),if(0<U<Pξ+),(ξ0,90°)⟹(ξ0-Δξ,90°),if(Pξ+<U<Pmax⁡),(ξ0,90°)⟹(ξ0,90°-Δη),if(Pmax⁡<U<1),
where
(16)Pmax⁡=Pξ++Pξ-.For the random walk computation along the line of symmetry(ξ0,0°), (8) becomes
(17)V(i,0)=Pξ+V(i+1,0)+Pξ-V(i-1,0)+Pη+V(i,1).
The radial parts of the transition probabilities (Pξ+ and Pξ-) in (17) are the same as those in (14a) and (14b). The transition probability for the η coordinate component of (17) is
(18)Pη+=([1/(Δη)2][(1/(Δξ)2)+(1/(Δη)2)]).
Note also that Pξ++Pξ-+Pη+=1. For dispersed particles at point (ξ0,0°), they have probabilities Pξ+,Pξ-, and Pη+ of moving to points (ξ+Δξ,0°),(ξ-Δξ,0°), and (ξ,Δη), respectively. The particle walking pattern of the particles is as follows:
(19)(ξ0,0)⟹(ξ0+Δξ,0),if(0<U<Pξ+),(ξ0,0)⟹(ξ0-Δξ,0),if(Pξ+<U<PT0),(ξ0,0)⟹(ξ0,Δθ),if(PT0<U<1),
where
(20)PT0=Pξ++Pξ-.
The potential at a specific point (ξ0,η0) is to be determined. We define the transition probability Pn as the probability that a random walk starting at the point of interest (ξ0,η0) ends at the boundary node (ξn,ηn) with prescribed potential (Vb(n)).Large particles are dispatched at the free node (ξ0,η0). The application of Exodus method begins by setting particles P(ξi,ηj)=0 at all other nodes (both fixed and free), except at the free node (ξ0,η0), where we assume a large value N. By scanning the mesh as in finite difference analysis (FDM), the particles are dispatched at each free node to its neighboring nodes according to the random walk transition probabilities described above. Detailed description of this process is given in [10]. If N and Nn designate the total number of particles dispersed and the number of particles that have reached the boundary n, respectively, the probability that a random walk terminates on the boundary is
(21)Pn=NnN.
If there are M boundaries or fixed nodes (excluding the corner points), the potential at the specified node (ξ0,η0) is
(22)V(ξ0,η0)=∑n=1MPnVb(n)=∑n=1MNnVb(n)N.
Since there are just two boundaries (the inner and outer oblate shell surfaces) in this potential problem, (22) simplifies to
(23)V(ξ0,η0)=(N1Vb(1)+N2Vb(2))N,
where Vb(1) and Vb(2) are the prescribed potentials at the inner and outer oblate spheroidal shells, respectively.The equation for the exact solution of the potential computation in oblate spheroidal shells is [9]
(24)V=Vb(1)+(Vb(2)-Vb(1))*([gd(ξ)-gd(ξ1)]gd(ξ2)-gd(ξ1)),
where Vb(1) and Vb(2) represent the nonzero potential (Dirichlet boundary conditions) at the inner and outer conducting oblate spheroidal shells, respectively, as shown in Figure 3. Also, gd(·) denotes the Gudermannian function [6] and is represented as
(25)gd(ξ)=sin-1(tanh(ξ)).


## 4. Results

Equations (10a)–(23) and (24)-(25) are used to implement the computation of potential distributions within two conducting oblate spheroidal shells numerically and analytically, respectively. The results obtained are as shown in Table 1. The results obtained from exact solution (analytical), the Exodus method, and the explicit finite difference solutions are compared in Table 1. Same step size was used for both the Exodus method and the FDM and this accounted for the closeness in the computed results obtained. The Exodus solution results were very close to the results obtained from exact solution because of the fact that though the Exodus method is a probabilistic method, its operation does not depend on random number generation which ultimately depends on the computation accuracy of the machine involved.Table 1
Comparing the Exodus solution with finite difference and exact solution.
Potential point
Exact potential (V)
Exodus methodpotential (V)
Finite difference (FDM)potential (V)
ξ
η
1.1
66°
57.1513
57.1542
57.1536
1.2
69°
63.7529
63.7577
63.7566
1.3
72°
69.8269
69.8328
69.8313
1.4
75°
75.3997
75.4060
75.4043
1.5
78°
80.5005
80.5067
80.5050
1.6
81°
85.1601
85.1656
85.1642
1.7
84°
89.4096
89.4141
89.4129
1.8
87°
93.2798
93.2830
93.2822
1.9
87°
96.8004
96.8021
96.8017


## 5. Conclusion

The use of the Exodus method to compute potential distribution inside two conducting oblate spheroidal shells maintained at two potentials has been implemented in this paper. The results obtained agreed with those obtained using finite difference (FDM) solution and the exact solution method. The Exodus method employed in this work can be said to be almost as accurate as the exact method when compared to the fixed random walk Monte Carlo method because the results obtained were as accurate as the exact method.



---

*Source: 126905-2014-03-17.xml*