7PO~,,,,,:<vvv@.xv~ *, The Finite Difference Method for a Circular Cylinder. Analytical techniques can be used to solve some diffusion problems. However, the solutions obtained by analytical methods are not always suitable for computing numerical values of the solution in cases where the geometry is complex and the material properties are not uniform. Also, in many cases the analytical solutions do not apply (e.g. the boundary conditions will not permit a separation of variables) for solving the problem. In such cases, the analytical solution can be used to describe some of the fundamental characteristics of the solution being studied but the problem is solved again by numerical methods for the purpose of generating an accurate solution. With today's high-speed computers, numerical techniques are important for this type of problem. The current density diffusion equation can be formulated as a boundary and initial value problem as n02J = s J/t s.t. J(r, t) = 0 for t<0 J(0, t)/r = 0 J(r, 0) = dR(r) dR(r) = 0 for 0r 0+. At t=0, the analytical solution for the step current input predicts that all the current passes on the outer surface (skin) of the conductor. This situation is impossible to model numerically. What can be done is that we will look at the analytical solution at some time t0+Dt and start off our numerical procedure by using the values of J(r, t0+Dt) as an initial condition. This approach is messy and leads to computational difficulties in getting the solution started. For example, the outermost element of the conductor is assigned a value which probably would result in the discrepency that the integral of the current density through the cross section of the conductor might not add up to be the total current in the conductor and thus all the values of the current density would be scaled accordingly. We will not use this approach and will reformulate our problem in another way, working directly with the magnetic vector potential. On the basis of the mathematical model described by equations (2.30) and (2.31), we can develop, using finite difference or finite element methods, an algorithm for the numerical calculation of Az(r, q, t). As we have stated in equations (30) and (31) above, we want to solve the system : (4.1) n0 2 A - s A/t = - Js = - s C(t) (4.2) s C(t) - s A/t = JT where C(t) = - f gradient of the scaler electric potential defined over the surface of the conductor at time tn Js = -s f the source current density in the conductor -s A/t the eddy current density in the conductor JT the total current density in the conductor s the conductivity of the conductor n0 the reluctivity of free space, in our case a constant This formulation is discussed in Garg and Weiss [4.1]. If we integrate this above equation over the cross-section F of the conductor, we obtain our current constraint (4.3) F s C(t) dS - F s A/t dS = I(t) where I(t) is the given current in the conductor at time t. Because the conductor is uniform in the z-direction, we have that f is a constant over the surface of the conductor: (4.4) C(t) F s dS - F s A/t dS = I(t) Solving this equation (4.4) for C(t) and substituting into equation (4.1) yields F s A/t dS I(t) (4.5) n0 2 A - s A/t = - - F s dS F s dS The problem to which we will apply the finite difference methods describes an infinitely long cylindrical conductor having a radius R. This conductor will be placed in free space and carry a current I(t) with a coaxial return at infinity. Inside the conductor, the conductivity will be that of copper at 273 degrees Kelvin, and outside the conductor it will be that of free space, namely zero. In cylindrical coordinates, the Laplacian is (4.5b) 2A = 2A/r2 + 1/r A/r + 1/r2 2A/q2 where r is the radial distance. With radial symmetry (A/q = 0), this reduces to (4.5c) 2A = 2A/r2 + 1/r A/r . In our problem F is a disc of radius R. Now a grid for the cross-section is defined by concentric circles at rk = k (r), 0 1/2, this formulation is unstable. By using the values at the previous time in the adjacent nodes in modeling the diffusion, the explicit scheme introduces characteristics into the behavior in that information at the jth node takes a time kDt to arrive at j+kth node. Implementation of an algorithm in which the time step is too large in relation to the space nodes results in instability. For the differential form of the diffusion equation and for the physical problem, the diffusion effects still occur at this same rate but the transfer of information is instantaneous (is that of the speed of light). The presence of our source term and the coupling to the current constraint do not permit us to solve the system by merely taking a combination of values of surrounding nodes from the previous time step. The Crank-Nicolson method has been the preferred method for the solution of the parabolic partial differential equations in one space dimension. Its advantage is higher precision for a given step size. The error term is of the form E2 = K2 Dr2 + K'2 2Dr + C2 Dt2, where C2 = (1/2!) 3A(k, n+l2 Dt)/3t, 0 < l2 < 1, K2 = (1/3!) 3A(k+m2 Dr, n)/r3, 0 < m2 < 1, K'2 = (1/2!) 2A(k+m'2 Dr, n)/r2, -1 < m'2 < 1. It also retains the parabolic character of the original equation because the solution at the point xj at time tn+1 is influenced not only by all the earlier values, but also by the solution at all other points at time tn+1 as well. The disadvantage is the need to solve a tri-diagonal matrix system at each time step. However, this is a straight forward procedure. The Larson totally implicit method has also been used to solve parabolic systems. The truncation error is of the form E3 = K3 Dr2 + K'3 2Dr + C3 Dt, where C3 = (1/2!) 2A(k, n+l3 Dt)/2t, 0 < l3 < 1, K3 = (1/3!) 3A(k+m3 Dr, n)/r3, 0 < m3 < 1, K'3 = (1/2!) 2A(k+m'3 Dr, n)/r2, -1 < m'3 < 1 . This scheme is of the same order as is the explicit scheme but one order lower in time than is the Crank Nicolson formulation at the same computational cost. Usually, C3 << C1. However, it should be pointed out that it is not only the magnitude of the order which is important, but also the constant Ki. One should be aware of the possibility that for a certain range of the parameters, t and r, the lower order method may be more accurate because the constant C3 may be less than C2 Dt. In essence, this may be trivial since it is very costly to calculate the constants Ci and Ki because of their dependance on so many parameters which occur in the geometry of each problem. When modelling our problem using relatively large time steps (that is for q 10), we avoid using the Crank-Nicolson scheme because of the oscillatory behavior exhibited even though it is stable. Zienkiewicz and Morgan [4.8] illustrate this. However, the most important reason we use the totally implicit method is for the advantages that occur in the special form of the forcing term. More will be stated about this in the section on finite elements. Now we would like to write our difference equation using a weighted average approximation. Considering the different cases - for q = 0, this is a Euler explicit forward difference formulation, for q=1/2, this is a Crank-Nicolson implicit central difference formulation, and for q=1, this is a fully implicit backward difference formulation. The differential equation (4.5) becomes the difference equation  For k = 0, 1, ... M ; n = 0, 1, ... N. where a(k, n) is the magnetic vector potential at the kth node in space and nth node in time, Fk the area of the region center at the node k, s(k, n, T) the conductivity of the region centered at the node at a temperature T, I(tn) the current magnitude at time tn t time step, r space increment, q variable weight parameter, 0q1 k=0 node corresponding to the center of the conductor, k=M node corresponding to the outer-most region of the conductor, n=0 node corresponding to the initial time, n=N node corresponding to the final time, . This equation as it stands has as unknowns s(k, n+1, T) and a(k, n+1). The coefficients s(k, n+1, T) can be estimated using s(k, n, T) and extrapolating the anticipated temperature rise. However, since this effect is negligible over one time step, we set s(k, n+1, Tn) = s(k, n, Tn) and update the value at the end of the time step. We make one more assumption - the reluctivity, n, is a constant since the material is non-ferromagnetic. If this were not so, we would have to specify its behavior. The boundary conditions for this problem are also straight forward. At the Neumann boundary we have A/r = 0. At r=0, we should use the condition that a(-1, n) = a(1, n). Hence in our finite differences, at k = 0, we must change the coding so that we do not address variables out of bounds. Also, we must change the form of the difference equation (4.6b) because of an apparent singularity of 1/r at the origin. This problem exists only because of the choice of the cylindrical coordinate system. If we take the limit as r tends to zero of equation (4.5c), we get by application of l'Hospital`s rule (4.11) 2A (r=0) = limit 2A/r2 + 1/r A/r = 2 2A/r2 . r -> 0 At the center of our grid, we change the computation molecule to use the above information. Notice that since s=0 in free space, our differential equation (4.5) reduces to Laplace's equation in this region. At the Dirichlet boundary, we have A equal to some prescribed value. We will not be trying to couple our finite difference method to the solution which holds outside the conductor, namely A(r, t) = -moI(t)/2p ln(r) + K, r > R (see equation 2.65) where K is an arbitrary constant. We choose to end our finite difference mesh at r=R, where R is the radius of the conductor. A(r, t) is a potential and, thus, it can be specified up to a constant. Since we have radial symmetry, we might as well make it zero at the end our mesh. Moreover, we are interested only in the change in A in time so that its absolute value is irrelevant. Through the use of Neumann's Fourier stability analysis, the above formulation is shown to be unconditionally stable when the source term is neglected, see Kreiss[4.9] and Khalig and Twizell [4.10]. We tried to show that this is also true when the source term is included, however the algebra becomes intractable and it is no longer possible to group the terms as was done for the simpler problem. Conclusion A method of calculating the eddy currents in a cylindrical conductor by a suddenly applied current to the system is given. It is shown that the analysis of the eddy currents in a cylindrical transients current carrying conductor leads to the well known parabolic differential equation with given initial and boundary conditions. However, its application to eddy current distribution in the above cylinder after applying a current ramp input function as well as the study of the error in the presence of the source current density term were not previously done. We illustrate the numerical techniques described in this section by presenting the results of numerical calculations where these techniques have been used to simulate the magnetic vector potential and the current density distribution modeled by equation (4.1) and (4.3) for several different times. The plots of the magnetic vector potential a(r, t) and the current density vector J(r, t) and da(r, t)/dt are shown in figures 4.3, 4.4, and 4.5 . They indicate the changes of the above field quantities in an infinitely long circular cylindrical conductor as a function of its radius at four different times t=0.5 msec, 1.0 msec., 1.5 msec., and 2.0 msec. The magnetic vector potential is given in units of Volt seconds / meter. The current density is given in Amps/meter2. Figure 4.3 Current density distribution vs. radial distance at four different times. Figure 4.4 Magnetic vector potential distribution vs. radial distance at four different times. Figure 4.5 The change in the magnetic vector potential distribution in time vs. radial distance at four different times. Figure 4.6. Finite difference results at time = 1.0 msec, 2.0 msec, ..., 10.0 msec, using Dt = 0.25 msec. and constant conductivity, s = 5.0 x 108. Figure 4.7 Close-up of previous diagram. References: [4.1] Garg, V. K. and Weiss, J., "Finite element solution of Transient Eddy-Current Problems in Multiply-Excited Magnetic Systems," IEEE Transactions on Magnetics, September 1986, Vol. MAG-22, No. 5, pp. 1257-1259. [4.2] Smith, G. D., Numerical Solution of Partial Differential Equations, Finite Difference Methods, Oxford University Press, 1978. [4.3] Lapidus, L. and Pinder, G. F., Numerical Solution of Partial Differential Equations in Science and Engineering, John Wiley and Sons, New York, (1982) [4.4] Vermuri, V. and Karplus, W. J., Digital Computer Treatment of Partial Differential Equations. Prentice - Hall, Inc., Englewood Cliffs, New Jersey (1981). [4.5] Vichnevetsky, R. , Computer Methods for Partial Differential Equations, Vol. I, Prentice - Hall, Inc, (1981) [4.6] Twizell, E. H., Computational Methods for Partial Differential Equations. Ellis Horwood Ltd. - John Wiley and Sons, New York (1984). [4.7] Evans, A., Lecture Notes For Numerical Solutions to Differential Equations, McGill University, Spring 1985. [4.8] Zienkiewicz and Morgan, Finite Elements and Approximations, John Wiley and Sons, (1983) [4.9] Kreiss, H. O. , Numerical Methods for Solving Time Dependent Problems for Partial Differential Equations. Les Presses de L'Universite de Montreal, C. P. 6128, Succ. "A", Montreal, (1978) [4.10] Khalig, A. Q. M. and Twizell, E. H., "L0-Stable Splitting Methods for the Simple Heat Equation in Two Space Dimensions with Homogeneous Boundary Conditions," SIAM Journal Numerical Analysis. Vol. 23, No. 3 June 1986, pp.473-484. th 4 FD page   OfOfVm Vm A` +:n ( n+1+7k )k+1 (Ok-1 8A8`" A6" e6" 6" 6" 6" 6" 6" =6" _6(O=0)%1)M (O4M-1TXT6;XT6;XT6@;DXT6d;hXT6;XTXTXTXTXTdhXT@DXT6;XT6;XT6<;@XT6^;bXT;?XT]aXfhrXhrHoy oy"$"6H"d$"6@H"d$"d$$"6<"H"6<"6<$ܠ"H +Ugq = 0)qq=1/2 )pq=1(:n ( n+1+Y6k )k+1 (N<k-1)k )k+1 (Nk-1)k )k+1 (Nk-1 )v 7)7v" A6" e6" 6" 6" 6" 6" 6" =6" _6 4 B 84 C8wddkz k + (4.9)zV({"a(k+1, n) - 2 a(k, n) + a(k-1, n))z+$ +z(.rz(.Dz ((2" f"z ( Tn ".g" ""-" " f k; NСz( 4(1-q)zD(Ua(k+1, n) - a(k-1, n)z+ 2 r rz)Dz (3}kzk (]m"a(k+1, n+1) - 2 a(k, n+1) + a(k-)1, n+1)z+* +zD(]Na(k+1, n+1) - a(k-1, n+1)z(srz(sDz (m2"Qf"z (eTn "sg"O""r"O"Qf bkb; bNbСz(fKqz(p 2 r rz)Dz (vkz1 (<a(k, n+1) - a(k, n)z+3tz(eD">z  (Yz + j=0z(IMz+jz (Wsz(Wqz )Q1-qz(V' ( )z +`)z (VIsz(V+z (V(z + (j, n, T)z(VQ (j, n+1, T)z(3 (k, n, T)z+#FaZoNu`^z1y` ffp fm,$CVx+Xz1^^ /`0D`1^^  `1^ w?ABCEde  !=!>!?!@!A!B!!!!!!!!$$$%$&$'$H$I$J$L$M$O$R$S$T%%% % &n&o&p&q0P&q&&&&&&&&&&&($(%()(+(,(-(.(1(2(3(5(6(7(;(=(>(?(@(Q(R(b(h(i(j(l(m(n(q(r(t(x(y(z(~(((((((((((((((((((((((((((((((((( @SLO(((((((((((((((((()))***-'-(-,---.-/-0-1-6-7-9-:-?-@-A-B-C-D-L-M-N-Y-Z-a-b-c-e-i-j-q-r-s-w-{-|--------------------------.>.@.IW.I.L../////////////////////////////00000000000000 0"0%0&0*0+0,0.06070=0>0?0H0I0M0N0O0Q0T0U0Y0Z0[0\0^0f0g0n0o0p0q1#1$1)1*112N2O2PV2P2b2d2e22223s3t45r5s55555566 666667 7 777N7O77777777889!9"9O9P9s9t99:::::::w:x;W;Y;Z;\=P=Q=b=c=d=f=z={=|============= @VS=======>$>%?M?N?O?P?Q?T?U?WBBGGHHHIIKILIIIIJ<JBJDJPJJK?KKL&LvLM MGMMNNHNNNO;OOP5PUP~PPPPPPPj^@X@?6@)k   K L p q x Ϳyk]kkOkA ! h* ! h* ! h* ! h* ! h* !@* !@* !@* !4@* !@* !@* !@* !@*! B  %-4,HI55ȸxhxxxxhXhxxJ !@*!h*!h*!h*!  *!  *!  *!  * ! h* ! h* ! h* ! h*56hi()GHI !#[ 4YZFG̾wiiii  *  * p$T*p!$T*p$T* @* @* !@* !@* !@*)[aE#h  X !E!!#z$$>$$%#&l&&&&&'-'.(d((,;-z-/K00F466666666677M777788Q88888ŵ ` * ` *@* @*  *A8:G:;7=S=T===> ?CABBBBCDEGGGGGGGGGGGGGGGGGGGHHHIHKHLHMHNHOHPHQHRHSHTHUHVHWHXHYHZH[H\HHHHHHHHHHHHHHHHHHHIKIL @* * * *NILIMINIOIPIQIRISITIUIVIWIXIYIZI[I\I]I^IIIIIJJJJJJJJJJ J J J J JJJJJ<J=J>J?J@JAJBJCJDJPJQK)K*KKLNLOLLMgMhMMNkNlNNOOP}P~PP   # * *JPPPPPP! `p # Regular OP -O234,P5 J&q(.I2P=j^6789:;<=>? 58ILPP@ABCDEFHP(HPP  '=/pJ2@P -:ImageWriter (K(