大気圏突入時の熱流束を知りたいー4. Sutton-Gravesの熱伝達モデルー
前のblogに書いた大気密度のモデルと宇宙船の運動モデルに加えて、Sutton-Gravesの熱伝達モデルを使うと大気から宇宙船への熱流束を評価することができる。
J.FayとF. Riddellが考案（Journal of Aeronautical Sciences, vol. 25, no. 2, 1958）し、K. SuttonとR. A. Gravesが簡易化（NASA TR-376, 1971）した便利な式があって、大気から宇宙船への熱流束は次の式で表せる。
これまでblogで書いてきたことを総動員すると、Apollo 6に対する各高度の熱流束は次の図のように計算できる。高度が高いところでは、船体の速度が速いが空気が薄いので熱流束は低い。また、高度が低いところでは、空気が濃いが船体の速度が遅いので熱流束は低い。ということで、どこか途中のところでピークがあるような分布となっている。私の計算では最大熱流束は260 (W/cm2)となっている。前回のblogで参照したApolloの簡単な資料によると、最大熱流束は240 (W/cm2)とあるので、簡易的な計算でも十分使えるのが分かる。
close all
clear variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial values
Ve_ms = 10e3; % Velocity at entrance to earth's atmosphere (km/s)
h0_m = 76.2e3; % Initial altitude (250 000 feet)
gamma0 = -5.9/180*pi; % entry angle
% Solve differential equation
t_range = [1, 3e2];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,y] = ode45(@R415, t_range, [Ve_ms; h0_m; gamma0;], opts);
wh = y(:,2)>=0;
t_s = t(wh)'; % time (s)
V_ms = y(wh,1)'; % velosity (m/s)
h_m = y(wh,2)'; % altitude (m)
gamma_deg = y(wh,3)'/pi*180; % angle (degree)
% Allen-Eggersのexponential atmosphere model
rho_kgm3 = exp_rho(h_m);
% Sutton-Gravesのconvective heating model
k = 1.7415e-4; % Earth
Rn_m = 3; % Nose radius (m)
qs = k*sqrt(rho_kgm3/Rn_m) .* (V_ms).^3 / 1e4; % W/cm2
%PLOT
subplot(3,1,1)
plot(rho_kgm3, h_m*1e-3)
xlabel('density (kg/m3)')
ylabel('Altitude (km)')
grid on
subplot(3,1,2)
plot(V_ms*1e-3, h_m*1e-3)
xlabel('Velocity (km/s)')
ylabel('Altitude (km)')
grid on
subplot(3,1,3)
plot(qs, h_m*1e-3)
xlabel('Convective heating (W/cm2)')
ylabel('Altitude (km)')
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = R415(t, y)
% y(1): velocity (m/s)
% y(2): altitude (m)
% y(3): angle (rad)
Cb_kgm2 = 395.8 ; % ballistic coefficient (kg/m2)
g= 9.8; % gravity acceleration(m/s2)
Re = 6.3781e6; % Earth radius (m)
dydt = [-0.5*exp_rho(y(2))*y(1)^2/Cb_kgm2 - g*sin(y(3));
y(1)*sin(y(3));
(y(1)/(Re+y(2)) - g/y(1))*cos(y(3))];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rho_kgm3 = exp_rho(h_m)
% h: altitude (m)
rho0_kgm3 = 1.752 ; % 0.0034 slugs/ft3
h0_m = 6705.6; % 22000 (1/ft)
rho_kgm3 = rho0_kgm3*exp(-h_m/h0_m);
end
大気圏突入時の熱流束を知りたいー3. 落下してくる宇宙船の運動モデルー
blunt body（Appolloの宇宙船カプセルみたいな形の物体）が落下してくる時の、速度を計算する方法を調べた。
Apollo Reentry GuidanceのAppendix Aにある微分方程式を少し修正して、以下の式を得る。
元の式の船体角度γの微分方程式には、lifting parameterを含んだ項があるのだが、上の式では無視している。根拠として薄弱なのは承知だが、International Journal of Aerospace Engineeringの"High-Fidelity Aerothermal Engineering Analysis for Planetary Probes Using DOTNET Framework and OLAP Cubes Database"という論文で、lifting parameterを省略しているのを参考にさせてもらった。現実問題としては、Apolloの簡単な資料にlifting parameterが載っていないので、同パラメータの項を無視できるのは助かった。
上述の式で、同Apolloの資料にあるパラメータ（Apollo 6、突入速度10 km/s、突入角度-5.9°、ballistic係数395.8 kg/m2）を用いて計算した結果が次の図である。大気密度は前回のblog記事の式を用いた。縦軸が高度、横軸が速度で、各高度での速度が分かるようになっている。ちなみに、この図は上述のInternational Journal of Aerospace Engineeringにある図と整合している。
close all
clear variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial values
Ve_ms = 10e3; % Velocity at entrance to earth's atmosphere (km/s)
h0_m = 76.2e3; % Initial altitude (250 000 feet)
gamma0 = -5.9/180*pi; % entry angle
% Solve differential equation
t_range = [1, 3e2];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,y] = ode45(@R415, t_range, [Ve_ms; h0_m; gamma0;], opts);
wh = y(:,2)>=0;
t_s = t(wh)'; % time (s)
V_ms = y(wh,1)'; % velosity (m/s)
h_m = y(wh,2)'; % altitude (m)
gamma_deg = y(wh,3)'/pi*180; % angle (degree)
%Plot results
figure();
subplot(3,1,1)
plot(t_s, h_m*1e-3)
xlabel('time (s)')
ylabel('Altitude (km)')
grid on
subplot(3,1,2)
plot(t_s, V_ms*1e-3)
xlabel('time (s)')
ylabel('Velocity (km/s)')
grid on
subplot(3,1,3)
plot(t_s, gamma_deg)
xlabel('time (s)')
ylabel('angle (deg)')
grid on
%Plot results
figure();
plot(V_ms*1e-3, h_m*1e-3)
xlabel('velocity (km/s)')
ylabel('Altitude (km)')
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = R415(t, y)
% y(1): velocity (m/s)
% y(2): altitude (m)
% y(3): angle (rad)
Cb_kgm2 = 395.8 ; % ballistic coefficient (kg/m2)
g= 9.8; % gravity acceleration(m/s2)
Re = 6.3781e6; % Earth radius (m)
dydt = [-0.5*exp_rho(y(2))*y(1)^2/Cb_kgm2 - g*sin(y(3));
y(1)*sin(y(3));
(y(1)/(Re+y(2)) - g/y(1))*cos(y(3))];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allen-Eggersのexponential atmosphere model
function rho_kgm3 = exp_rho(h_m)
% h: altitude (m)
rho0_kgm3 = 1.752 ; % 0.0034 slugs/ft3
h0_m = 6705.6; % 22000 (1/ft)
rho_kgm3 = rho0_kgm3*exp(-h_m/h0_m);
end
大気圏突入時の熱流束を知りたいー2. 大気のモデルー
宇宙船の運動を考える前に、大気をモデル化する。具体的には、任意の高度での空気の密度をモデル化する。
ここでは、AllenとEggersの論文に出てくる指数モデルを紹介する。
同論文では、この単純なモデルで高度20,000から180,000フィート（6から55 km）の密度をよく表せるとしている。
大気圧のデータとしてはU.S. standard atmosphere, 1976というのがある。手持ちの理科年表にもこのUS標準大気のデータが載っているので、Allen-Eggersの指数モデルと比較してみた。下図では縦軸を高度にとってプロットしており、それなりに一致しているのが分かる。
close all
clear variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allen-Eggersのexponential atmosphere model
rho0_kgm3 = 1.752 ; % 0.0034 slugs/ft3
beta_km = 1/6.7056; % 1/22000 (1/ft)
AE.Z_km = linspace(0,60);
AE.rho_kgm3 = rho0_kgm3*exp(-beta_km*AE.Z_km);
% 1976, U.S. 標準大気
tmp= xlsread('./Z_vs_rho.xlsx');
US.Z_km = tmp(:,1)';
US.rho_kgm3 = tmp(:,2)';
%PLOT
semilogx(AE.rho_kgm3, AE.Z_km, 'Displayname', 'Allen-Eggers model')
hold all
semilogx(US.rho_kgm3, US.Z_km, '+', 'Displayname', 'US standard atmosphere, 1976')
xlabel('density (kg/m3)')
ylabel('Altitude (km)')
legend('show')
ylim([0,60])
grid on
Z_vs_rho.xlsxに保存したUS標準大気のデータは以下のとおり。
| Z (km) | rho (kg/m3) |
|---|---|
| 0 | 1.2250E+00 |
| 1 | 1.1117E+00 |
| 2 | 1.0066E+00 |
| 3 | 9.0925E-01 |
| 4 | 8.1935E-01 |
| 5 | 7.3643E-01 |
| 6 | 6.6011E-01 |
| 7 | 5.9002E-01 |
| 8 | 5.2579E-01 |
| 9 | 4.6706E-01 |
| 10 | 4.1351E-01 |
| 11 | 3.6480E-01 |
| 11.1 | 3.5932E-01 |
| 12 | 3.1194E-01 |
| 13 | 2.6660E-01 |
| 14 | 2.2786E-01 |
| 15 | 1.9476E-01 |
| 16 | 1.6647E-01 |
| 17 | 1.4230E-01 |
| 18 | 1.2165E-01 |
| 19 | 1.0400E-01 |
| 20 | 8.8910E-02 |
| 21 | 7.5715E-02 |
| 22 | 6.4510E-02 |
| 23 | 5.5601E-02 |
| 24 | 4.6938E-02 |
| 25 | 4.0084E-02 |
| 26 | 3.4257E-02 |
| 27 | 2.9298E-02 |
| 28 | 2.5076E-02 |
| 29 | 2.1478E-02 |
| 30 | 1.8410E-02 |
| 33.2 | 1.3145E-02 |
| 35 | 8.4634E-03 |
| 40 | 3.9957E-03 |
| 45 | 1.9663E-03 |
| 47.4 | 1.4187E-03 |
| 50 | 1.0269E-03 |
| 51 | 9.0690E-04 |
| 55 | 5.6810E-04 |
| 60 | 3.0968E-04 |
| 65 | 1.6321E-04 |
| 70 | 8.2829E-05 |
| 72 | 6.2374E-05 |
| 75 | 3.9921E-05 |
| 80 | 1.8458E-05 |
| 86 | 6.9580E-06 |
| 90 | 3.4160E-06 |
| 91 | 2.8600E-06 |
| 100 | 5.6040E-07 |
| 110 | 9.7080E-08 |
| 120 | 2.2220E-08 |
| 130 | 8.1520E-09 |
| 140 | 3.8310E-09 |
| 160 | 1.2330E-09 |
| 180 | 5.1940E-10 |
| 200 | 2.5410E-10 |
| 250 | 6.0730E-11 |
| 300 | 1.9160E-11 |
| 350 | 7.0140E-12 |
| 400 | 2.8030E-12 |
| 450 | 1.1840E-12 |
| 500 | 5.2150E-13 |
| 550 | 2.3840E-13 |
| 600 | 1.1370E-13 |
| 650 | 5.7120E-14 |
| 700 | 3.0700E-14 |
| 750 | 1.7880E-14 |
| 800 | 1.1360E-14 |
| 850 | 7.8240E-15 |
| 900 | 5.7590E-15 |
| 1000 | 3.5610E-15 |
大気圏突入時の熱流束を知りたいー1. はじめにー
とあるキッカケで大気圏突入時の宇宙船への入熱を調べたので、何回かに分けてblogに記録を残す。
宇宙船への熱の出入りは、influxとしては「大気との熱伝達」と「大気からの輻射熱」があって、outfluxとしては「宇宙船からの輻射熱」と「宇宙船表面の蒸発」（最近のは蒸発しないかも・・・）がある。「大気からの輻射熱」というのは、宇宙船によって断熱圧縮された高温の大気からの輻射熱だ。これら4つの熱束のバランスで最終的な入熱がきまるが、大気との熱伝達については1950年代くらいから研究があっていろいろなモデルがあるようなので、少し調べてみた。
Thermal & Fluids Analysis Workshop 2012のAerothermodynamics Courseというトークのスライドが参考になった。このトークの趣旨は、「今はシミュレーションでいろいろ計算できるが、その結果を鵜呑みにしてはだめ。シミュレーションの妥当性を確認するために、（いろいろな仮定は入っているが）理論的なアプローチで得られるモデルと比較することが大事なので、そういったモデルを勉強しましょう」というものだ。
それで、一番簡単な「blunt bodyが落下してくる時の、大気との熱伝達を評価する方法」をまとめることにした。blunt bodyというのは日本語にすると鈍頭（どんとう）物体とか言うようだが、Appolloの宇宙船カプセルみたいなやつだ。本当はスペースシャトルが落下してくるときのEquilibrium Glideモデルを知りたかったが、そこまで調査しきれなかった。
大気との熱伝達を評価するモデルは色々あるようだが、簡便な方法で評価するため、以下の3つのモデルを把握すればよい。それぞれのモデルについて、blog記事を分けて書いていくことにする。
MATLABでplotの線を透明にしたい
MATLABで配列xとyをプロットする時はplot(x,y)ですよね。plotのオプションを指定する時は、 Name,Value の引数ペアをコンマ区切りで加えていきます。例えば、点線で線の太さを2ポイントにしたければ、plot(x, y, ’LineStyle', ':', 'LineWidth', 2)です。
plotの線を透明化するときは、この手法は使えません。plot関数はchart line オブジェクトとよばれるオブジェクトを出力することができるので、そのオブジェクトを使って線を透明化します。例えば、次のようにやると透過度が0.2の線（0で透明、1で非透明）となります。
p=plot(x, y); p.Color(4)=0.2;
ちなみに、p.Color(1:3)にはRGBの3成分の値が格納されています。
別のやり方
chart lineオブジェクトを使わなくても実はできる。赤い線で半透明な線を引きたければ次のようにする。
plot(x, y, 'Color', [1, 0, 0, 0.2]);
ただ、このやり方はこちらが色のRGB値を指定しないので、自分としては煩わしいやり方です。特に複数の線をプロットするときに、色を自分で指定するのは煩わしいので、色は自動でプロットさせといて、プロット後にchart lineオブジェクトを使って透明化しています。
例
2021年1月1日から2021年12月17日までの水戸市の気温データをプロットしてみます。データは気象庁のホームページからダウンロードします。「過去の地点気象データ・ダウンロード」から、データの種類を「時間別」、項目を「気温」としてデータをダウンロードすると、data.csvというファイルを入手できます。このデータを何も考えずにプロットすると下図のようになって見にくいです。
そこで、1時間ごとのデータは半透明にしつつ、1週間の移動平均をとったデータを重ね書きすると下のようになります。
このグラフを作るためのコードは以下のとおりです。
clear variables
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = './data.csv';
HEADER_LINES = 5;
DELIMITER = ',';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+++ファイルの行数を調べる+++
fid = fopen(filename);
alldata={};
while ~feof(fid)
thisline = fgetl(fid);
if ~ischar(thisline); break; end
alldata{end+1,1} = thisline;
end
nol=numel(alldata); % nol = number of line
fclose(fid);
%+++データの数を決定し、配列を用意する+++
data_num = nol - HEADER_LINES;
time_string = cell(1,data_num);
T_deg = zeros(1, data_num);
%+++データを配列に格納する
for i=1:data_num
newStr = split(alldata{HEADER_LINES+i}, DELIMITER);
time_string{i} = newStr{1};
T_deg(i) = str2double(newStr{2});
end
%+++データをtimeseries オブジェクトにする
ts = timeseries(T_deg, time_string, 'Name', 'temperature (deg)');
ts.TimeInfo.Format = 'mmmdd';
%+++移動平均データの算出して、timeseries オブジェクトにする
windowSize = 24*7; %元のデータが1時間毎のデータなので、これで1週間の平均
b = (1/windowSize)*ones(1,windowSize);
a = 1;
ts_ave = timeseries( filter(b,a,T_deg), time_string, 'Name', 'averaged temperature (deg)');
ts_ave.TimeInfo.Format = 'mmmdd';
%+++描画
p=plot(ts);
p.Color(4) = 0.2;
hold all
plot(ts_ave, 'color', p.Color(1:3));
ylabel('Temperature (deg)')
核融合ベンチャーその４：Commonwealth Fusion System
今、核融合ベンチャーの中で一番熱いのがCommonwealth Fusion Systems (CFS) じゃないでしょうか。このベンチャーは、MITのPlasma Science and Fusion Centerが立ち上げたのですが、その経緯をAlcator C-Mod（アルカトール・シー・モド）のシャットダウンまで遡って説明したいと思います。
Alcator C-modの話
Alcator C-ModはMITにあった中型のトカマク装置なのですが、磁場が他の装置と比べものにならないくらい強いということで、とても変わった実験成果を挙げていました。他の装置で検証できないような結果もたくさんありますが、素晴らしい成果をたくさん挙げています。それなのに、この装置はシャットダウンされることになります。2006年のITER協定署名以降、ITER関連で多くの予算が必要になった結果、米国内のITER以外の核融合予算が減っていきました。その影響はAlcator C-modにもおよび、2013年にシャットダウンすることが決まります。ただ、2013年以降にも実験を再開したり、再開したものの優秀なエンジニアがいなくって困っているらしいとか、やっぱりまたシャットダウンするらしいとか、常に噂が流れて不安定な状況が続きます。そのころのMITは、『学生の教育のためにAlcator C-Modは必要だ！』という運動を盛んに行っていました。うる覚えですが、学生主導で署名サイトのようなものが立ち上げられていた気がします。そういう状況の中で、MITが「学生が新しい超電導材を使って画期的な炉設計をしたぞ」とプレス発表をします。炉設計というのはまだ絵に描いた餅なので、当時は「なぜ炉設計ぐらいでプラス発表するの？？？」という感じでした。炉の名前はアイアンマンのARCでふざけているし、学会誌がFusion Engineering and Designという敷居の低い学会誌だし。。。おそらくプレス発表としては「学生が」というところを強調したくて、「こんな素晴らしい教育機関の実験装置をシャットダウンするなんてバカなことは止めてよね」というメッセージなのかなぁと思って記事を読んでいました。ちなみに、その学生の炉設計はこの論文です（arXivにも論文あり）。まぁ、そんなこんながあったものの、結局Alcator C-modはシャットダウンします。
ARCの特徴
そして2018年にCommonwealth Fusion Systemsが創業され、「学生が設計したというあのARCを作るぞ」と言って登場するわけです。上述の通り、ARCのプレス発表はAlcator C-Mod延命戦略の一貫だと私は思っていたので、「ベンチャーを立ち上げるなんて、ARCの設計は本気だったのかぁ」ととても驚きました。それにしても、『国からの予算がもらえないなら、投資マネーを集めて自分たちで研究を続けるぜ』ってのは素晴らしいですよね。もちろんTAEやGeneral Fusionが市場を開拓してくれていたというのもあるのでしょうが、それにしても素晴らしいです。このARCの特徴は以下の通りです。
- 強磁場
- I-mode
- 分解可能
- 内部鎖交コイル
- 溶融塩ブランケット
ARCはプラズマの中心で9.2 Tもあります。強磁場ですね〜、MIT魂をつよく感じます。ITERだと5.3 Tですし、現役の最大の装置JETでは3.45 T、今日本で作っているJT-60SAでは2.25 Tです。先ほど述べた、他に類を見ないほど強磁場だといったAlcator C-Modでも7 Tです。ARCでは9.2 Tですから、めちゃくちゃ強力です。この強力な磁場を可能にしているのが、 REBCO超伝導体です。WikipediaによるとREBCO線材は2009年から流通しだしたということなので、新しい線材です。ITERの基本設計をしていた頃に、この線材があったらITERの設計もだいぶ変わっていたかもしれませんね。
そして、2番目の特徴のI-modeですが、これはAlcator C-Modで見つかったプラズマの閉じ込め状態です。ここでもMIT魂を感じますね。L-modeとH-modeの話は、Tokamak Energyのblogで話を超簡単にしましたので、気になるひとはそちらを見てみてください。このI-modeというのは、L-modeのようにELMはない（壁への負荷が低い）のに、H-modeのようにエネルギー閉じ込め効率が高い（なんだったらH-modeよりも高効率）というプラズマの状態です。I-modeになるには、高磁場であることが条件のようで、他の装置で再現できているのかが分からないのですが、ARCはこのI-modeを採用しています（私が知らないだけで、再現できているかもしれません）。
３番目と４番目は同じようなことなのですが、このARCは磁場コイルから真空容器まで、分解可能で部品を交換して再組み立てできる設計になっています。下に論文のgraphical abstractの図を載せました。オレンジ色でDの形をしたコイルが半割りになっていて、真空容器を取り出せるようになっています。これは画期的なことで、ITERや従来の原型炉の設計ではそういったことができません。Dの形をしたコイルはトロイダル磁場コイルと呼ばれる超重要なコイルなのですが、このコイルは真空容器を壊さないかぎり取り出せません。なので、ITERではものすごいR&Dを重ねて、コイルの信頼性を確認して設計を進めてきました。ARCのようにあとで分解して修正できるとなると、開発はかなり楽になるとおもいます。いくらでもトライ・アンド・エラーができますからね。
そして、分解可能だからこそ、ARCではトロイダル磁場コイルと鎖交するように、色々なコイルを設置できています。ITERや原型炉ではそんなことできないので、トロイダルコイルの外側にPFコイルというのをおきます。こうすると、プラズマとPFコイルとの距離ができてしまうので、PFコイルを強力なものにしないとプラズマを制御できないのです。下の図をみてください、ARCではトロイダル磁場コイルと鎖交するように内部にPFコイルを配置できるので、プラズマの近くにコンパクトなPFコイルを設置しています。これが本当に作れれば、かなり経済的な炉になると思います。
そして、最後の溶融塩ブランケットについてです。ARCでは真空容器を二重壁にして、二重壁の間を液体金属である溶融塩で満たします。溶融塩の中にはリチウムが入っていて、核融合反応で生じた中性子をリチウムにあてることでトリチウムを生成して燃料増幅を実現します。もちろん、中性子を受け止めることによって核融合で発生した熱の回収もできます。また、溶融塩は冷媒の役目もあって、真空容器を冷却しています。ITERや原型炉の設計では、ブランケットという交換可能なモジュールを数百個用意して、これを真空容器の内側に設置します。真空容器はけっして壊れないように気をつけて一生使い続け（なにせ分解できませんから）、ブランケットだけを数年ごとに交換していきます。ARCではそんなことしないで、溶融塩で真空容器が腐食してしまっても、数年で交換してしまえば問題ないと考えているようです。ITERの場合、真空容器はかなり大きくて、地震やディスラプションに耐えられるようにゴツく作ってあります。それに、ITERでは冷却水が40気圧もあるので、それに耐えるためにも作りがゴツイです。ARCでは、そこまで真空容器の設計をちゃんとしていないとは思いますが、かなりペラペラな作りに見えます。だから、交換してもそんな大したコストでないと考えているのだと思います。
ARCの疑問点
核融合反応は、プラズマのサイズが大きい方が実現しやすいのですが、ARCはコンパクトなんです。これって、とても不思議なことです。プラズマの性能（どれだけ核融合炉に近づいているか）を考えるときに、ローソン図というのを使います。ARCの性能をローソン図で見てみましょう。下の図はwikipediaに載っていたローソン図で、縦軸が密度n、エネルギー閉じ込め時間τ、温度Tという三つのパラメータの三重積nτTで、横軸が温度Tです。JT-60UとJETという装置はQ=1というライン上にのっていますので、『加熱のために入力したエネルギーと同程度の核融合反応を引き起こすことに成功している』というのがわかります。ITERはQ=1よりは上側にきていて、あとちょっとでreactor conditionに届きそうです。さて、ARCのパラメータは論文によるとn=1.3x1020 m-3、τ=0.64 s、T=14 keVなので、三重積nτTが11.6x1020 (m-3s keV)です。ローソン図に重ねてみるとITERよりもJETに近い状態です。どうして、ARCが「小型なのにITER並みの性能がある！」と報道されているのか分からないです。
ARCがコンパクトなのが不思議だと言いました。それはなぜかと言いますと、三重積のnτTを大きくするために、通常の原型炉は装置を大きくすることによって、τを大きくしているからです。装置が小さい原型炉というのは、かなり夢のまた夢の話のはずなんです。ARCがコンパクトなのは、『MIT方式の超強磁場方式によって、τが大きくなったのかな』とも思ったのですが、論文を読むとARCのエネルギー閉じ込め時間τは0.64 sなんです。これは、個人的には期待はずれでした。ITERは3.7 sもあって、現状の装置にくらべてメチャクチャ長いのですが、ARCは現状の実験装置程度なんです。密度も普通ですし。。。核融合炉では、『核融合反応で発生したエネルギーでプラズマ自身を加熱する自己加熱』というのをさせたいのですが、ARCではほとんど自己加熱できないんじゃないですかねぇ。これでは原型炉とは呼べないと思います。ARCは原型炉ではなくて、「コンパクトながらもJETを少し発展させたプラズマを作れる装置」と思えば、不思議ではなくなりますね。
Commonwealth Fusion Systemsの疑問点
CFSは、この前「REBCO製のD形のコイルで超強力な磁場生成に成功したぜ！」ってプレス発表をしました。この「D形形状で強磁場を発生させる」って、とても大変なことなんです（下手に作ると、コイルがフープ力でぶっ壊れます）。だから本当に凄いなとは思うのですが、発表の写真みると分割可能なコイルになっていないんですよ。上で述べたようにARCの特徴って分解可能ってところがかなりのウェイトを占めているので、分解可能なコイルを実証して欲しいですよね。コイルの製作に成功したので、2021年からSPARCという装置を作って技術検証し、2025年からいよいよARCを作り始めるようなのですが、このSPARCも分解可能ではないんですよ。。。どこで『分解可能かどうかの検証』をするのでしょうか。SPARCって、ただ磁場が強いだけのトカマクという感じでARCとはだいぶギャップがある（2025年は厳しい）と思うのは私だけでしょうか。ARCの真空容器とか、溶融塩中に浮いているコイルとか、ちゃんと設計するとかなり難しいと思う。数年で設計できるのかなぁ。。。
核融合ベンチャーその３：Tokamak Energy
核融合ベンチャー企業３社目は、イギリスのTokamak Energy社です。これまでに紹介した、TAE社、General Fusion社の後にはいくつかのベンチャー企業が立ち上がるのですが、その中で成果を挙げている企業（というか、資金を集められている企業）となるとTokamak Energy社だと思います。
Tokamak Energy社は、球状トカマクという方式の核融合炉を目指しており、TAE社、General Fusion社と比べるとかなり王道寄りの方式を採用しています。王道はトカマク方式です。球状トカマクはそのトカマクの一種なのですが、原型炉・商業炉を目指すと工学的に厳しい問題がたくさんあるので、核融合研究の王道にはなっていないです。球状トカマクの一般的な問題は最後で述べるとし、まずはTokamak Energy社に関する現状と疑問点をまとめてみます。
Tokamak Energy社は、現在ST40という装置で実験をしています。この装置を一言で言ってしまうと『STARTという球状トカマク装置の強磁場版』です。STARTというのは1990年代のUKAEA（英国原子力研究所）の装置で、磁場圧力で規格化した熱圧力（いわゆるベータ値）がトカマク装置の3倍以上も高いという驚異的な結果を叩き出しました。それで、これは凄いということで、UKAEAはSTARTをアップグレードしてMASTという装置を作り、PPPL（米国プリストン大学プラズマ研究所）はNSTXという装置を作ります。その後、MASTはMAST Upgradeに、NSTXはNSTX-Uへとアップグレードしていきます。つまり、球状トカマクというのはずっと研究が続いていて今は第３世代なのですが、Toakamk Energy社はまだ球状トカマク第１世代のSTART規模という状態です。それなのに、こんなに注目を得ているというのは不思議でならないです。
彼らの成果で学術誌にちゃんと発表されている最新のものは、2019年のPlasma Physics and Controlled Fusionの論文（On the confinement modeling of a high field spherical tokamak ST40）だと思います。それを見ると、ST40はまだL-modeプラズマという『エネルギー閉じ込め性能があまりよくないプラズマ』しか生成できていないです。L-modeプラズマで核融合炉を作ると装置が大きくなってしまうので、トカマクの人たちはH-modeプラズマという『エネルギー閉じ込め性能が高いプラズマ』で核融合炉を作ろうとしています。しかし、H-modeプラズマは、ディスラプションというプラズマが突然崩壊する現象が起きます。ディスラプションが発生すると、大きな電磁力が炉心にかかり炉心が壊れてしまう恐れがあるので、核融合炉ではディスラプションを1発でも起こしてはならないのです。ただ、それがとても難問なのです（現在、世界中で研究者が研究中です）。また、H-modeプラズマでは、ELM（エルム）と呼ばれる現象が発生し、トカマク周辺部のプラズマが間欠的に外に吐き出されます。エルムがあると、トカマク内の不純物を吐き出せるというメリットがある一方で、炉壁に大きな熱負荷が生じるという大きなデメリットがあります。なので、エルムの強度や発生周期を制御しながら、H-modeプラズマを維持したいのですが、その制御方法というのは今世界中で研究中で、いくつかの案はあるもののまだ決定打が出ていないという状況です。トカマクはこのH-modeと向き合った研究をしているのですが、Tokamak Energy社はL-modeプラズマなので最先端の研究の入り口に立っていません。今後、Tokamak Energy社の研究が進んでH-modeプラズマを相手にするかもしれませんが、そうなると既存の先行研究との違いがよく分からなくなってきます。つまり、Tokamak Energy社はH-modeプラズマの核融合炉をどのように成立させるかの道筋を見せていません。それなのに、ベンチャー企業として評価されて資金を集められているというのが不思議です。
これは、詳しい事情を全く知らない人間が推測だけで思っていることですが、Tokamak Energy社ってプラズマ合体の実験をしたいだけの装置の気がするんですよねぇ。STARTとMASTというのはプラズマ合体という変わったプラズマのつけ方をします。このプラズマ生成法を推しているのが、STARTの立ち上げに貢献したM. Gryaznevich氏とA. Sykes氏で、彼らはTokamak Energy社のメンバーです（Gryaznevich氏は取締役員の一人です）。このプラズマ合体法は、球状トカマク第３世代のMAST Upgradeでは採用されていないんです。本当に球状トカマクの最先端の研究をしたければ、MAST Upgradeで研究をすれば良いのであって、START規模の装置に戻って研究することないと思うんですよねぇ。だから、「球状トカマクの研究をしたい」とか、「球状トカマクの核融合炉を実現したい」とかじゃなくて、「プラズマ合体法をやりたい」っというのが本音の動機じゃないのかなと思うのです（あくまで個人的な推測です）。
プラズマ合体法を使うと、トカマク生成直後にそれなりに高温のプラズマが出来ます。Tokamak Energy社のホームページをみると、『2021年のアメリカ物理学会でST40が球状トカマクの分野では最高となる温度を達成したことを発表した』とありますので、数keVぐらい出したのかもしれません。それはとても良いのですが、その高温のプラズマをどうやって維持するのかが難しいところなので、「最大瞬間風速ですごい値出しました」というのは、やはり原型炉のことまで考えていない1990年代の研究っぽいなぁという気がしてしまいます。
A.Sykes氏が2000年代に「L-modeで核融合炉を作れば良いじゃないか」と言っていたのを聞いたことがあります。L-modeの問題は炉が大きくなってしまうことなのですが、球状トカマクならベータ値が高いのでL-modeでも炉がそこまで大きくならないだろうし、何よりディスラプションが無いのが良いという主張でした。ここからはまたまた推測でしかないのですが、『Tokamak Energy社にはA. Sykes氏が絡んでいる』ということと、上述のように『同社は、H-modeプラズマをどのように制御して核融合炉を実現させるかを明示していない』ということを考慮すると、Tokamak Energy社はL-modeの球状トカマク核融合炉を目指しているのかもしれませんね。（あくまで個人的な推測です。根拠薄弱でございます。）
トカマクも球状トカマクも、プラズマはドーナッツの形状（トーラス形状）をしています。球状トカマクでは、ドーナッツの穴の部分がとっても狭いです。そのため、「中性子遮蔽・エネルギー回収・トリチウム増殖」という重要な機能をはたすブランケットという部品を設置するスペースがないのです。何か画期的なブランケットのアイディアが必要です。また、球状トカマクはトカマクにくらべてコンパクトな炉になる（材料費がかからない炉になる）という主張がありますが、コンパクトになると単位面積あたりの受熱量が高くなるので、熱構造的には設計がむずかしくなります。Tokamak Energy社が、L-modeプラズマを想定しているのか、H-modeプラズマを想定しているのかはわかりませんが、いずれにしても『球状トカマクに対してブランケットをどう配置し、どう保守するのか、また熱構造的な成立性をどう考えているのか』をちゃんと発表して欲しいと思います。核融合研究者の目から鱗が落ちるような炉設計があれば、それは球状トカマクだけの話ではなくて核融合業界全体としての福音になるはずです。
