Numerical Simulation of a TWDEC for an ARTEMIS Reactor by G3Li2Hfi

VIEWS: 0 PAGES: 23

									Numerical Simulation of a TWDEC for an ARTEMIS Reactor
                                              Kevin N. Jarboe
                                                 12/16/99
                                                 ECE 272


Abstract
Numerical Simulations of proton behavior in the TWDEC of an ARTEMIS fusion reactor were carried out
in the MATLAB programming language. Studies show that proton velocity decreases in the modulator.
Protons attain a maximum number density approximately 105 m from the modulator exit. Velocity
decreases in the decelerator as protons lose kinetic energy. A.C. power generation is similar to a simple
amplifier and is possible, in principle.




1. Introduction
The search for an efficient, non-polluting means of energy production is an ongoing struggle. Current
methods of energy production include burning fossil fuels and forced nuclear fission. Both processes are
inefficient and wasteful. Furthermore, neither energy source is expected to last for more than a century.
Thus, a new source of energy must be developed in the near future. Controlled thermonuclear fusion shows
promise as a new, desirable energy source. Work must be done to realize this new energy source.


2. Total Scheme

The ARTEMIS fusion reactor, developed by Hiromu Momota et al., burns a D-3He plasma in a Field

Reversed Configuration (FRC):

                           D + 3He  p +  + 18.3 MeV                                      (1)

         ARTEMIS, shown in Fig. 1, burns D-3He plasma in the formation and burning chambers. The

Direct Energy Converter (DEC) segments of ARTEMIS convert energy from fusion products into useful

A.C. power [1].




                                                     1
         ARTEMIS uses two types of DECs: the CUSPDEC to separate electrons and fusion protons, and

a cylindrical Travelling Wave Direct Energy Converter (TWDEC) [2].



                                                                                  Superconducting
                                              Cryopump                                 coils




                                                                                         TWDEC
       CUSPDEC
                            Fig. 2 Schematic of ARTEMIS Direct Energy Convertors



         The CUSPDEC is shown on the left of the diagram. The TWDEC, on the right, is the focus of this

study. Superconducting coils form a magnetic field in the DEC to guide protons. Cryopumps maintain the

necessary vacuum within the DEC. Cylindrical grids, the vertical lines in Fig. 2, carry potential energy to

and from the TWDEC.



         The TWDEC has other parameters [3]:

Total length of TWDEC             30 m
Wave length of modulator ()      2 m
Number of modulator grids         5
Number of decelerator grids       9 to 13
Length of modulator               2 m
Length of decelerator             2 2.8 m
                                              .

Table 1 TWDEC parameters




                                                      2
        Fig. 3 is a schematic diagram of one-quarter of one grid. Each grid is a coaxial set of conducting,

water-cooled pipes [4].




                                                     3
3. Modulator

         In this one-dimensional study, it is sufficient to visualize a "stream" of 15 MeV protons, evenly

spaced over 2, approaching the modulator in the TWDEC.

         Grids in the modulator excite a standing potential wave:

                                                      
                           V ( z, t )  Vm sin(kz) cos( t )
                                         2                                                          (2)
                            where k 
                                         Lm

Vm is the voltage magnitude, Lm is the modulator length, and  is the angular frequency of the potential

wave.



                                                                                                   Modulator
                                                                                                   entrance



 ith Proton             4th Proton          3rd Proton          2nd Proton            1st Proton



                                                2

        Fig. 4    15 MeV protons approaching the modulator


         Upon entering the modulator, protons accelerate according to the momentum equation for protons:

                            d 2z    q V ( z, t )
                                                                                                   (3)
                            dt 2
                                    mp  z

                            dv    qV
                                                
                                 m cos(kz) cos( t )                                                (4)
                            dt    mp

where q is electron charge and mp is the proton rest mass.


         Velocity changes due to the potential wave is determined at any known z(t), t. Because z(t) is

unknown in the modulator, the Leap Frog method is used to provide a close approximation. The Leap Frog

method assumes that the velocity differential over a small time interval, dt, is approximately constant:




                                                         4
                                       qVm
                             dv                     
                                           cos(kz) cos( t )  dt
                                       mp

         (5)

To simplify calculations, let v0 = z0 = t0 = 0.




         To initialize the process, the velocity change over a small interval, dt/2, is determined and t is

incremented accordingly:

                                                     qVm                        dt
                               dv( z 0 , t 0 )                      
                                                         cos(kz0 ) cos( t 0 )                         (6)
                                                     mp                         2

                                     dt 
                             z  t    z 0  v0  dv( z 0 , t 0 )  
                                                                         dt
                                     2                                 2
                                                                                                       (7)
                                   dt
                             t
                                   2
These coordinates specify an approximate value of velocity at time dt:

          v(t  dt )  v0  dv
                        dt    dt   qV               dt     dt          (8)
          where dv z t  , t     m cos k  z  t    cos       dt
                                                                       
                        2     2    mp               2    2 

Given these initial parameters, one begins to construct a function z(t). Using the v(t) from (8), calculations

give the position at time 3dt/2 and t is incremented:

                                   3dt 
                             z t        z  v z t  dt , t  dt   dt
                                    2 
                                                                                                       (9)
                                          3dt
                             t  t  dt 
                                           2
Velocity at time 2dt is approximated as in Eq. (8):

          v(t  2dt )  v0  dv
                          3dt       3dt     qV                3dt    3dt            (10)
          where dv z  t 
                               , t         m cos k  z  t 
                                                                        cos  
                                                                                        dt
                           2         2      mp                 2            2 

The new v(t) gives the position at 5dt/2,




                                                            5
               5dt                 3dt       3dt 
         z t        z  v z  t 
                                          , t        dt
                2                   2         2                                             (11)
                      5dt
         t  t  dt 
                       2
t is incremented and velocity is approximated at time 3dt,

v(t  3dt )  v0  dv
               5dt       5dt     qV                5dt    5dt 
where dv z t 
                    , t         m cos k  z  t 
                                                             cos       dt
                2         2      mp                 2  
                                                                         2 

         (12)

and so on.



The Leap Frog method is summarized pictorially:

Vo,z0    z(dt/2)           v(dt)             z(3dt/2)          v(2dt)            z(5dt/2)          v(3dt)




t=0      dt/2              dt                3dt/2             2dt               5dt/2             3dt


  Fig. 5 Illustration of the Leap Frog Method


As shown in Fig. 5, the leap frog method utilizes the relationship between dv and z to slowly construct an

approximation of v(t) and z(t).




                                                        6
         The MATLAB program listed in Appendix A-1 determines the velocity and position of N protons

passing through the modulator. A WHILE loop allows the program to consider each proton separately, and

the cumulative results are shown on one graph.




   Fig. 6   Velocity vs. position for a "stream" of 20 protons in a modulator


Fig. 6 shows that proton velocity follows a general trend in the modulator. No proton exhibits positive

acceleration in the modulator, and all protons have constant velocity for a short time.




                                                      7
          The modulator slows every proton that passes through, with maximum deceleration when i = 1.

As i in Fig. 4 gets large, proton velocity at the modulator exit increases slightly as shown in Fig. 7:




Fig. 7 Velocity of i-th proton after modulation



          Length 2, shown in Fig. 4, constrains t in the standing wave of Eq. (2) to 0 < t < (2/v0). Thus,

proton deceleration approaches an upper bound as i increases. Simulations for large N accurately define

this limit.




                                                       8
         After exiting the modulator, protons travel with constant velocity and "bunch" together. Fig. 8,

below, illustrates the "bunching" of two proton ions:



  Proton 2                         Proton 1                          Proton 2         Proton 1
  Velocity v2                      Velocity v1                       Velocity v2      Velocity v1



                          t=0                                                      t = tb


    Fig. 8   Model of protons 1 and 2 at time t after modulation

         The 1st proton to exit, p1, has velocity, v1, and the 2nd proton, p2, has velocity, v2, v2 > v1. At some

time tb, p2 will have nearly the same position as p1.

         Such "bunching" is measured in terms of number density, N/dz, where N = (# of protons) and dz =

(zmax-zmin). Proton position is computed for some t after modulation. Phase space plots and number density

plots are generated as shown in Fig. 9:




      Fig. 9 Number density and phase space plots of protons at time, t, after modulation




                                                        9
For the case of 20 protons, maximum bunching occurs approximately 105 m from the modulator, as shown

in Fig. 10:




     Fig. 10 Number density vs. position for modulated protons




Results of the modulator simulation show that maximum "bunching" occurs approximately 105 m from the
modulator exit and is relatively high at nearly 350 protons per meter.




                                                  10
 4. Decelerator

 Modulated protons enter the decelerator, spaced z meters from the modulator. Grids in the decelerator

 excite a travelling potential wave:

                    V ( z , t )  Vm cos(kz  t   )                                        (13)

 where k is the wave number and  is the phase shift.

          Protons "bunches" passing through the decelerator get "trapped" in this potential wave. Protons

decelerate and induce a potential on decelerator grids. A.C. power is generated as this potential builds up.

          Appendix A-2 lists a MATLAB program that constructs a velocity vs. position function for a

single proton in a decelerator. The Leap Frog method described earlier is used again in the decelerator, but

with a different equation of motion.


 dv  qVm k                    
          coskz - t   
 dt  m p
     
                               
                               
                                                         1
                                                 3                                          (14)
                                                 
                       - 3                      
 where k                                        
             3qVm
                  sin( ) * (z - z 0 ) - v 0  
                                                3
            m                                    
               p                                 
 Note that this decelerator simulation uses a 'Renormalization' technique to avoid round off errors.
 .

 If t, v, and  are units of real time, velocity, and angular frequency, then T, V, and  are scaled units of
 time, velocity, and angular frequency. In this simulation,


 t  10 9 T  T
       V      V
 v  9                                                                                      (15)
     10       
              
   9 
      10       

      V 
     d 
       
       2
 dv          1 dV
 dt d T   dT
                                                                                              (16)




                                                         11
Results of this 1 proton simulation are shown in Fig. 11:




 Fig. 11 Velocity vs. position for a proton in the TWDEC decelerator

After an initial rise in velocity, the proton velocity decreases in the decelerator as expected. A large
applied voltage magnitude may account for the increase in velocity seen in the latter portion of the
decelerator.




                                                      12
     5. A.C. Power Generation


     Due to time constraints, no A.C. power generation study was conducted. However, a brief description of
     the concepts follows.
               To illustrate the TWDEC operation, first consider an oscillator circuit, as shown in Fig. 12.




Input
from thermal
fluctuations




                                                        Output

                     Fig. 12   RLC oscillator circuit
               Input electrons from thermal fluctuations excite a triode, which drives a connected RLC circuit.

     Signals from the triode provide feedback and strengthen the output signal, which drives a load. Several

     circuits in parallel yield a large output power.

               The TWDEC circuitry is similar to the oscillator of Fig. 12. Schematics for one segment of the

     TWDEC circuitry are depicted in Fig. 13 below:


                   Modulator                      Decelerator
                     grid                            grid




                                 RL




               Fig. 13 Schematic of circuitry connected to TWDEC



          Protons exiting the modulator provide input to the "amplifier," a decelerator grid. Stray capacitance
     develops between a passing proton "bunch" and a decelerator grid, allowing electric potential to develop



                                                             13
between the grid and ground. This potential build-up causes high frequency (~30 MHz) power to oscillate
in a connected RLC circuit. The RLC circuit in the TWDEC drives an electrical load, as does the oscillator
in Fig. 3. However, a transformer in the TWDEC is typically coupled to a rectifier or an AC-AC power
converter and an electrical load. A fraction of the output power provides feedback to the modulator of the
TWDEC and the feedback strengthens certain signal components. As with the oscillator, several
decelerator grids connected in parallel give a proportionately higher power output.




                                                     14
6. Concluding Remarks


         The preceding studies examine the behavior of proton ions in a TWDEC. Studies of the
modulator were more extensive than those of the decelerator. The proton "bunching" produced by the
modulator is essential for an efficient decelerator. Simulations show that proton "bunches" can be
produced. Thus, theoretical studies show that, in principle, direct energy conversion is possible.




                                                     15
Appendix A-1 Listing of Modulator Simulation Code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     Distribution of modulated protons in the modulator
%     segment of the TWDEC in the ARTEMIS d-3He fusion
%     reactor
%
%     Kevin Jarboe
%     12/09/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

     f        =   8.54E+6;           % frequency
     q        =   1.6E-19;           % electron charge
     m_p      =   1.67E-27;          % proton rest mass
     pi       =   3.141592;
     Vm       =   0.42E+6;           % max potential
     qMVm     =   q*Vm/m_p;          % 4.02E+13;
     z0       =   0.;                % initial position
     t0       =   0.;                % initial time
     omega    =   2.*pi*f;           % 2*pi*f
     dt       =   1E-09;       % incremental time change
     v0       =   5.4E+7;            % initial velocity
     np       =   20;                % Number of protons
     l_m      =   2*pi;        % length of modulator
     k        =   2*pi/l_m;          % wave number in modulator
     lambda   =   1/(2*pi*k);        % wavelength
     pmod     =   150;               % length of post-modulator region
     step     =   10;
     z_s      =   (lambda)/(np);     % distance between any 2 protons
                                     % before modulation
     t_s      =   z_s/v0;            % time interval between any 2
                                     % protons before modulation

     i       =    1;
     plt     =    [0   0];
     j       =    0;
     mod     =    [0   0 0 0 0];
     v_f     =    [0   0 0 0];
     num_dens=    [0   0];
     phase_0 =    [0   0 0 0];
     phase   =    [0   0];
     row     =    0;

figure;
while i<=np
%       disp([' Proton # ',int2str(i)]);
%            Initial vel, z, t
             vel = v0;
             z = z0;
             t = t0+(i-1)*(t_s);
             V = Vm*(cos(omega*t))*(sin(k*z));
             dv = (-1*qMVm*k)*cos(omega*t)*cos(k*z)*dt/2;
             vel = vel + dv;
             z = z + vel*(dt/2);
             vel = v0;


                                    16
           t = t + dt/2;

           while z <= l_m
%          Write results to array, mod
           row = row+1;
           j=j+1;
           mod(row,:)= [t z V dv vel];
           plt(j,:)=[z vel];
           dv = (-1.*qMVm*k)*cos(omega*t)*cos(k*z)*dt;
           V = Vm*(cos(omega*t))*(sin(k*z));
           z = z + vel*dt;
           vel = vel + dv;
           t=t+dt;
           end

%          Approximate velocity at end of modulator
           dz_f = l_m-z;
           delta_vel = mod(row,5)-mod(row-1,5);
           delta_z = mod(row,2)-mod(row-1,2);
           row = row + 1;
           j= j+1;
           dv = dz_f*(delta_vel/delta_z);
           vel = vel + dv;
           z = l_m;
           t = t + dz_f/vel;
           V = Vm*(cos(omega*t))*(sin(k*z));
           mod(row,:)= [t z V dv vel];
           plt(j,:)=[z vel];
%          Record final velocity of i-th proton
           v_f(i,:) = [i t z vel];
           i = i+1;

%     Plot proton motion in a modulator
      plot(plt(:,1), plt(:,2),'.', plt(:,1), plt(:,2),'-');
      hold on;
      j=0;
      plt=[0 0];
end

% Label plot
xlabel('z [m]');
ylabel('v [m/s]');
title('Proton motion in a modulator');
orient landscape;
grid on;

% Plot velocity of i-th proton after modulation
figure;
plot(v_f(:,1),v_f(:,4),'.');
orient landscape;
grid;
xlabel('i-th proton');
ylabel('v_f [m/s]');
title('Velocity of i-th proton after modulation');




                                   17
%      Plot vel. vs. position for i-th proton in post-modulator region
z=0;
i=1;
t_max=v_f(1,2);

while i <= np

      if t_max < v_f(i,2)
             t_max = v_f(i,2);
      end
      i=i+1;
end


% Find z0 for each proton after all have exited modulator
i=1;
while i <= np
      z = (t_max-v_f(i,2))*v_f(i,4);
      t = v_f(i,2)-t_max;
      phase_0(i,:)=[i t z v_f(i,4)];
      i=i+1;
end


% Generate phase space plots every pmod/(step*v_f(1,4))
i=0;
t=0;


      v_max = v_f(np,4)
      v_min = v_f(1,4)

      j=1;
      while j <= np
            if v_max < v_f(j,4)
                  v_max = v_f(j,4)
            end

            if v_min > v_f(j,4)
                   v_min = v_f(j,4)
            end
            j=j+1;
      end

figure;
dt_pm=pmod/(step*v_f(1,4));
v = v_f(:,4);
while i < step
      z = phase_0(:,3)+v_f(:,4)*t;
        j=1;
      z_min = t*v_max;
      z_max = t*v_min;
        while j <= np
                if z_max < z(j,1)
                        z_max = z(j,1)
                end


                                      18
                     if z_min > z(j,1)
                             z_min = z(j,1)
                     end
                     j=j+1;
          end

        a = (np*i+1);
        b = (np*i+np);
        z_i = z_min;
        z_f = z_max;
        dz = z_f-z_i;
        n = np/dz;
        num_dens(i+1,:) = [(z_i+z_f)/2 n]

%     Graphical output
         figure;
         plt=[0 0];
         plt(1,:)= [z_i n];
         plt(2,:)= [z_f n];
         subplot(2,1,1),plot(plt(:,1),plt(:,2),'.',plt(:,1),plt(:,2),'-');
         grid on;
         title('Number density in post-modulator region');
         ylabel('N/dz [1/m]');
         xlabel('z [m]');

        hold on;
        subplot(2,1,2), plot(z,v,'.');
        title('Phase space plot in post-modulator region');
        ylabel('v [m/s]');
        xlabel('z [m]');
        hold on;
        orient landscape;
        grid on;
        phase(a:b,:)=[z v_f(:,4)];
        i=i+1;

        t=t+dt_pm;
end

% Graphical output
figure;
plot(num_dens(:,1),num_dens(:,2),'.');
grid on;
zoom;
orient landscape;
ylabel('N/dz [1/m]');
xlabel('z [m]');
title('Number density in post-modulator region');


figure;
plot(phase(:,1),phase(:,2),'.');
grid on;
zoom;
orient landscape;
ylabel('v [m/s]');


                                        19
xlabel('z [m]');
title('Phase space plot in post-modulator region');

figure;plot(phase_0(:,3),phase_0(:,4),'.');
grid;
orient landscape;
xlabel('z [m]');
ylabel('v [m/s]');
title('Phase space plot of protons at modulator exit (t=0)');




                                   20
Appendix A-2 Listing of Decelerator Simulation Code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     Motion of a proton in the decelerator segment of the TWDEC
%     in the ARTEMIS d-3He fusion reactor and uses scaling factor,
%     alpha = 1E-9
%
%     K. Jarboe
%     11/30/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%     declare variables
%        f  = frequency
%     phi   = phase angle
%     q     = electron charge
%     m_p   = proton rest mass
%     omega = 2*pi*f
%     phi   = phase angle of traveling wave
%     x0    = initial position
%     V0    = scaled initial velocity
%     dT    = scaled incremental time change
%     pi    = pi
%     xf    = final position
%     T0    = scaled initial time
%     Tf    = scaled final time
%     i     = counter used in Do loop...equivalent to t/dt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%


     f     = 8.54E-3
     q     = 1.6E-19
     m_p   = 1.67E-27
     pi    = 3.1416
     Vm    = 0.42E+6
     x     = 0.027
     x0    = 0.
     xf    = 17.59                        % 2.8*2pi
  qMVm     = q*Vm/m_p%4.02E+13            % q*Vm/m_p
     T     = 0.
     psi   = -3.*pi/4.
     omg   = 0.054                        % 2pi*f
     dT    = 1.
     V0    = 0.054
     V     = V0
     dV    = 0.
     i     = 1
     x     = V0*dT/2;




                                   21
k = 1.E+9*(-1*omg^3/((3E+9*omg*qMVm)*(sin(psi))*(x-x0)-
(V0)^3))^(1./3.);
      dV= (-1.0E-18)*((qMVm*k)*cos(k*x-omg*(T)+psi))*dT/2;
      x = x + (V+dV)*dT/2;
      T = T + dT/2

while (x <= xf)
       decel(i,:)=[T x k dV V];
k= 1.E+9*(-1*omg^3/((3E+9*omg*qMVm)*(sin(psi))*(x+(V*dT/2)-
    (V0)^3))^(1./3.);
       dV= (-1.0E-18)*((qMVm*k)*cos(k*x-omg*(T)+psi))*dT;
       V = V + dV;
       x = x + V*dT;
       T = T+dT;
       i = i+1;
end


figure;plot(decel(:,2),decel(:,5),'-',decel(:,2),decel(:,5),'.');
grid;orient landscape;
title('Proton in decelerator, phi = -3pi/4')
xlabel('z (m)')
ylabel('v (m/s)')



diary decel
diary on
disp('t(ns)   z(m)   k   dV/10^7(m/s)   V/10^7(m/s)')
out=[decel(:,1)*1E+09 decel(:,2) decel(:,3) decel(:,4)/1E+07
decel(:,5)/1E+07]
diary off




                                   22
References
1
    Momota et al., "Conceptual Design of the D-3He Reactor ARTEMIS," Trans. on Fusion Technology, 21,
          (1992), 2309.

2
    Momota et al., "A Traveling Wave Direct Energy Converter for D-3He Fueled Fusion Reactor," Trans. on
          Fusion Technology, 35, (1998), 60.
3
    Ishikawa et al., "Basic Numerical Simulation of TWDEC for a D-3He Fusion Reactor," Fusion
           Engineering and Design, 41, (1998), 543.
4
    Momota et al., "Loss Estimation of TWDEC for a D-3He FRC Fusion Reactor," Proc. of the US-Japan
          Workshop and the Satellite Meeting of ITC-9, NIFS-PROC-41, (Apr. 1999), 73.




                                                     23

								
To top