【2024年_Y Yokoyama】High-speed photoelastic tomography for axisymmetric stress fields in a soft material: temporal evolution of all stress components

High-speed photoelastic tomography for axisymmetric stress fields in a soft material: Temporal evolution of all stress components - ScienceDirect  

日本語訳

本研究では、フォトエラスティック・トモグラフィー（PT）と高速偏光カメラを用いて、軟材料の動的な軸対称応力場におけるすべての応力成分を再構築する新しい手法を提示する。 研究対象として、一時的な応力場再構築の例として静的および動的なヘルツ接触を扱った。

• 静的ヘルツ接触（固体球をゲルブロックに押し当てる場合）では、弾性率47.4 kPaのウレタンゲル内の全応力成分を、測定された光弾性パラメータを用いてPTにより再構築した。その結果を理論解と比較したところ、良好な一致を示した。
• 動的ヘルツ接触（球がゲルに衝突する場合）では、高速偏光カメラを用いてゲル内部の一過的応力場を再構築した。PTにより、せん断応力波および軸方向応力波を定量的に測定し、基盤上で異なる伝播速度を示すことが確認された。

本手法により、O(10⁻¹)～O(10¹) kPaの範囲にわたる応力場を、大変形を伴う状況でも同時に測定できることが可能となり、動的シナリオにおける急速に変化する応力テンソル成分を高精度に捉えられることが実証された。 さらに、計算された衝撃力のスケーリング則が理論予測と一致し、PTによる軟材料内の動的な軸対称応力場測定の正確性が検証された。  

（Introduction）

軟材料における応力を評価することは、バイオメディカル工学や細胞プリンティング応用の分野で非常に重要である。たとえばバイオメディカル工学では、脳動脈瘤の周囲における応力分布を明らかにすることが、その破裂に至るメカニズムを解明するために不可欠である［1–3］。また、無針注射によって生じる応力場を調べて疼痛レベルを評価するには、組織（やはり軟材料である）の応力測定が必要である［4,5］。さらに、液滴や液体ジェットの衝突によって生じる材料内部の応力分布を理解することは、多様な工学プロセスに関連しており、研究対象として大きな関心を集めている［6–10］。 **光弾性法（Photoelasticity）**は、材料内の応力分布を推定するためのよく知られた手法であり［11–14］、上述のようなシナリオで有用であると期待されている。光弾性法は、材料に生じる複屈折の原理に基づいている。複屈折とは、応力を受けた材料の局所的な屈折率が方向によって異なるという光学的異方性である。偏光光が応力を受けた物体を通過する際の強度を観察することで、互いに直交する二つの屈折率の差とその軸の方向を測定し、材料内部の応力状態を推定できる。 **応力光学則（Stress-optic law）**は、材料中の二次主応力差 σd の積分と透過光の位相遅れ Δ の間に比例関係があることを示す［12,15］。 $$\Delta = C \int \sigma_d(y) \, dy$$   ここで、係数 C（応力光学係数）は材料固有の値であり［12］、y はカメラの光軸方向を表す。二次主応力とは、y軸に垂直な断面内に現れる見かけ上の主応力である［16］。さらに、透過偏光の方向 φ は、二次主応力の方向に関連している。本研究では、この位相遅れと偏光方向を光弾性パラメータと呼ぶ。 光弾性法は、圧力センサなどの他の測定手法に比べて、非侵襲的に材料全体の応力場を測定できる利点を有する。デジタル画像相関法（DIC）［7,17］も、局所的な材料変位を解析し、その構成方程式を適用することで全体応力場を取得可能である。ただし、その手法では粒子をマーカーとして材料に混合する必要があり、また極薄のレーザーシートを生成する光学系が必要となる。それに対して光弾性法では、応力光学則が有効である限り、マーカーを導入して材料を汚染することなく、応力値の測定が可能である。また、実験装置も円偏光光源と偏光カメラがあればよい。 材料内の三次元応力場を測定するために設計された手法は**積分光弾性法（Integrated Photoelasticity）**と呼ばれる［14,18–20］。ここで「積分」とは、記録された位相遅れが、光線の通過経路に沿った各点での応力状態に依存する光弾性パラメータの積分値を表すことを意味する。従来、この積分光弾性法はガラスのような硬い材料に主に適用され、残留応力の可視化に用いられてきた［11–14,21,22］。一方、我々の先行研究では、軟材料における軸対称応力場に対して積分光弾性法の適用可能性を検証し［16,23］、ガラスの約100万分の1程度の剛性しか持たないゲルに対しても有効であることを示した。 上述のように、応力を受けた材料を透過する偏光光の光弾性パラメータ（位相遅れと偏光方向）は、その内部応力場の積分値と関連している（Fig.1参照）。したがって、材料内の応力場を得るには、測定された光弾性パラメータから内部応力場を再構成する必要がある。これは**光弾性トモグラフィー（Photoelastic Tomography）**と呼ばれ、Abenらによって集中的に研究された［24–26］。トモグラフィーは三次元物体の内部構造を非破壊で解析する手法として広く知られており、たとえばX線を用いたCTスキャンでは、透過X線の減衰から内部構造を決定する。この場合、各点の内部構造を特徴づける量はスカラー値であるため、比較的容易に実現できる。 一方、応力は各点において最大で6成分を持つテンソルである。得られる光弾性パラメータは（位相遅れと偏光方向の）2値のみであり、再構成すべきテンソル成分数に比べると著しく少ない。このため、一般に応力場のトモグラフィーは難しいが、応力場が軸対称である場合にはテンソル成分数が減るため、再構成が可能となる［27,28］。すなわち、円筒座標において $$\sigma = \begin{bmatrix} \sigma_{rr} & 0 & \sigma_{rz} \\ 0 & \sigma_{\theta \theta} & 0 \\ \sigma_{rz} & 0 & \sigma_{zz} \end{bmatrix}$$   ここで、σrr は半径方向応力、σθθ は周方向応力、σzz は軸方向応力、σrz はせん断応力である。 しかし、光弾性トモグラフィーはまだ軟材料の内部応力の全視野測定には適用されていない。加えて、光弾性法は主に定常的な応力場の測定に用いられてきた。これは、位相シフト法を用いて光弾性パラメータを測定する際に、同一応力場に対して複数枚の画像を取得する必要があるためである。しかし、カメラセンサ内に異なる角度の偏光子を組み込んだ偏光カメラの開発により、光弾性パラメータをワンショットで測定することが可能になった［29］。さらに、高速偏光カメラ［30］を用いることで、非常に高速な現象においても一過的な光弾性パラメータを測定できるようになった。ただし、これまで高速偏光カメラは、軟材料における動的三次元応力場を「可視化」する（すなわち光弾性パラメータを測定する）ために利用されてきたものの［31］、光弾性トモグラフィーに基づいて測定された光弾性パラメータから軟材料の応力を定量的に再構築した例は知られていない。 本研究では、上記の課題に対応するため、高速偏光カメラを用いて軟材料における動的な軸対称応力場のすべての応力成分を測定する光弾性トモグラフィーを適用する。光弾性トモグラフィーと高速偏光カメラを組み合わせることで、軟材料の動的軸対称応力場に含まれるすべての成分を再構成することが可能となる。応力テンソルの全成分を測定できれば、二次主応力ではなく主応力およびその方向の解析が可能になるだけでなく、材料の降伏点評価に用いられるフォン・ミーゼス応力の測定も可能になる［23］。 本研究で扱う対象は、静的および動的なヘルツ接触問題である。その理由は二つある。第一に、外力や接触条件が既知であれば、材料内部の応力場や表面に作用する衝撃力は理論的に導出可能である［32–36］。これにより、実験的に得られた光弾性パラメータから再構成された応力場の確からしさを評価できる。第二に、動的な軸対称応力場は多くの現象に現れる。たとえば、液滴が軟表面に衝突する場合［37,38］、液体ジェットがゲルに浸透する場合［4,5,31］、ゲル内でキャビテーションが発生する場合［39］などである。これら二つの理由から、本研究は光弾性トモグラフィーが動的軸対称応力場を測定するための有効な手法であることを示す。我々は、光弾性トモグラフィーが軟材料の高速現象を理解する手法として広く応用されると信じている。

2 方法論

2.1 積分光弾性法（Integrated photoelasticity）

円偏光光源からの光が応力を受けた材料に入射すると、その光は楕円偏光となり、材料内の応力状態に対応した**光弾性パラメータ（位相遅れ Δ と偏光方向 ϕ）**をもって放出される［18］。Δ と ϕ は高速偏光カメラを用いることで時間的に測定可能である［30］。 材料内部の応力場と光弾性パラメータとの間には次の関係がある［40–42］： $$V_1 \equiv \Delta \cos 2\phi = C \int_{-\infty}^{\infty} (\sigma_{xx} - \sigma_{zz}) \, dy, \tag{2}$$ $$V_2 \equiv \Delta \sin 2\phi = 2C \int_{-\infty}^{\infty} \sigma_{xz} \, dy, \tag{3}$$   ここで C は材料固有の応力光学係数であり、σxx, σzz, σxz は直交座標系での応力成分（y 軸をカメラの光軸とする）である。Δ が λ/4 より小さく、主応力方向の回転が π/6 を超えない場合には式 (2)、(3) が成り立つ［12,19］。本研究の全実験データはこの仮定の範囲内である。 これらの積分式からわかるように、偏光カメラで測定される光弾性パラメータは、カメラ光軸方向に沿った内部応力場の積分値である。そのため、材料内部の各点での応力テンソル成分を得るには、光弾性パラメータを応力場に変換するトモグラフィーが必要となる。 ただし応力場のトモグラフィーは、各点がスカラー量で特徴づけられる通常のトモグラフィーとは異なる。応力はテンソルで表されるため、従来のスカラートモグラフィーをそのまま応用するのは難しい。 Aben らは、応力場トモグラフィーの問題を「ある一つの応力テンソル成分のスカラートモグラフィー」に帰着できれば解決できると提案した［18,24］。材料内部の応力場を光弾性トモグラフィーで決定するには、x–y 平面に平行で距離 Δz だけ離れた上下二つの断面で光弾性パラメータを測定する必要がある［28,43］。一般的には、z 軸周りに対してカメラ光軸の向きを多方向から測定することが望ましい。 任意の三次元応力場の一部要素（Fig.2 参照）について、x 軸方向の力の釣り合いは次式で表される： $$\Delta z \int_A^B \sigma_{xx} \, dy = T_u - T_l, \tag{4}$$   ここで Tu, Tl は要素の上下表面に作用するせん断力、A, B は y 軸方向の応力場境界を示す。これらは式 (3) を用いて記述できる： $$T_u = \frac{1}{2C} \int_{x}^{X'} V'_2 dx, \quad T_l = \frac{1}{2C} \int_{x}^{X} V_2 dx. \tag{5}$$   これを整理すると： $$\int_A^B \sigma_{xz} dy = \frac{1}{2C} V_2, \tag{6}$$ $$\int_A^B \sigma_{zz} dy = \frac{1}{2C\Delta z} \left( \int_x^{X'} V'_2 dx - \int_x^{X} V_2 dx \right) - \frac{V_1}{C}. \tag{7}$$   この式により、応力場のテンソルトモグラフィーをスカラートモグラフィーとして扱うことができる。なお、式 (7) による軸応力の再構成には力の釣り合い仮定が必要であるのに対し、式 (6) によるせん断応力の再構成にはその仮定は不要である。

---

2.2 光弾性トモグラフィー（Photoelastic tomography）

応力場が軸対称である場合（例：剛体球が平板に垂直接触する問題）、式 (6)、(7) を用いて実験的に応力場を決定できる［27,28］。このとき応力テンソルは以下の 4 成分を持つ：

• 半径方向応力 σrr
• 周方向応力 σθθ
• 軸方向応力 σzz
• せん断応力 σrz

すなわち式 (1) に対応する。 このとき、カメラ光軸に垂直な平面は r–z 平面、対称軸は z 軸となる。σrz および σzz は測定された Δ と ϕ から積分関係式を用いて再構成され、残りの σrr, σθθ は線形弾性理論の式を用いて再構成される［27］。 以下、この小節では一般的な軸対称スカラー場に対するトモグラフィーとそのアルゴリズムを説明し、その後に応力場への適用を述べる。

---

2.2.1 Abel変換とオニオン・ピリング法（onion-peeling method）

軸対称場を再構成する代表的手法が逆 Abel 変換である。 一般に、軸対称場 p(r) と、その積分から得られる投影関数 P(x) の関係は Abel 変換と呼ばれ、次式で表される： $$P(x) = \int_{-\infty}^{\infty} p\left(\sqrt{x^2 + y^2}\right) \, dy, \tag{8}$$   書き換えると： $$P(x) = 2 \int_x^{\infty} \frac{p(r) r}{\sqrt{r^2 - x^2}} \, dr. \tag{9}$$   その逆変換は次式で与えられる： $$p(r) = -\frac{1}{\pi} \int_r^{\infty} \frac{dP}{dx} \frac{dx}{\sqrt{x^2 - r^2}}. \tag{10}$$   しかし実験で得られる投影データは離散化されており、上式を直接適用すると x=r に特異点が生じるほか、x 方向の P の正確な変化が事前に既知でなければならない。さらに、Abel 変換は無限範囲を前提とするが、実験データは有限範囲であり、外側では p, P が厳密にゼロでなければならない。このため、対象は観測範囲内に完全に含まれる局所的な物体でなければ精度が落ちる。 この課題を解決するために用いられる数値アルゴリズムの一つがオニオン・ピリング（OP）法である［44,45］。OP 法では、軸対称場を厚さ Δr の N 層のリング（Fig.3）としてモデル化する。 最外層リングの値 p(rN) は、投影データ P(rN) と光路長 2WN,N を用いて次式で求められる： $$P(r_N) = 2 W_{N,N} p(r_N). \tag{11}$$   次に、N−1 層目の値 p(rN−1) は、投影データ P(rN−1)、既に求めた最外層 p(rN)、およびそれぞれの光路長を用いて次式で求められる： $$P(r_{N-1}) = 2 W_{N-1,N-1} p(r_{N-1}) + 2 W_{N-1,N} p(r_N). \tag{12}$$   一般に、i 層目の投影データは次式で表される： $$P(r_i) = 2 \sum_{j=i}^{N} W_{i,j} p(r_j), \tag{13}$$   ここで ri = iΔr は中心軸からの距離であり、Wi,j は次式で与えられる： $$W_{i,j} = \begin{cases} 0 & j<i, \\ \frac{\Delta r}{2} \left(\sqrt{(2j+i)^2 - 4i^2}\right) & j=i, \\ \frac{\Delta r}{2} \left(\sqrt{(2j+i)^2 - 4i^2} - \sqrt{(2j-i)^2 - 4i^2}\right) & j>i. \end{cases} \tag{14}$$  

2.2.2 せん断応力と軸方向応力の決定

せん断応力 σrz は、光弾性パラメータから式 (3) と OP（onion-peeling）法を基に決定する。前節で述べた OP 法は軸対称場を対象とするが、直交座標系 (x–y 平面) で分布する場合、せん断応力 σrz はもはや完全に軸対称ではない。このとき直交座標系におけるせん断応力 σxz と円筒座標系における σrz の関係は次式で表される： $$\sigma_{xz} = \sigma_{rz} \cos\theta \tag{15}$$   したがって、光弾性パラメータの積分値とせん断応力の関係を導く際には角度 θ を考慮する必要がある。これは、前節の単純な OP 法との大きな違いである。つまり、再構成すべき応力 σrz と積分される応力 σxz は角度 θ を介して結び付けられる［28］。これは、光弾性法がカメラ光軸に垂直な平面上に投影された応力成分しか測定できないためである。 N 番目のリングにおける V2(N) は、リングに作用するせん断応力 σrz(N) の積分であり、次式で表される： $$V^{(N)}_2 = 4 C W_{N,N} \sigma^{(N)}_{rz} \cos\theta_{N,N} \tag{16}$$   ここで WN,N は、カメラ光軸が N 番目のリングを通過する際の光路長である。この式より、リングに作用する σrz(N) を決定できる。N−1 番目のリングについても、既に決定された σrz(N) を用いて同様に計算できる： $$V^{(N-1)}_2 = 4C W_{N-1,N-1} \sigma^{(N-1)}_{rz}\cos\theta_{N-1,N-1} + 4C W_{N-1,N} \sigma^{(N)}_{rz}\cos\theta_{N-1,N} \tag{17}$$   同様に、i 番目のリングに対して一般化すると： $$V^{(i)}_2 = 4C \sum_{j=i}^N W_{i,j}\sigma^{(j)}_{rz}\cos\theta_{i,j} \tag{18}$$   ここで θi,j は次式で与えられる： $$\theta_{i,j} = \begin{cases} 0 & j=i \\ \cos^{-1}\left(\frac{2i}{2j-1}\right) & j>i \end{cases} \tag{19}$$   軸応力 σzz も同様に決定できる。式 (7) より、i 番目のリングに作用する σzz(i) は： $$V^{(i)}_1 - \frac{1}{2\Delta z}\sum_{j=1}^i \big(V'^{(j)}_2 - V^{(j)}_2\big)\Delta x = 2C \sum_{j=1}^i W_{i,j}\sigma^{(j)}_{zz} \tag{20}$$   文献［28］によれば、式 (18)、(20) を用いて外側のリングから内側へと再帰的に計算すれば、全リングにおけるせん断応力と軸応力を決定できる。 さらに計算効率を高めるため、これらの総和式を行列形式に変換し、係数行列と応力ベクトル・積分値ベクトルの積として表すことができる。たとえば式 (18) は： $$V_2 = 4C \alpha_{rz} \sigma_{rz} \tag{21}$$   同様に、軸応力についても： $$V_1 - \frac{\Delta x}{2\Delta z}\beta V_2 = 2C \alpha_{zz}\sigma_{zz} \tag{24}$$   このとき、逆行列 αrz⁻¹, αzz⁻¹ をあらかじめ計算しておけば： $$\sigma_{rz} = \frac{1}{4C}\alpha_{rz}^{-1}V_2, \quad \sigma_{zz} = \frac{1}{2C}\alpha_{zz}^{-1}\left(\frac{\Delta x}{2\Delta z}\beta V_2 - V_1\right) \tag{28,29}$$   これにより、外側から逐次的に計算する方法に比べて大幅に計算時間を短縮できる。

---

2.2.3 半径方向応力と周方向応力の決定

半径方向応力 σrr と周方向応力 σθθ は、線形弾性理論の方程式を用いて決定する［25,27］。これらは以下の釣り合い方程式および適合条件式を満たす必要がある： $$\frac{\partial\sigma_{rr}}{\partial r} + \frac{\sigma_{rr}-\sigma_{\theta\theta}}{r} + \frac{\partial\sigma_{rz}}{\partial z} = 0 \tag{30}$$ $$\frac{\partial}{\partial r}\{\sigma_{\theta\theta} - \nu(\sigma_{rr} + \sigma_{zz})\} - (1+\nu)\frac{\sigma_{rr}-\sigma_{\theta\theta}}{r} = 0 \tag{31}$$   ここで ν はポアソン比である。これを変形すると： $$\frac{\partial\sigma_{rr}}{\partial r} = -\frac{\sigma_{rr}-\sigma_{\theta\theta}}{r} - \frac{\partial\sigma_{rz}}{\partial z} \tag{32}$$ $$\frac{\partial\sigma_{\theta\theta}}{\partial r} = \frac{\sigma_{rr}-\sigma_{\theta\theta}}{r} - \nu \frac{\partial\sigma_{rz}}{\partial z} + \nu \frac{\partial\sigma_{zz}}{\partial r} \tag{33}$$   すでにせん断応力 σrz と軸応力 σzz の分布が 2.2.2 で得られているため、これらを既知項として数値的に差分解を行うことで σrr, σθθ を求めることができる。 一次前進差分を用いると： $$\sigma^{(i)}_{rr} = \sigma^{(i-1)}_{rr} + \Delta r \left( \ldots \right) \tag{34}$$ $$\sigma^{(i)}_{\theta\theta} = \sigma^{(i-1)}_{\theta\theta} + \Delta r \left( \ldots \right) \tag{35}$$   境界条件を適切に設定すれば、半径方向応力と周方向応力を決定できる。

---

2.2.4 準静的近似

これまでのトモグラフィーアルゴリズムは、せん断応力 σrz の再構成を除き、力の釣り合いを仮定している。したがって、測定された光弾性パラメータの時系列データから応力を再構築するには**準静的近似（quasi-static approximation）**を採用する必要がある。 この近似の妥当性を確認するために、衝突速度 V と材料内の圧力波速度 vp を比較する。圧力波速度 vp は次式で与えられる［32］： $$v_p = \sqrt{\frac{E}{\rho_s}} \tag{36}$$   ここで E は基板材料のヤング率、ρs は密度である。球の速度 V が圧力波速度 vp より十分に小さい低速衝突の場合、準静的近似は有効であると考えられる。

2.3 実験

静的および動的ケースにおける実験検証のための装置構成図をそれぞれ Fig.4(a)、Fig.4(b) に示す。実験装置は光源、測定対象物、および偏光カメラで構成される。 測定対象には、ウレタンゲルブロック（ウレタンゲルファントム、株式会社エクシール［岐阜、日本］）を使用した。サイズは 50 × 50 × 50 mm³、密度は ρs = 1064 kg/m³ である。ゲルのヤング率 E は、電子天秤を用いて測定した z 軸方向の表面変形量 δzmax と荷重垂直力 Fn から、ヘルツ接触理論 [32] を用いて推定した： $$\delta_{zmax} = \left(\frac{9F_n^2}{16E^{*2}R}\right)^{1/3}, \tag{37}$$   ここで $E^* = \frac{E}{1-\nu^2}$ である。本研究では、ゲルのポアソン比 ν を 0.5 と仮定し、ヤング率は 47.4 kPa と決定された。ヘルツ接触理論は接触半径（球と基板の接触領域）が R より十分小さいことを前提とするが [32,47]、実際の表面変形結果は式 (37) による曲線と良好に一致した（Fig.5）。このヘルツ接触問題に関する詳細は、著者らの先行研究 [48] に記載されている。

---

静的実験 (Fig.4(a))

直径 2R = 14.7 mm のスチロール球をゲル表面に垂直に押し当て、荷重垂直力 Fn を作用させた。この実験は「球と半無限体のヘルツ接触」[32] を模擬している。荷重垂直力 Fn は Fn = mg により算出され、m は電子天秤で測定した見かけ質量、g は重力加速度 9.81 m/s² である。 球を押し付けた際の位相遅れ（Δ）と偏光方向（ϕ）の分布は、偏光カメラ（Photron, CRYSTA PI-5WP）および光源（Thorlabs, SOLIS-565C）を用いて記録した。光源の波長帯域幅による測定精度の低下を抑えるため、カメラセンサ前に半値幅約 20 nm のバンドパスフィルタ（Photron製）を設置し、中心波長 λ = 540 nm の光を使用した。 測定画像の解像度は 191 × 208 ピクセルで、観測領域は 1.37R × 1.50R。これを用いて応力再構成を行った。空間分解能は 52.7 µm/pix である。

---

動的実験 (Fig.4(b))

半径 R = 2.98 mm、密度 ρ = 997.3 kg/m³ のプラスチック球を自由落下させ、ゲルに衝突させた。球の衝突速度 V は落下高さを調整することで 0.3～0.7 m/s の範囲に設定した。 衝突時の位相遅れおよび偏光方向分布は、高速偏光カメラ（Photron, CRYSTA PI-1P）と光源（Thorlabs, SOLIS-525C）を用いて記録した。カメラレンズ前にバンドパスフィルタを設置し、波長 λ = 520 nm の光を使用した。 カメラの時間分解能は 20,000 fps であり、高速応力測定が可能である。測定画像の解像度は 115 × 171 ピクセル、観測領域は 3.28R × 5.14R。空間分解能は 84.9 µm/pix である。

---

光学系とデータ処理

光源からの入射光（波長 λ）は、角度 0° の直線偏光子を通過した後、角度 45° の 1/4 波長板を通して円偏光となる。応力を受けたゲルは、この円偏光を透過させ、位相遅れ Δ と偏光方向 ϕ をもつ楕円偏光として放出する。偏光カメラを用いることで、透過光の Δ と ϕ を同時に測定できる。 偏光カメラのイメージセンサ上には、角度 0°、45°、90°、135° の 4 枚の線偏光子が隣接ピクセルごとに配置されており、これを「スーパー・ピクセル」[50] と呼ぶ。各偏光子を通して得られる強度をそれぞれ I0°, I45°, I90°, I135° とすると、4段階位相シフト法［15,16,30,51］により以下の式で Δ と ϕ を求められる： $$\Delta = \frac{\lambda}{2\pi} \sin^{-1} \frac{\sqrt{(I_{90°}-I_{0°})^2 + (I_{45°}-I_{135°})^2}}{I/2}, \quad \phi = \frac{1}{2}\tan^{-1}\frac{I_{90°}-I_{0°}}{I_{45°}-I_{135°}},$$   ここで $I = I_{0°} + I_{45°} + I_{90°} + I_{135°}$ である。 得られたデータの例を Fig.6(a), (b) に示す。これらは専用ソフトウェア（Photron Ltd., CRYSTA Stress Viewer）で処理され、さらに Matlab のメディアンフィルタ関数 (medfilt2) を用いてノイズ低減を行った。測定原理の詳細は文献 [16,30,52] にも記されている。

---

理論値との比較

理論的な位相遅れおよび偏光方向分布は、ヘルツ接触理論 [16,23,32,33] に基づく応力場の理論解を積分光弾性法 [16] に適用することで得られる。例を Fig.6(c), (d) に示す。 実験条件と同一の設定で計算を行い、応力光学係数 C = 1.14 × 10⁻⁹ Pa⁻¹ を用いた。この値は、計算された位相遅れと実験データとの差（RMSE）を最小化することで決定された（決定手順は著者らの先行研究 [16] に記載）。結果は Fig.7 に示す通りであり、応力光学係数は広い荷重範囲 Fn にわたってほぼ一定であった。

3 結果と考察

3.1 静的ヘルツ接触における光弾性トモグラフィーの実験的検証

本節では、剛体球を弾性半無限体に押し付ける静的ヘルツ接触問題を対象とし、大変形を示す軟材料において、光弾性トモグラフィーが軸対称応力場の全成分（σrr, σθθ, σzz, σrz）を再構成できることを示す。 Fig.8 では、左側に理論解（ヘルツ接触問題を m, R, E の条件で解いたもの）、右側に実験的な光弾性パラメータから再構成した応力場を示す。すべての応力成分において、ゲル表面近傍で誤差が見られるものの、全体的に両者は良好に一致している。これは、実験と理論計算で球とゲル表面の接触条件が必ずしも一致しないためであり、その差が誤差の原因であると考えられる。 特に、せん断応力 σrz の分布におけるノイズの原因として以下が挙げられる：

• z 軸近傍で偏光方向の値がゼロに近づき（Fig.6(b)参照）、負の値を取ることがある。理論的には z 軸右側の方位角はゼロ以上であるため、再構成前のデータフィルタリングや z 軸位置の決定がノイズに大きく影響する。
• 軸応力 σzz の分布には水平縞状のノイズが見られるが、せん断応力には見られない。これは σzz の再構成に V2 の z 軸方向勾配（V'2 − V2）が利用されるためであり、V2 の分布が完全に滑らかでないことから、不連続性が z 軸方向に残り、それが縞状ノイズを生じさせる。結果として、このノイズは σrr および σθθ にも波及する。

すべての応力成分が決定できたことから、延性材料の破壊基準や塑性変形の指標として用いられるミーゼス応力 σMも評価可能となる［23］： $$\sigma_M = \sqrt{3 J_2}, \tag{40}$$   ここで $$J_2 = \tfrac{1}{6}\big[(\sigma_{rr}-\sigma_{\theta\theta})^2 + (\sigma_{\theta\theta}-\sigma_{zz})^2 + (\sigma_{zz}-\sigma_{rr})^2\big] + \sigma_{r\theta}^2 + \sigma_{\theta z}^2 + \sigma_{rz}^2. \tag{41}$$   Fig.8(e) は理論値と実験値のミーゼス応力を比較したものであり、良好な一致を示している。 さらに Fig.9 では、z/R = 0.7 の位置における再構成応力場と理論解の線プロファイルを比較した。分布の形状はよく一致しており、ノイズは主に光弾性パラメータ（位相遅れ・偏光方向）のノイズに起因している。特に σzz の計算においては、式 (27) で示された積分値ベクトルが z 方向隣接画素間の変化量に敏感であるため、V2 のノイズが σzz に強く影響する。より適切なフィルタ処理を施すことで、このノイズは低減できると考えられる。 これらの結果は、光弾性トモグラフィーが大変形を伴う軟材料における軸対称応力場の全成分を再構成可能であることを示している。これにより、材料の降伏評価に重要なミーゼス応力の推定や、主応力およびその方向の解析が可能となる。

---

3.2 動的ヘルツ接触における高速応力場測定

本節では、剛体球が弾性半無限体に衝突する動的ヘルツ接触問題を対象とし、光弾性トモグラフィーを用いた軟材料の高速応力場測定の適用可能性を検討する。 まず、準静的近似の妥当性を確認するために、球の衝突速度 V と基板の圧力波速度 vp（Sec.2.2.4 参照）を比較した。実験では V の最大値は 2.7 m/s、vp ≃ 6.7 m/s（E = 47.4 kPa, ρs = 1064 kg/m³）であり、最大比 V/vp ≃ 0.4 であった。V/vp < 1 であるため、準静的近似が成立すると判断した。 Fig.10 は半径 R = 2.9 mm の球が V ≃ 2.2 m/s で衝突した際のゲル中の位相遅れ・偏光方向分布を示す。

• 衝突前 (t = 0 ms)：応力がないため位相遅れはゼロ、偏光方向は −90°～90°のランダム値。
• 衝突直後 (t ≃ 0.5–0.75 ms)：球の変位が最大となり、位相遅れも最大となる。
• その後：球の上方運動に伴い、位相遅れは減少。

Fig.11, Fig.12 は、Fig.10 のデータから再構成したせん断・軸・半径・周方向応力分布の時間変化を示す。最大応力は 0.25～0.5 ms の間に発生し、その後は球の反発運動に伴い減少した。接触端部近傍では、付着力により σrz, σzz が負値を示す領域も確認された。σrr では z 方向に広がった応力分布が観測され、負の応力領域も現れた。σθθ は衝突直後に正負両方の値を示し、その後負の応力が減衰した。 Fig.13 は基板表面 (z=0) におけるせん断応力 σrz(r,0,t) および軸応力 σzz(r,0,t) のキモグラフを示す。カラーバー範囲を変えて弱い応力も可視化すると、せん断応力と軸応力の波が異なる速度で伝播することがわかる。これは偏光カメラで得られる Δ・ϕ の可視化データだけでは得られず、応力場を再構成して初めて定量化可能となる。 さらに、O(10⁻¹)～O(10¹) kPa の広い範囲の応力場を同時に測定できることも示された。測定可能範囲は光源波長 λ および材料の応力光学係数 C に依存する。式 (38) より位相遅れ Δ は 0～λ/4 に制限され、式 (2)(3) より Δ の大きさは C によって決定される。したがって、使用するゲルを変更（C を変化）することで測定感度を調整可能である。たとえばゼラチンゲルは本研究のウレタンゲルよりも C が大きいため、小さい応力に対してより敏感である［12,16,53,54］。また、カラー偏光カメラを用いることで測定範囲をさらに拡張できる［55］。 これらの結果は、光弾性トモグラフィーが応力波伝播のような高速かつ時間変動する応力場においても有効に適用可能であることを示している。

3.3 最大衝撃力のスケーリング則

本節では、ゲル表面に作用する法線方向の衝撃力を推定し、その最大値に関するスケーリング則を動的ヘルツ接触理論からの理論予測と比較することで、光弾性トモグラフィーによる時間変化する応力場再構成の精度を評価する。 法線方向の衝撃力 Fn は、ゲル表面に作用する軸応力 σzz を積分することで求められる： $$F_n(t) = 2\pi \int_0^{\infty} \sigma_{zz} \, r \, dr. \tag{42}$$   Fig.14 は異なる衝突速度に対する Fn の時間発展を示す。衝撃後に Fn は急激に増加して最大値に達し、その後は球の鉛直運動に伴って減少する。これはゲルの粘性減衰効果による。衝突速度が高い場合、球が上昇に転じると Fn が負になることもあり、これは球とゲル表面間の付着力によると考えられる。 ここで、ピーク力 Fn,max とその発生時刻 tFn,max に注目する。ピーク力と時刻のスケーリング則は、動的ヘルツ接触モデル［35,36］から理論的に導かれる。球の運動方程式は次式で表される： $$m \frac{d^2 z}{dt^2} + F_n = 0, \tag{43}$$   初期条件は t=0 で z=0, dz/dt=V である。z は球の鉛直変位である。式 (43) にヘルツ接触理論の関係： $$F_n = \frac{4}{3} E^* R^{1/2} z^{3/2}, \tag{44}$$   を代入すると： $$\frac{4}{3}\pi \rho R^3 v \frac{dv}{dz} + \frac{4}{3}E^* R^{1/2} z^{3/2} = 0, \tag{45}$$   ここで v = dz/dt。これを解くと： $$v^2 = V^2 - \frac{4}{5}\frac{E^*}{\pi \rho R^{5/2}} z^{5/2}. \tag{46}$$   最大変位 zmax は v=0 のときに達し、そのスケーリングは： $$z_{max} \propto E^{-2/5} \rho^{2/5} R V^{4/5}. \tag{47}$$   これを式 (44) に代入すると、ピーク力のスケーリング則は： $$F_{n,max} \propto E^{2/5} \rho^{3/5} R^2 V^{6/5}. \tag{48}$$   さらに、Fn,max を慣性力 ρV²R² で無次元化すると： $$\frac{F_{n,max}}{\rho V^2 R^2} \propto \left(\frac{\rho V^2}{E}\right)^{2/5}. \tag{49}$$   Fig.15(a) は実験値と式 (49) の比較を示し、実験結果は理論予測とよく一致している。 一方、ピーク時刻 tFn,max については、zmax ∝ tFn,max V と式 (47) を用いることで： $$t_{F_n,max} \propto E^{-2/5} \rho^{2/5} R V^{-1/5}. \tag{50}$$   これを衝突時間 T=R/V で無次元化すると： $$\frac{t_{F_n,max}}{T} \propto \frac{\rho V^2}{E}^{2/5}. \tag{51}$$   Fig.15(b) は実験結果と式 (51) を比較したものであり、ピーク時刻についても理論予測と一致している。 なお、動的ヘルツ接触理論も含め、ヘルツ接触理論は基本的に基板の小変形を仮定している。しかしながら、今回の大変形を伴う実験においてもスケーリング則は理論と良好に一致した。

---

4 結論

本研究では、光弾性トモグラフィーを用いて軟材料の動的軸対称応力場の全成分を再構成した。静的および動的ヘルツ接触に焦点を当て、手法の有効性を定量的に検証した。

• 静的ヘルツ接触：剛体球をゲルに押し付け、偏光カメラで測定した光弾性パラメータから全応力テンソル成分を再構成した。その結果をヘルツ接触理論による応力場と比較したところ、両者は良好に一致し、光弾性トモグラフィーが軟材料の応力場測定に適用可能であることが示された。さらに全応力成分を得られるため、材料の降伏値評価に用いられるミーゼス応力などの解析も可能となった。
• 動的ヘルツ接触：剛体球をゲルに衝突させた場合について、偏光カメラで測定した光弾性パラメータから動的応力場を再構成した。その結果、O(10⁻¹)～O(10¹) kPa に及ぶ広範囲の応力を大変形下で同時に測定できた。さらに、応力波伝播のように時間的に急激に変化する応力場の測定にも成功した。また、再構成した応力場から表面に作用する法線衝撃力を算出し、その最大値に関するスケーリング則を理論予測と比較したところ、良好な一致を示した。

これらの結果は、高速光弾性トモグラフィーが軟材料における動的応力場の高精度測定に有効であることを確認するものである。本手法は、剛体球の衝突だけでなく、液滴や液体ジェットの衝突など多様な工学的プロセスにおける応力分布の理解に寄与する可能性がある。

---

謝辞

本研究は以下の助成を受けて実施された：

• JSPS 科研費 (JP20H00223, JP22J13343, JP22KJ1239, JP23KJ0859)
• JST PRESTO (JPMJPR21O5, Japan)