VIEWS: 0 PAGES: 23 POSTED ON: 9/13/2012 Public Domain
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 coskz - t dt m p 1 3 (14) - 3 where k 3qVm 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