偏微分方程式

偏微分方程式 傾きがわかった関数の軌跡を求める．変数は二つ以上

微分方程式の種類 • 常微分方程式→独立変数が一個単振動 m:mass, c:dash pot, k:spring, t:time, x:axis • 偏微分方程式→独立変数が２個以上　　　　　　　　　　　　　　　　　　　一次元熱伝達 T:temperature, t:time, x:axis 二階の常微分方程式 d2x dx m + c + kt = 0 dt2 dt 二階の偏微分方程式 ∂T ∂2T cr = l ∂t ∂x2

微分方程式を解くという意味 dy dx わかっている各点の傾きが既知 y この線を求める初期値初期値 x

変数が２個の場合 ∂u ∂x ∂u ∂y ，　　　　わかっている u 初期値 y 傾きが分っているこの面を求める x

変数が２個の場合 ∂u ∂x ∂u ∂y ，　　　　わかっている u 初期値 y Dx 次の点を求める方法 x

変数が２個の場合 ui - ui-1 = Δx ∂u ∂x u 次の点 x = xi x = xi-1 前の点 y yj-2 yj-1 yj yj+1 yj+2 yj+3

二つの変数の例（偏微分方程式の前に） h1, h2, h3, h4, h5の変化？初めは，すべて h1～h5 = 1.0 時間に対する高さの変化？ h1 h2 h3 h4 h5 高さh Dt 時間の増加量 Dt(h2-h3)-Dt(h3-h4) 断面積 A = 0.5 単位時間当たりの流量 F = h1 – h2

Ｅｘｃｅｌで計算 初期値流出量/断面積流入量/A 流出量/A 時間刻みDt

Ｅｘｃｅｌで確認 Dt Dt A 前の値ー流出前の値ー流出＋流入 F = h1 – h2

計算結果 1秒後 2秒後

二つの変数の例（偏微分方程式の前に） 一秒間に１ｍ３ h1, h2, h3, h4, h5の変化？初めは，すべて h1～h5 = 0.0 十分時間がたったときの h1～h5　？ h1 h2 h3 h4 h5 高さh 断面積 A = 0.5 単位時間当たりの流量 F = h1 – h2

Ｅｘｃｅｌで計算 初期値時間刻み当たりの流出量前回と同じ

Ｅｘｃｅｌで計算 初期値時間刻み当たりの流出量/断面積前回と同じ

熱伝導の基礎 ∂T ∂x 温度T x ∂T ∂x x Q = - l A l:熱伝導率 A:断面積伝わる熱量は温度勾配に比例

熱伝導 ∂T ∂t ∂T ∂t ∂T ∂x ∂T ∂x ∂ ∂x ∂T ∂x Q = - l A x x+ Δx x l cr ∂2T ∂x2 l:熱伝導率 c:比熱 r:密度 A:断面積 Q - Q’ = l ADx = cr ADx Q’ = - l A + (- l A)Dx ∂2T ∂x2 = xに関する二階微分

の式ここの温度Tj x x+Δx x T t = ti t = ti-1 前の点ここの２階微分 Tj+1- Tj Tj- Tj-1 - h h h h x xj-2 xj-1 xj xj+1 xj+2 xj+3 ∂2T ∂x2 ∂2T ∂x2 iは時間変化 jは座標変化 Tj+1- 2Tj+ Tj-1 = h2

Tiの式 ∂T ∂t ∂T ∂t Ti- Ti-1 = Dt l cr l cr Tj+1- 2Tj+ Tj-1 h2 Ti= Ti-1 + Dt ∂2T ∂x2 ∂2T ∂x2 = Tj+1- 2Tj+ Tj-1 = h2

Tiの式 Ti T t = ti l cr Tj+1- 2Tj+ Tj-1 h2 t = ti-1 Ti= Ti-1 + Dt Ti-1 前の点 x xj-2 xj-1 xj xj+1 xj+2 xj+3

テキスト問題 t = 0 1 温度u l cr 0 w 1 l:熱伝導率 A:断面積 = 1 両端の温度は常にゼロ

#include <stdio.h> #define N 20 int main( ) { double u[N+1], w[N+1]; double k=0.001; double h, r, s; int i, j; h = 1.0/(double)N; r = k/(h*h); s = 1.0 - 2.0*r; 今の温度次の温度時間刻み w[j] w, u t = ti u[j+1] u[j] ｘの区間 t = ti-1 u[j-1] x xj-2 xj-1 xj xj+1 xj+2 xj+3

for(i=0;i< N;i++) u[i]=1.0; for(i=0;i<=N;i++) w[i]=0.0; u[0] = 0.0; u[N] = 0.0; for(j=1;j<=200;j++){ if((j%10)==0){ printf("%5.3lf ", (double)j*k); for(i=0;i<=N;i+=2) printf("%5.3lf ", u[i]); printf("\n"); } for(i=1;i<N;i++) w[i] = r*(u[i+1]+u[i-1]) + s*u[i]; for(i=1;i<=N;i++) u[i] = w[i]; } return 0; } 初期条件温度＝１（両端以外）温度＝０（両端）剰余演算子１０で割ったあまり出力 l cr Tj+1- 2Tj+ Tj-1 h2 Ti= Ti-1 + Dt r = k/(h*h); s = 1.0 - 2.0*r; 1

for(i=1;i<N;i++) w[i] = r*(u[i+1]+u[i-1]) + s*u[i]; for(i=1;i<=N;i++) u[i] = w[i]; h = 1.0/N r = k/(h*h) s = 1.0 - 2.0*r = 1 - 2*k/h2 Tj+1- 2Tj+ Tj-1 h2 Ti= Ti-1 + Dt

微分方程式を考えない方法 t = 0 1 温度u l cr 0 w 1 l:熱伝導率 A:断面積 = 1 両端の温度は常にゼロ

微分方程式を考えない方法 t = 0 l:熱伝導率 A:断面積 r:密度 c:比熱 1 温度u 分割する h h 0 w 1 u[i-1] u[i] u[i+1] (u[i] - u[i-1])/h 温度勾配？ w[i]

熱伝導の基礎 ∂T ∂x 温度T x ∂T ∂x x Q = - l A l:熱伝導率 A:断面積伝わる熱量は温度勾配に比例温度勾配

微分方程式を考えない方法 l:熱伝導率 A:断面積 r:密度 c:比熱 k:時間刻み t = 0 1 温度u 分割する h h 0 w 1 u[i-1] u[i] u[i+1] w[i] Q = - (u[i] - u[i-1])/h*lAk Q’ = - (u[i+1] - u[i])/h*lAk w[i] = u[i] +(Q-Q’)/(hArc) = u[i] +(u[i+1]-2u[i]+u[i-1])lk/rch2

微分方程式を考えない方法 l:熱伝導率 A:断面積 r:密度 c:比熱 t = 0 1 温度u 分割する h h 0 w 1 u[i-1] u[i] u[i+1] w[i] Q = - (u[i] - u[i-1])/h*lAk Q’ = - (u[i+1] - u[i])/h*lAk w[i] = u[i] +(Q-Q’)/(hArc) = u[i] +(u[i+1]-2u[i]+u[i-1])lk/rch2

エクセルでトライ！ 時間刻み B2,C2,D3 から導出

エクセルでトライ！ 時間刻み B2,C2,D2 から導出 w[i] = u[i] +(u[i+1]-2u[i]+u[i-1])lk/rch2 lk/rch2=k*(l/rc)/h2=0.001*1/0.05^2=0.4

今日の課題弾性波の伝播（鉄の音速） u:変位 E:弾性係数 r:密度 A:断面積 x 衝撃を与える L u(x,t) x x

弾性波の伝播 u:変位 E:弾性係数 r:密度 A:断面積 x Dx 速度を与える L 質量：rDxA e x Dx ∂ (EeA)Dx ∂x EeA EeA +

弾性波の伝播 質量：rDxA e x Dx ∂ (EeA)Dx ∂x EeA EeA+ ∂2u ∂ F = rDxA= (EeA)Dx ∂t2 ∂x ∂u e = ∂x • ∂2u ∂2u • = E • ∂t2 ∂x2

弾性波の伝播 E:弾性係数 r:質量 A:断面積速度を与える L バネ定数：k = EA/h n+1個 k h F=kx λ=Ph/EA P=(EA/h)λ 質量：m = rAh

Kumamoto University 質点ばねモデル u:変位 E:弾性係数 r:密度 A:断面積 x 衝撃を与える L k1 m1 m2 k2 衝撃 v1 v2

n+1個 k u[i-1] u[i] u[i+1] 変位今の状態 v[i-1] v[i] v[i+1] 速度 w[i-1] w[i] w[i+1] 変位次の状態 s[i-1] s[i] s[i+1] 速度

n+1個 質量：m k L/n u[i-1] u[i] u[i+1] k*(u[i+1] - u[i]) k*(u[i] - u[i-1]) F = m*(s[i] - v[i])/t = k*(u[i+1] - 2*u[i] + u[i-1]) 次の速度今の速度時間刻み今の変位 s[i] = v[i] + t*k/m*(u[i+1] - 2*u[i] + u[i-1]) w[i] = u[i] + t*v[i]

エクセルでトライ！ t < h/(音速) 0.00001 二個セット m = rAh k = EA/h 二つ空ける時間刻み今の変位今の速度次の速度 s[i] = v[i] + t*k/m*(u[i+1] - 2*u[i] + u[i-1]) w[i] = u[i] + t*v[i] 次の位置 Fe: E=206×109, r=7800

エクセルでトライ！ 時刻に対する X20の変化をプロット時刻に対する v20の変化をプロット v20 ココの逆数が音速 t

分割の注意点 クーラン条件（時間刻み幅）時間の刻み幅Dtは，次の時間ステップで弾性波が次の要素に届かない時間 V*Dt < h 要素の寸法は最少波長の三分の一以下 h < λ/3

プログラムを作る！ 20m/s 衝撃時間　t = 0.001 s　（この時間 v=20とおく） L L= 1 m r = 7800 kg/m3 E = 206 GPa 鉄の音速を出すプログラム

#include <stdio.h> #define N 20 //領域分割 #define M 400 //時間分割 int main( ) { FILE *fp = fopen("data.csv","w"); double u[N+1], w[N+1]; //現変位，次変位 double v[N+1], s[N+1]; //現速度，次速度 double t=0.0000005; //時間刻み double E = 206e+9; //ヤング率 double r = 7800; //密度 double h, k; //区間幅，バネ定数 double L = 1.0; //長さ double ct = 0.001; //接触時間 double v0 = 20; //衝撃速度 double m; //質量 int i, j; V*Dt < h=1.0/20=0.05 V*Dt=5900*0.0000005=0.029

h = L/(double)N; k = E/h; m = h*r; for(i=0;i<=N;i++) { u[i]=0.0; w[i]=0.0; v[i]=0.0; s[i]=0.0; }

for(j=0;j<M;j++) { for(i=0;i<=N;i++) { if((double)j*t<ct) v[0] = v0; w[i]=u[i]+t*v[i]; if(i==0) s[i]=v[i]+k*t/m*(u[i+1]-u[i]); if(i>0 && i<N) s[i]=v[i]+k*t/m*(u[i+1]-2.0*u[i]+u[i-1]); if(i==N) s[i]=v[i]+k*t/m*(-u[i]+u[i-1]); } 質量：m i = N k i = 0 L/n

fprintf(fp, "%6.3lf ", (double)j*t*1000); for(i=0;i<=N;i+=2) fprintf(fp,"%7.3lf ", v[i]); fprintf(fp,"\n"); for(i=0;i<=N;i++) { u[i] = w[i]; v[i] = s[i]; } } fclose(fp); return 0; } 2回に1回

偏微分方程式

偏微分方程式

Presentation Transcript