009_HeatInduction.mw

誘導加熱された円筒

イントロダクション

> restart;

このワークシートでは、以下の絵のような長い金属の円筒について考えます。

円筒は軸に対して平行な磁場を持っており、特定の深さで電流密度の制約があります。

今回は、時間と非常に長い円筒の表面にからの深さに温度がどのように依存するかを電流密度、パワー密度、および偏微分熱伝導方程式を使用し説明を行います。

[Inserted Image]

加えて、加熱は渦流が交互に磁場を上昇されることで引き起こされます。

表皮効果や磁場の周波数への依存は高電流密度に作用して材料の表面を加熱します。

深さの増加に従ってそれは急減します。表面の有効な磁場強度が与えられた時間で与えられた温度に達するにどれくらい必要であるかを調べます。

問題の定義

この例では、直径 6cm のアルミニウムの円筒の熱誘導について下記の問題を考えます。

表面の有効磁場強度の値 (HS, [A/m]) がいくつのとき表面の温度が 60 秒で 800K に達するか?

さらに、表面から 0.5cm の距離の電流密度 (J、[A/m^2]) は表面と比べ少なくとも半分の値を持たなければなりません。

アルミニウムの材料特性

  電気伝導率:sigma - 3.82*10^7 [1/(Ohm*m)]
  相対浸透率:mur - 1

  熱伝導率:lambda - 237 [W/(K*m)]

  熱の容量:c - 888 [J/(kg*K)]

  密度:rho - 2700 [kg/m^3]

  シリンダの半径:R - 0.03 [m]

このワークシートでは、記号パラメータより数値を使用することを示すために変数名の終わりに A をつけて使用します。

> sigA:=3.82*10^7: murA:=1: lambA:=237: cA:=888: rhoA:=2700: mu0:=4*Pi*10^(-7): RA:=0.03:

解析

円筒が長いの場合、円筒の長さが無限であると同等に終端効果のない近似を行います。また、コイルの周囲に発生する磁場は、円筒の長さに全体で均一であるとします。この場合、円筒の実行電流密度(J)は次の公式で与えられます。

> J := (r, delta) -> HS*(-KelvinBer(1,(sqrt(2)*r)/delta)-KelvinBei(1,(sqrt(2)*r)/delta)-I*KelvinBei(1,(sqrt(2)*r)/delta)+I*KelvinBer(1,(sqrt(2)*r)/delta))/(delta*(KelvinBer(0,(sqrt(2)*R)/delta)+I*KelvinBei(0,(sqrt(2)*R)/delta)));

J := proc (r, delta) options operator, arrow; HS*(-KelvinBer(1, sqrt(2)*r/delta)-KelvinBei(1, sqrt(2)*r/delta)-I*KelvinBei(1, sqrt(2)*r/delta)+I*KelvinBer(1, sqrt(2)*r/delta))/(delta*(KelvinBer(0, sqr...

delta は通常、"進入の深さ" として以下の様に定義されます。

> delta1 := (omega, sigma, mur) -> sqrt(2/(sigma*omega*mur*mu0));

delta1 := proc (omega, sigma, mur) options operator, arrow; sqrt(2/(sigma*omega*mur*mu0)) end proc

ここで磁場内の円筒の周波数が omega [1/sec] の時、

> PD := (r, delta) -> (1*evalc(abs(J(r,delta,R))^2))/sigma;

PD := proc (r, delta) options operator, arrow; evalc(abs(J(r, delta, R))^2)/sigma end proc

> simplify(PD(r, delta));

2*(KelvinBei(1, 2^(1/2)*r/delta)^2+KelvinBer(1, 2^(1/2)*r/delta)^2)*HS^2/(sigma*(KelvinBer(0, 2^(1/2)*R/delta)^2+KelvinBei(0, 2^(1/2)*R/delta)^2)*delta^2)

一時的な温度分布において、極座標における熱伝導のための偏微分方程式の解決が必要です。

diff(T(r, t), t) = lambda*((diff(T(r, t), r, r))+(diff(T(r, t), r))/r)/(rho*c)

T(r,t): 温度 [K]

境界条件が近似的な断熱条件に近づく(すなわち、熱は全く円筒の表面から周囲に伝わらない)と仮定すると、全ての熱は円筒の表面を通って入るので、解析解を見つけることができます。

> T := (t,r) -> T0 + (2*PA*t)/(R*c*rho) + (PA*R*((r/R)^2-1/2-4*(sum((exp(-(BesselJZeros(1,i)^2*lambda*t)/(c*rho*R^2))*BesselJ(0,BesselJZeros(1,i)*(r/R)))/(BesselJZeros(1,i)^2*BesselJ(0,BesselJZeros(1,i))),i=1..n))))/(2*lambda);

T := proc (t, r) options operator, arrow; T0+2*PA*t/(R*c*rho)+1/2*PA*R*(r^2/R^2-1/2-4*(sum(exp(-BesselJZeros(1, i)^2*lambda*t/(c*rho*R^2))*BesselJ(0, BesselJZeros(1, i)*r/R)/(BesselJZeros(1, i)^2*Bess...

T0: t=0 秒における初期温度

PA: 表面上の熱入力 [W/m^2]

R: 円筒の半径

総和は i=1..infinity から求めるにもかかわらず、増加する時間に応じて指数関数が非常に小さくなります。従って、大抵の場合、総和は初めの約50 (n=50) の要素に制限しても十分である。表面発熱体 (PA) は浅い表層の中に表皮効果で熱を発生させるので妥当であると想定します。PA を求めるには、パワー密度は r=0..R の間隔で統合されるに違いない。

解決法

まず初めに、適切な周波数 (omA) を見つけなければいける為に、0.5cm の深さの電流密度が表面の値の半分とすると:

> omA := fsolve(
        evalf(
         simplify(
          subs(
           R=RA, (1*abs((1*J(RA, delta1(omega,sigA,murA)))/HS))/2
          )
         )
         =evalf(
          subs(
           R=RA, abs((1*J(RA-0.005,delta1(omega,sigA,murA)))/HS)
          )
         )
        ), omega=1000);

omA := 1004.131532

次に、表面のパワーロス (PAA) を計算します。

> PAAI := int(simplify(subs({R=RA, sigma=sigA}, PD(r, delta1(omA, sigA, murA)))), r=0..RA);

PAAI := int(0.2290944224e-5*abs(HS^2*(-1.*KelvinBer(1., 219.5492289*r)-1.*KelvinBei(1., 219.5492289*r)-1.*I*KelvinBei(1., 219.5492289*r)+1.*I*KelvinBer(1., 219.5492289*r))^2), r = 0 .. 0.3e-1)PAAI := int(0.2290944224e-5*abs(HS^2*(-1.*KelvinBer(1., 219.5492289*r)-1.*KelvinBei(1., 219.5492289*r)-1.*I*KelvinBei(1., 219.5492289*r)+1.*I*KelvinBer(1., 219.5492289*r))^2), r = 0 .. 0.3e-1)PAAI := int(0.2290944224e-5*abs(HS^2*(-1.*KelvinBer(1., 219.5492289*r)-1.*KelvinBei(1., 219.5492289*r)-1.*I*KelvinBei(1., 219.5492289*r)+1.*I*KelvinBer(1., 219.5492289*r))^2), r = 0 .. 0.3e-1)

印加磁場 (HSA) が円筒を 60 秒で 800K まで温度を加熱するのにどれ位の大きさの値が必要かを fsolve コマンドで計算されます:

> HSA := fsolve(
        subs({c=cA, rho=rhoA, lambda=lambA, R=RA, T0=293.15, PA=PAAI, n=50}, T(60,RA))=800,
       HS=200000);

HSA := 268874.6938

磁場には 268875 A/m の値があることがわかります。これを磁気誘導 B の約 0.48T と等しいとします。
表面のパワーロスは次の様になります。

> PAA := evalf(int(simplify(subs({R=RA, HS=HSA, sigma=sigA},PD(r, delta1(omA, sigA, murA)))),r=0..RA));

PAA := 298150.4660

解の可視化

温度増加を時間の関数としてプロットします。

> Tn:=simplify(subs({c=cA, rho=rhoA, lambda=lambA, R=RA, T0=293.15, PA=PAA, n=50}, T(t,r))):
plot3d(Tn,r=0..RA, t=0..60, labels=["r in m","Time sec","Temp. K"], axes=boxed);

[Plot]

> plot3d(Tn,r=0..RA,t=0..4,labels=["r in m","Time sec","Temp. K"], axes=boxed);

[Plot]

このプロットは、熱が初め 'ゆっくり' と円筒の中央へ向かっていき消散するかを明確に表します。

初めに、温度は表面だけ上昇し、次に数秒後から中央で上昇します。全ての時間において中央と表面の少しの温度差が確認できます。この温度差は一定の仮定が終了した後の集中一点に集まります。

要約

誘導効果による温度上昇の計算は多くの産業分野で重要です。この例では、Maple ではどのように計算ができるかに関して短く紹介をしています。時間に関する温度と非常に長い円筒の中の位置は、熱伝導の電流密度、パワー密度、および偏微分方程式を使用することでどう計算されるかを示します。

参考文献

Carslaw, H.S.; Jaeger, J. C.: Conduction of heat in solids; 2nd. ed., Oxford University Press

Joos, Georg: Lehrbuch der theoretischen Physik; 8. Auflage, Leipzig: Akademische Verlagsgesellschaft Geest & Portig KG