竜田揚助の数理科学解説所

数学、物理、について書いてます。

流れがわかる流体力学part7(円柱に働くマグヌス力)

今回は前に導入した複素速度ポテンシャルを用いて、回転する円柱に現れるマグヌス力を求めてみます。マグヌス力とは一般に回転する円柱や球に働く揚力のことで、物体が回転するのにつられて周りの流体が回転し、渦をなすことが原因です。このことからわかる通り、流体に粘性があることが必要です。ですが今回は簡単のため次の近似を行います。

「物体と流体の間の摩擦は考慮するが、流体間の摩擦は考慮しない。」

こうすることでベルヌーイの定理を用いることができ、議論か難しくならずに済むのです。この仮定の下でマグヌス力を導出していきましょう。

 

【目次】

 

今回取り扱う系の確認

今回は、

非圧縮性/完全流体/二次元

で考えていきます。ただし上でも述べた通り、今回は簡単のため完全流体でありながらも物体と流体の間の摩擦はあって流体は物体につられて回転しようとするものとします。

具体的にはような状況を考えます。

f:id:Hannak:20211014212610p:plain

つまり、一様に流れる流体中に回転する円柱が入っていて、円柱自身は並進運動しない状況を考えていくのです。

 

マグヌス力を求める

一様流れ中の球の速度ポテンシャルを求めた時と同様の方法を用いて複素速度ポテンシャルを求めていきます。基本は同じことの繰り返しですのでここではすこし簡便に求めていきます。

まず

f:id:Hannak:20211014212617p:plain

という複素速度ポテンシャルに注目すると、

f:id:Hannak:20211014212625p:plain

が流速ベクトルの各成分になります。そこから、動径方向に垂直な成分の速度は、

f:id:Hannak:20211014212633p:plain

であるとわかりますから、これが静止流体中にx軸方向に−Uで動く円柱がある時の複素速度ポテンシャルとなるためには、その円柱の中心が原点に一致しているとして、

f:id:Hannak:20211014212639p:plain

となれば良いので、結局静止流体中にx軸方向に−Uで動く円柱がある時の複素速度ポテンシャルは、

f:id:Hannak:20211014212648p:plain

であれば良いです。さらに円柱の回転につられて動く流体の複素速度ポテンシャルを表すには、

f:id:Hannak:20211014212652p:plain

という渦を表す複素速度ポテンシャルを用いれば良いでしょう。

さらに、流体の一様流れなども考慮して考えると、結局この系の複素速度ポテンシャルは、

f:id:Hannak:20211014212656p:plain

で与えられます。ここで

f:id:Hannak:20211014212701p:plain

であることに注意して式を変形していくと、

f:id:Hannak:20211014212710p:plain

となりますから、速度ポテンシャルは、

f:id:Hannak:20211014212714p:plain

となり、これを偏微分して流速ベクトルの各成分を求めると、x成分は、

f:id:Hannak:20211014212720p:plain

となります。ここで、

f:id:Hannak:20211014212728p:plain

であることを用いました。この調子でy成分のほうも計算すると

f:id:Hannak:20211014212733p:plain

となります。ここで、

f:id:Hannak:20211014212742p:plain

を用いました。以上から、r=aとなる位置(つまり円柱の表面)での流速は、

f:id:Hannak:20211014212749p:plain

です。ここで今回考えている系ではベルヌーイの定理が成り立ちますので、これを用いると、

f:id:Hannak:20211014212814p:plain

となり、円柱表面での流体の圧力がもとまりました。状況を図に整理すると、

f:id:Hannak:20211014212818p:plain

となります。今注目したいのは、圧力由来の上向きの力なので、積分によって、

f:id:Hannak:20211014212828p:plain

となることが分かります。この上向きの力がマグヌス力と呼ばれるもので、密度、渦の強度、流体が流れてくる速度が大きいほど強いマグヌス力が生じることが読み取れます。

 

以上でマグヌス力についての解説は終わりです。次回は渦糸について見ていきます。

 

【次】

準備中

 

【前】

hannak.hatenablog.com

 

流れがわかる流体力学part6(複素速度ポテンシャルの導入と具体例)

今回は、前回までで見てきた速度ポテンシャルと流れ関数を組み合わせて、複素速度ポテンシャルというものを定めます。そして複素速度ポテンシャルに具体的な関数を与えることで、様々な流れを見ていくことにいたします。

【目次】

 

複素速度ポテンシャルの導入

複素速度ポテンシャルとは、速度ポテンシャルΦと流れ関数ψを用いて、

f:id:Hannak:20211014120133p:plain

と定義される関数のことです。この関数は、速度ポテンシャルと流れ関数の性質より、

f:id:Hannak:20211014120139p:plain

という式を満たしますが、これは複素解析で出てくるコーシー・リーマンの方程式にほかなりません。複素解析によると、複素関数wがコーシー・リーマンの方程式を満たすというのは、wが正則であるということと同値であり、

f:id:Hannak:20211014120146p:plain

をzで微分した結果は、xで微分した結果に一致するという事実が知られています。

ここで述べた複素解析の知識は、過去の記事に簡潔に述べられてますので必要に応じてご覧ください。(複素微分(コーシーリーマン方程式/複素偏微分))

今述べた正則関数の性質から、

f:id:Hannak:20211014120152p:plain

と計算されるので、実部と虚部を比べて、

f:id:Hannak:20211014120158p:plain

がわかります。この式は、複素速度ポテンシャルから流速ベクトルの情報を引き出すことができることを言っています。ですので複素速度ポテンシャルが与えられれば微分することで流速ベクトルがわかり、流体の流れを知ることができます。そこで次節では具体的な関数を与えて流れを求めていくことにします。

 

複素速度ポテンシャルの具体例

本コラム第3回で、基本的な流れを表す速度ポテンシャルを具体的に見ましたが、今回はそこで扱った具体例を、複素速度ポテンシャルの場合に拡張していきます。

 

一様流れ

流体が一様に流れていく様子を表す複素速度ポテンシャルは、

f:id:Hannak:20211014120207p:plain

です。実際に計算してみると、流速ベクトルは

f:id:Hannak:20211014120214p:plain

となっていますので、確かにこれは(a,b)方向へ一様に流れる流体を表しています。

 

 

 

十字のついたてがある場合の流れ

f:id:Hannak:20211014120247p:plain

上図のようにx軸y軸と位置する十字のついたてがある系の流れを表す複素速度ポテンシャルは、

f:id:Hannak:20211014120220p:plain

です。実際に微分してみると、

f:id:Hannak:20211014120225p:plain

となりますから、流速ベクトルのx成分,y成分は、

f:id:Hannak:20211014120231p:plain

となり、これを整理すると下のような微分方程式になります。

f:id:Hannak:20211014120236p:plain

上の第二式より、

f:id:Hannak:20211014120241p:plain

というのが、流線の方程式なので確かに上図の流れをきちんと表していそうです。

 

湧き出し

原点にある湧き出し口から放射状に流体が流れ出てくる場合の複素速度ポテンシャルは、

f:id:Hannak:20211014120251p:plain

です。今回も実際に微分してみると、

f:id:Hannak:20211014120258p:plain

となりますので、半径rの円から出てくる流体の総量は、

f:id:Hannak:20211014120304p:plain

になり、mは流体が出てくる量を示すパラメータになっていてこれが正なら湧き出しで負なら吸い込みです。この時の流れを図にすると下のような感じでしょうか。

f:id:Hannak:20211014120310p:plain

 

二重湧き出し

上で考えた湧き出し吸い込みを下図のように二つ用意してみます。

f:id:Hannak:20211014120316p:plain

この時の複素速度ポテンシャルは湧き出しの複素速度ポテンシャルを二つ足し合わせて、

f:id:Hannak:20211014120321p:plain

と計算されます。このとき、吸い込み口と湧き出し口が非常に近いとして近似を二か所にわたって行いました。詳細は上の赤い字を見てみてください。

そして、

f:id:Hannak:20211014120327p:plain

というふうな極限を取ると、

f:id:Hannak:20211014120333p:plain

であることが分かります。速度ポテンシャルでこの式に対応するものは流体中を動く球が作り出した流れを求めるときに活躍したように、この複素速度ポテンシャルは後で同様の活躍をしてくれます。

 

湧き出しの流速ベクトルを一斉に90°だけ回転させると渦を巻いて回転することが分かると思います。このことを式の上で表現するには、複素速度ポテンシャルにi(虚数単位)を掛け合わせて、

f:id:Hannak:20211014120339p:plain

とすれば良いです。実際に微分してみると、

f:id:Hannak:20211014120345p:plain

となりまして、これを座標に書き込むと、

f:id:Hannak:20211014120351p:plain

となり、図中のように動径方向成分とそれに直行する成分とに分けることで、流体は円周を反時計回りに回転していることが分かります。上図によると、その流体の速度は、m/rであり、これを一周にわたって周回積分すると、

f:id:Hannak:20211014120356p:plain

であります。ここで最後のΓは、一般のに

f:id:Hannak:20211014120401p:plain

で定義され、循環と呼ばれています。渦の強さを表す量だととらえるとよいでしょう。

二つ上の式からmを循環を用いて表すと、

f:id:Hannak:20211014120408p:plain

となり、これを用いて初めの式を書き換えると、

f:id:Hannak:20211014120412p:plain

となります。

これが渦を表す複素速度ポテンシャルで、原点以外の点では

f:id:Hannak:20211014120426p:plain

のように渦度は0ですから、たしかに速度ポテンシャルΦを定義することができます。

(このように一部分に渦があっても渦度0の場所を考える限りは速度ポテンシャルを定義できますから、複素速度ポテンシャルを使って議論を進めることが可能なのです)

 

以上でざっと具体例を見終えました。次回はこれを用いてマグナス効果を見ていきます。マグナス効果といえば野球の変化球でしょうか。

 

【次】

準備中

 

【前】

hannak.hatenablog.com

流れがわかる流体力学part5(流れ関数の導入とその意味)

今回からしばらくは二次元平面上の流体の動きを考えていきます。二次元の流体力学においては流れ関数という関数が非常に重要な役割を果たしますので、今回はその流れ関数というものを導入し、直感的な意味合いを紹介いたします。

 

【目次】

 

 

今回扱う系の確認

今回は以下の条件を満たす系を考えていきます。

非圧縮性/渦なし/二次元

渦なしとは渦がない状況で、流速ベクトルのローテーションが0であることが条件でした。

 

流れ関数の導入

まず今回は非圧縮性の流体なのでρ=const.となります。そのため、連続の式は、

f:id:Hannak:20211014115519p:plain

となります。(ρは密度、は流速ベクトル)この式を満たすためには、ある関数ψを用いて、

f:id:Hannak:20211014115526p:plain

と表せればよいです。これは実際に代入して確かめられます。さらに今回は渦度=0の場所を考えていますから、

f:id:Hannak:20211014115534p:plain

となり、速度ポテンシャルと同じようにラプラス方程式を満たします。ここで頭の整理のために速度ポテンシャルと流れ関数を比較してみます。

f:id:Hannak:20211014115541p:plain

この表からわかる通り、この二つはかなり似ていますね。

今のところはこの流れ関数のご利益を受けていませんが次回に複素速度ポテンシャルという、Φ+iψで表される関数を考えていくときっとありがたみがわかります。

 

流れ関数の意味

このままでは流れ関数についてあまりつかみどころのないまま終わってしまいますので、最後になかなか面白い性質を見てみます。

まず二次元平面上に任意の点A,Bをとり、二点を結ぶ任意の曲線を取ります。この状況を図にするとちょうど下のようになります。このとき、曲線ABを上から下へ横切っていく単位時間当たりの流体の量を求めていきます。

f:id:Hannak:20211014115549p:plain

分かりやすいように線素をとって拡大すると、

f:id:Hannak:20211014115559p:plain

となり、黄色い部分が微小時間Δtの間で通過する流体です。ですので、これらを足し合わせてΔtで割れば、求める値は

f:id:Hannak:20211014115604p:plain

となります。驚くべきことに通過する流体の体積は点BとAの流れ関数の差で表されています。

このことから、「流体は流れ関数が一定となる曲線を横切らない」ことが分かるので、

「流れ関数が一定の曲線=流線」

が成り立つとわかります。このことが分かればだいぶ馴染みやすいように思えますね。以上で流れ関数についての話は終わりましたので、次回は複素速度ポテンシャルを導入しそれに具体的な関数を与えてみることにします。

 

【続き】

準備中

【前の記事】

hannak.hatenablog.com

ルベーグ積分part2(ルベーグ積分の直感的導入)

 前回はルベーグ測度を直感的に導入しましたが、今回も引き続き直感的にルベーグ積分のほうを見ていきたいと思います。(より厳密な議論は次回以降きちんと行います。)

まずは、積分のイメージです。ルベーグ積分について簡単に述べているウェブサイトやyoutubeを見たことがある方はご存じかと思いますが、ルベーグ積分の場合は積分区間を縦ではなく横に割っていきます。図にするとちょうど下のような感じですね。(今回は簡単のため一変数関数を扱っていきます。)

f:id:Hannak:20210907170817p:plain

縦軸を適当な区間に分けると、その各々の区間に対応する横軸の領域はいくつかの区間にまたがって存在していますが(上図で黄色く塗ったA₃などがそうです。)これらをまとめてひとつの集合だとします。

そして、横に分割している赤い線を目安にして下図のようなデジタルっぽいかくかくのグラフ(青線)を取ります。

f:id:Hannak:20210907170821p:plain

このかくかくのグラフをφ(x)と名付けると、

f:id:Hannak:20210907170824p:plain

という関数を用いて、

f:id:Hannak:20210907170827p:plain

で書き表すことができます。(a₁a₂a₃…は上図のように定めました。)

そして肝心な積分はというと

f:id:Hannak:20210907170830p:plain

で定義されます。これがルベーグ積分です。

 この式で気になるのは、m(dx)が出てきている点です。このmはルベーグ測度のことです。リーマン積分できないような関数では上図におけるA₁などの領域の長さを普通の意味で測ることはできませんが、ルベーグ測度を用いることでこの困難を解決することができます。このことから、ルベーグ積分はリーマン積分の兄弟分的なノリで作られたものではなく、その根底にはもっと深い理論が存在しているのだということを感じることができます。

 これから続く記事ではこの理論をできるだけ簡潔かつ分かりやすく書いていくつもりですので、ぜひともご覧ください。

 

続き

準備中

 

前の記事

hannak.hatenablog.com

【読めばわかる一般相対論part13】水星の近日点移動

今回は一般相対論の検証となる現象として有名な水星の近日点移動を取り扱っていきます。水星の近日点移動というのはその名の通り水星の近日点がどんどんとずれて移動していくことですが、その移動の要因としては大半が地球の自転軸のぶれによるもので、後の要因は水星の周囲の惑星からの引力であるということがわかっています。

 しかし、この二つの要因から導き出される計算結果は実際の近日点移動よりも値がすこしばかり小さく、第3の要因があるのではないかと考えざるを得ません。

 実は、その第3の要因というのが一般相対論的な効果なのです。本稿ではこのことについてみていきます。

 

【目次】

 

水星の近日点移動とは

まずは水星の近日点移動とはなにかを簡潔に述べますと、それはつまり水星の近日点が下図のように、少しずつ回転していく現象のことです。

f:id:Hannak:20211011175853p:plain

高校では惑星の軌道は左のような楕円であると習いましたが、実際には右のように動いているのです。(上の図はかなり誇張しています。)

 

相対論的効果による水星の近日点移動

太陽を球対称静的質量分布とみなせばシュワルツシルト計量を考えることで現象を解くことができます。

シュワルツシルト計量は

f:id:Hannak:20211011175623p:plain

と表せました。(太陽の質量をMとしています。)

一般に計量テンソルと線素は

f:id:Hannak:20211011175633p:plain

という関係にあります。ですから、シュワルツシルト解より、

f:id:Hannak:20211011175636p:plain

が計量テンソルであるとわかります。

今回は水星は重力の作用しか受けないので一般相対論においては外力0となります。(一般相対論において重力の効果は時空間のゆがみに含まれるからです。)

ですから、ニュートン力学運動方程式の相対論バージョンである、測地線の方程式は、

f:id:Hannak:20211011175641p:plain

となり、今後はこれを考えていけばよさそうです。計算を楽にするために、この現象がθ=π/2なる面上で起こるとすると、θ=π/2をシュワルツシルト解に代入して、

f:id:Hannak:20211011175700p:plain

となります。

そして、測地線の方程式を求めていくためにクリストフェル記号(上式のΓ)を求めていきます。クリストフェル記号は、

f:id:Hannak:20211011175645p:plain

で定まっていますので、丁寧に計算していくと、まずはi=0のときについて、

f:id:Hannak:20211011175703p:plain

なので、(j,k)=(0,1), (1,0)以外は0となり、

f:id:Hannak:20211011175649p:plain

となります。よって、

f:id:Hannak:20211011175653p:plain

となり、すなわち、

f:id:Hannak:20211011175657p:plain

という保存量をえました。この調子でi=3のときも考えていくと、

f:id:Hannak:20211011175710p:plain

なので、(j,k)=(1,3), (3,1)の時以外は0であり、

f:id:Hannak:20211011175714p:plain

になりますので、測地線の方程式は、

f:id:Hannak:20211011175735p:plain

と変形できます。よって、

f:id:Hannak:20211011175754p:plain

になります。以上で求めた二つの保存量αとβを用いて、シュワルツシルト解を変形していくと、

f:id:Hannak:20211011175759p:plain

となります。ここでの計算は、(万有引力運動方程式から楕円軌道を導くときのように)rとφの関係式に持っていこうという気持ちでやっていきます。上の式を少し整理していくと、

f:id:Hannak:20211011175804p:plain

になります。そしてこれをφで微分することでシンプルな式を得ます。:

f:id:Hannak:20211011175808p:plain

これが水星の軌道が満たすべき式になりますがまともに解こうとするとかなり大変ですので(そもそも解くことができないように思えます。)近似を2段階に分けて行っていきます。

【最初の近似】rは太陽から水星までの距離とは違いますが、太陽から水星までの距離と同じように、rなどと比べるとかなり大きいことは確かです。そこで①の右辺第三項を無視するという近似を行うと、

f:id:Hannak:20211011175812p:plain

となりこれの解は、単振動と同様にして

f:id:Hannak:20211011175817p:plain

になることが分かります。

【第二の近似】このままだとただの楕円軌道のままですから、そこで実際には余弦の位相がφではなく、δφだったとします。(このδによって近日点移動の分を表そうという算段です)

f:id:Hannak:20211011175822p:plain

式にするとちょうど上のようになります。これを①に再び代入して、δを求めていきます。

f:id:Hannak:20211011175826p:plain

上式において、余弦の係数を見比べると、

f:id:Hannak:20211011175837p:plain

と計算されますから、余弦の位相がちょうど2πだけ変化するのは、

f:id:Hannak:20211011175841p:plain

となるときです。これはつまり、

f:id:Hannak:20211011175844p:plain

であることをしめしており、つまり水星の近日点移動は一周するたびにGM/lc²だけずれます。

 実はこの式に実際の値を代入すると、本稿の初めに述べた水星の近日点移動のうち原因がわかっていなかった分の数値にぴたりと一致します。このことから、この現象は一般相対論の正当性を示すものの一つとして非常に有名であるのです。

 

以上が水星の近日点移動の一般相対論的効果の内容です。これで「読めばわかる一般相対論」と題したシリーズでの相対論の解説は終わりになりますが、ほかにも重力波宇宙論などの応用的な内容についても今後別の機会に見ていこうと思います。

 

前の記事

hannak.hatenablog.com

 

ルベーグ積分part1(ルベーグ測度/可測集合の直感的な導入)

今回から何回かにかけてルベーグ積分を扱っていきます。なかなか難しい内容ですのでまずは直感的な説明を行い、次回以降その詳細や証明を見ていこうと思います。

 

【目次】

 

ジョルダン測度の考え方

ルベーグ積分の前にまず、ジョルダン測度の考え方をざっと見ていきます。(リーマン積分という高校で習った積分ジョルダン測度の考えに基づいています。)

下図のような赤い部分(集合Sとします。)の面積を求めたいとき、ジョルダン測度の考え方にのっとると、この図形を’有限個の長方形の面積の和’で近似していくことになります。

f:id:Hannak:20210907170712p:plain

f:id:Hannak:20210907170717p:plain

のように長方形をI(アイ)で表すとして、その面積(| I |とする。)は小学校で習った公式を思い出すと、

f:id:Hannak:20210907170722p:plain

であります。このような長方形のうちSに含まれるものだけを集めて、面積の和を取れば、Sの面積の大まかな近似になりまして、

f:id:Hannak:20210907170726p:plain

になります。このときの長方形の取り方をいろいろに工夫してみると、上式の値は様々に変わるでしょう。このようにしてできる様々な値の上限を

f:id:Hannak:20210907170731p:plain

のように※付きで表し、ジョルダン内測度と呼びます。定義の仕方よりジョルダン内測度は’’実際の面積’’より大きくなることはないでしょう。(なぜなら、Sに完全に含まれている長方形の面積の和を考えているからです。)

次に’’実際の面積’’より少し大きめに見積もったものも考えます。そのためには、大盤振る舞いにSと少しでも共通部分を持った長方形をすべて足し合わせてしまえばよいでしょう。式にすると下です。

f:id:Hannak:20210907170735p:plain

これはいくらでも大きくはなりますが、下には有界で、その下限を

f:id:Hannak:20210907170739p:plain

のような下付きの※を用いて表し、ジョルダン外測度と呼びます。

もしも、小さめに見積もったジョルダン内測度と大きめに見積もったジョルダン外測度が一致すれば、この見積もりは’’正しい’’ことになり、

f:id:Hannak:20210907170742p:plain

として、ここで求めた内測度か外測度のどちらかを(両方等しいのでどっちでもいいのですが)Sの面積とすればよいでしょう。

以上がジョルダン測度のおおよその考え方です。しかしこの考え方には欠点があり、例えば図形の内部に可算無限(無限だが数えられる量です。例えば自然数。)個の穴がまんべんなく偏在していればジョルダン内測度は0となってしまい、外測度は0以外の有限の値になってしまうので面積は定義できなくなってしまいます。その欠点を解消する考え方がルベーグ測度です。以降はルベーグ測度についてみていきましょう。

 

ルベーグ測度の考え方

まず、前節の最初の段階で図形を有限個の長方形で近似したかと思いますが、今回は近似の精度を上げるため使う長方形を可算無限個に増やします。そして、下図のように図形を覆ってしまいます。(絵では可算無限個は表現できないため有限個になっておりますがご勘弁を。)

f:id:Hannak:20210907170746p:plain

そして、作った長方形のうち、Sと少しでも共通部分を持つものの面積ををすべて足し合わせていきます。(つまりちょっと大きめに値を計算するのです。)そして、ジョルダン測度の場合と同様にこれの下限を取ってみます。

f:id:Hannak:20210907170752p:plain

これをルベーグ外測度といいます。

この調子でルベーグ内測度も定めたいところですが、この方法を愚直に内測度にも適用していくとジョルダン測度の欠点を引き継いでしまいます。そこで、「図形の内部から面積を求めていく」のはあきらめて、「Sの外枠の面積(下図で赤色の部分)をやや大きめに求めて、長方形の面積から引く」ことでSの小さめな面積を求めていきます。

f:id:Hannak:20210907170753p:plain

この赤い外枠の面積を大きめに求めるには、この部分の外測度を求めればよいでしょう。

f:id:Hannak:20210907170757p:plain

そして、長方形の面積から外測度m*を引いた値を内測度として、上式のように定めます。以上で内測度と外測度がわかりました。そうしましたら、ジョルダン測度の時と同様に、内測度と外測度が一致するときが面積を上手く求めたということになります。そして、この値をルベーグ測度といい、Sの面積になります。

f:id:Hannak:20210907170804p:plain

のように、内測度=外測度という条件は、上の黒く囲った式がSを完全に含む任意の長方形Iについて成り立つときであると言えます。

そして、ルベーグ外測度(今後外測度といったときはルベーグ外測度のこととします)は

f:id:Hannak:20210907170810p:plain

という性質を満たしてくれます。証明したいところですが今回は直感性に重きを置きたいので割愛します。

以上によってルベーグ測度の導入が済んだのですが、実はカラテオドリさんによって、

f:id:Hannak:20210907170814p:plain

が上で導出した「外測度=内測度」であるための必要十分条件であることが証明されております。上式において「長方形の面積は’’縦×横’’」という小学校の先生が言っていた謎の決まりは使わずに済みますので、むしろこちらを外測度=内測度となる条件式として定義するべきです。そしてこの条件を満たしてくれる集合Sを可測集合と定義します。

 

今回は以上で終わりとします。次回はざっくりとルベーグ積分を導出し、その次の回からはいよいよある程度の厳密性を保ちつつもわかりやすいようにルベーグ積分の解説を始めていきたいと思います。どうぞ引き続きご覧ください。

 

次の記事

hannak.hatenablog.com

 

サイトマップ

サイトマップ - 竜田揚助の数理科学解説所