第5章 連続体の力学 ５．１　連続体と偏微分方程式 【リスト5-1】カラーマップ表示関数 def chkC(R): #■R/G/B値のチェック if R>255: return 255 if R<0: return 0 return int(R) def dspCell(canvas, U,KX,KY,CS): #■結果表示 UMax = U[1][1]; UMin = U[1][1] #最大・最小値を求める for i in range(KX): UMax=max(UMax, max(U[i])) UMin=min(UMax, min(U[i])) S="最大値 %8.4f 最低値 %8.4f" % (UMax, UMin)#最大・最小表示 canvas.create_text(110,10, text=S, font=("ＭＳ　ゴシック",10),fill="blue") UDX= UMax - UMin; B=0 for j in range(1,KY): jj=j*CS+20 for i in range(1,KX): ii=i*CS+10 V=int(((U[i][j]-UMin)/UDX)*765)#表示色の決定 if V<256: G=int(55+200*V/255); R=0 elif V<512: V-=256; R=int(V) ; G=255 else : V-=512; R=255 ; G=255-int(V) C="#%02X%02X%02X" % (chkC(R),chkC(G),B) canvas.create_rectangle(ii,jj,ii+CS,jj+CS, fill=C, outline="") ５．２　定常状態 【リスト5-2】ラプラス方程式とポアソン方程式(カラーマップ表示については【リスト5-1】参照) from tkinter import * import random, math, time def array(N1, N2=0, N3=0): if N2 == 0 and N3==0: return [0 for i in range(N1)] elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)] else: return [[[0 for i in range(N3)] for j in range(N2)] for k in range(N1)] KX=100; KY=100; EPS=0.01; NN=1000 U=array(2,KX+1,KY+1) F=array(KX+1,KY+1) poisson=True def initCond(KX, KY, MX, MY, DX, DY):#■初期条件の設定 global F, poisson; CX=KX/2; CY=KY/2 for j in range(KY): for i in range(KX): F[i][j]=0 if poisson: F[i][j] = 100 * (DX * (MX - i) + DY * (MY - j)) def boundaryCond(KX, KY, DX, DY): #■境界条件設定 global U for k in range(2): for i in range(KX): U[k][i][0] = 0; U[k][i][KY] =16 for j in range(KY): U[k][0][j] = 0; U[k][KX][j] =4 def calc(KX, KY, NN, EPS): #■計算 global root, ID1, ID2 MX = KX + 1; MY = KY + 1; DX = 1.0/ KX; DY = 1.0 / KY F1 = 1 / (DX * DX); F2 = 1 / (DY * DY); F3 = 0.5 / (F1 + F2) initCond(KX, KY, MX, MY, DX, DY); boundaryCond(KX, KY, DX, DY) # 収れん計算 ID1 = 1; ID2 = 0 for N in range(NN): ID = ID1; ID1 = ID2; ID2 = ID; ER = 0 for j in range(1,KY): for i in range(1,KX): U[ID2][i][j]= F3*(F1*(U[ID1][i+1][j]+U[ID1][i][j+1])+ F2*(U[ID1][i-1][j]+U[ID1][i][j-1])+ F[i][j]) if abs(U[ID2][i][j])>EPS: E = abs((U[ID2][i][j]-U[ID1][i][j])/U[ID2][i][j]) if E > ER : ER = E if ER < EPS : return #収れんしたとき終了 if (N % 50)==0: root.title("%d %10.6f" % (N, ER)) root.update() ・ ・(中略)カラーマップ表示(chkC，dspCell)は【リスト5-1】参照 ・ #■メイン処理 root=Tk(); root.title("ポアソン方程式/ラプラス方程式") canvas=Canvas(root, width=220,height=230, bg="white");canvas.pack() calc(KX, KY, NN, EPS) if poisson: root.title("ポアソン方程式") else : root.title("ラプラス方程式") dspCell(canvas, U[ID2],KX,KY, 2); root.update() ５．３　熱拡散方程式 【リスト5-3】熱拡散方程式(カラーマップ表示については【リスト5-1】参照) from tkinter import * import random, math, time def array(N1, N2=0, N3=0): if N2 == 0 and N3==0: return [0 for i in range(N1)] elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)] else: return [[[0 for i in range(N3)] for j in range(N2)] for k in range(N1)] KX=50; KY=50; DT=0.0001;Ndisp=1;MX=KX+1; MY=KY+1 T=array(2,MX,MY); U=array(MX,MY) ; V=array(MX,MY) DX = 1 / KX; DY = 1 / KY R1 = 2*DT/DX; R2 = 2*DT/DY;R3 = DT/(DX*DX); R4 = DT/(DY*DY) N=0; ID1=0; ID2=1 for j in range(KY):#■初期条件の設定 for i in range(KX): YY=DY*j; U[i][j]=50*YY*(1-YY); V[i][j]=0 def boundaryCond(): #■境界条件設定 global T,KX,KY,MX,MY for j in range(MY): T[ID1][0][j]=0; T[ID1][KX][j]=T[ID1][KX-1][j] #左側冷却 for i in range(MX): T[ID1][i][0]=0; T[ID1][i][KY]=0 #上・下冷却 for i in range(int(MX/8-1), int(MX / 8 + 1)): for j in range(int(MY/8-1), int(MY / 8 + 1)): T[ID1][i][j] = 1 #上部中央部に熱源あり def calc(): #■計算 global root, T, ID1, ID2,R1,R2,R3,R4, KX,KY for j in range(1,KY): for i in range(1,KX): T[ID2][i][j]=(T[ID1][i][j] -R1*U[i][j]*(T[ID1][i+1][j]-T[ID1][i-1][j]) -R2*V[i][j]*(T[ID1][i][j+1]-T[ID1][i][j-1]) +R3*(T[ID1][i+1][j]-2*T[ID1][i][j]+T[ID1][i-1][j]) +R4*(T[ID1][i][j+1]-2*T[ID1][i][j]+T[ID1][i][j-1])) ID2 =ID1; ID1 = 1-ID1 ・ ・(中略)カラーマップ表示(chkC，dspCell)は【リスト5-1】参照 ・ def simulate(): global canvas, root,N,Ndisp for k in range(Ndisp):N+=1; boundaryCond();calc() canvas.delete(ALL) dspCell(canvas, T[ID2],KX,KY, 4) root.title("計算回数 %d " % N); root.update() #■画面コピー用マウスイベント def leftMouseDown(event): simulate() #■メイン処理 root=Tk(); root.title("ポアソン方程式/ラプラス方程式") canvas=Canvas(root, width=220,height=230, bg="white");canvas.pack() N=0; simulate() #左ボタン #canvas.bind("",leftMouseDown) #画面コピー用マウス #root.mainloop() #画面コピーするときは以下をコメントにする while True: simulate() ５．４　障害物がある場合の流体 【リスト5-4】障害物指定の処理 #■障害物の設定 flowExist=[[1 for j in range(NY)] for i in range(NX)] # 障害物の部分=0 def setObstacle(X1, Y1, X2, Y2):#障害物の設定 global flowExist for i in range(X1,X2+1): for j in range(Y1, Y2+1): flowExist[i][j]=0 #■障害物位置 setObstacle( 8, 12, 10, 16); setObstacle(16, 17, 20, 18) setObstacle(18, 18, 22, 22); setObstacle(20, 12, 24, 14) setObstacle(30, 16, 32, 18); setObstacle(44, 7, 46, 9) setObstacle(42, 13, 52, 15); setObstacle(38, 20, 44, 20) setObstacle(44, 21, 48, 25); setObstacle(38, 26, 44, 26) setObstacle(30, 8, 32, 12) 　 【リスト5-5】圧力境界条件の初期化 #■圧力境界条件初期化 boundLine=array(NX,NY) # 境界線 boundType=array(NX,NY) # 境界種類 for j in range(1,NY):#境界線の抽出 for i in range(1,NumX): if (flowExist[i+1][j]-flowExist[i][j])==-1:boundLine[i][j]=1 for i in range(NumX,1,-1): if (flowExist[i][j]-flowExist[i-1][j])== 1:boundLine[i][j]=1 for i in range(NX): for j in range(1,NumY): if (flowExist[i][j+1]-flowExist[i][j])==-1:boundLine[i][j]=1 for j in range(NumY,1,-1): if (flowExist[i][j]-flowExist[i][j-1])== 1:boundLine[i][j]=1 for i in range(1, NumX): for j in range(1,NumY): boundLine[i][j] = boundLine[i][j] * flowExist[i][j] for j in range(1,NumY):#パターンよる圧力境界条件の設定 for i in range(1,NumX): QQ = (boundLine[i+1][j]+boundLine[i-1][j]+ boundLine[i][j+1]+boundLine[i][j-1]) QQ *= (1-flowExist[i][j]); QQ=int(QQ) if QQ == 1: boundType[i][j] = 12 elif QQ == 2: boundType[i][j] = 6 elif QQ == 3: boundType[i][j] = 4 elif QQ == 4: boundType[i][j] = 3 【リスト5-6】各種パラメータの設定と初期化 from tkinter import * import random, math, time def array(N1, N2=0, N3=0): if N2 == 0 and N3==0: return [0 for i in range(N1)] elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)] else: return [[[0 for i in range(N3)] for j in range(N2)] for k in range(N1)] dspType=5 # 表示データ 1: 速度絶対値 2:X方向圧力　3:Y方向圧力 # 上記に4を加算すると速度(矢印)を重ね表示 cndDt=0.005 # 時間刻み cndRE=40 # レイノルズ数 cndLoop=1000 # 実行回数 pbIter=200 # 収れん繰返し回数 NumX=61; NumY=35 # セル数 SizeX=10; SizeY=10 # セルサイズ NX = NumX+1; NY=NumY + 1 # セル数+1 #■座標設定 pbX=array(NX); pbY=array(NY);pbX[1]=0; pbY[1]=0 for i in range(1, NX): pbX[i] = pbX[i-1]+ SizeX for i in range(1, NY): pbY[i] = pbY[i-1]+ SizeY #■計算因子初期化 pbAbsV=array(NX,NY)#(速度絶対値格納用) pbXG=array(NX); pbXA=array(NX); pbXB=array(NX) pbXC=array(NX); pbXD=array(NX) pbYG=array(NY); pbYA=array(NY); pbYB=array(NY) pbYC=array(NY); pbYD=array(NY) for i in range(1,NumX):#X方向計算用 pbXG[i] = 2/(pbX[i+1]-pbX[i-1]); pbXA[i] = pbXG[i]*pbXG[i] pbXB[i] = (pbX[i+1]-2*pbX[i]+pbX[i-1])*math.pow(pbXG[i],3) pbXC[i] = pbXA[i]-pbXB[i]*0.5; pbXD[i] = pbXA[i]+pbXB[i]*0.5 for i in range(1,NumY):#Y方向計算用 pbYG[i] = 2/(pbY[i+1]-pbY[i-1]); pbYA[i] = pbYG[i]*pbYG[i] pbYB[i] = (pbY[i+1]-2*pbY[i]+pbY[i-1])*math.pow(pbYG[i],3) pbYC[i] = pbYA[i]-pbYB[i]*0.5; pbYD[i] = pbYA[i]+pbYB[i]*0.5 #■初期条件 pbU=array(NX, NY);pbV=array(NX, NY); pbP=array(NX, NY) pbQ=array(NX, NY);pbD=array(NX, NY); pbS=array(NX, NY); pbUL=array(NX) for j in range(1,NY): for i in range(1,NX):pbU[i][j]=1 　 【リスト5-7】障害物がある場合の流体 #■圧力境界条件初期化 ・ ・（前略）【リスト5-6】,【リスト5-4】，【リスト5-5】の順に続ける。 ・ def boundVel():#■速度境界条件 global NumX, NumY, pbU, pbV J1=NumY; J2=J1-1; J3=J1-2; I1=NumX; I2=I1-1; I3=I1-2 for i in range(1,I1): pbU[i][J1] = pbU[i][J3]; pbV[i][J1] = -pbV[i][J3] pbU[i][J2] = pbU[i][J3]; pbV[i][J2] = 0 pbU[i][ 1] = pbU[i][ 3]; pbV[i][ 1] = -pbV[i][ 3] pbU[i][ 2] = pbU[i][ 3]; pbV[i][ 2] = 0 for j in range(1,J1): pbU[I1][j] = pbU[I3][j]; pbV[I1][j] = 0 　　　　pbU[1][j]=pbU[3][j];pbV[1][j] = 0 pbU[I2][j] = pbU[I3][j]; pbV[I2][j] = 0 　　　　pbU[2][j]=pbU[3][j];pbV[2][j] = 0; def tempVel():#■仮の速度 global NumX,NumY,cndRE, pbU, pbV,cndDT,pbXA, pbXB,cndDt global pbYA, pbYB, pbD, pbS,pbXG, pbYG, flowExist NY2 = NumY - 2; NX2 = NumX-2; R1 = 1 / cndRE for j in range(2,NY2+1): for i in range(2,NX2+1): U0 = pbU[i ][j ]; V0 = pbV[i ][j ] UP0 = pbU[i+1][j ]; UM0 = pbU[i-1][j ] U0P = pbU[i ][j+1]; U0M = pbU[i ][j-1] VP0 = pbV[i+1][j ]; VM0 = pbV[i-1][j ] V0P = pbV[i ][j+1]; V0M = pbV[i ][j-1] UNX = U0*(UP0-UM0)/2 - abs(U0)*(UP0-2*U0+UM0)/2 UNY = V0*(U0P-U0M)/2 - abs(V0)*(U0P-2*U0+U0M)/2 VNX = U0*(VP0-VM0)/2 - abs(U0)*(VP0-2*V0+VM0)/2 VNY = V0*(V0P-V0M)/2 - abs(V0)*(V0P-2*V0+V0M)/2 UV = ((UP0-2*U0+UM0)*pbXA[i] + (U0P-2*U0+U0M)*pbYA[j] - (UP0-UM0)*0.5*pbXB[i] - (U0P-U0M)*0.5*pbYB[j]) VV = ((VP0-2*V0+VM0)*pbXA[i] + (V0P-2*V0+V0M)*pbYA[j] - (VP0-VM0)*0.5*pbXB[i] - (V0P-V0M)*0.5*pbYB[j]) pbD[i][j] = U0 + cndDt*(-(UNX*pbXG[i]+UNY*pbYG[j])+R1*UV) pbS[i][j] = V0 + cndDt*(-(VNX*pbXG[i]+VNY*pbYG[j])+R1*VV) for j in range(2,NY2+1): for i in range(2,NX2+1): pbU[i][j] = pbD[i][j]*flowExist[i][j] pbV[i][j] = pbS[i][j]*flowExist[i][j] def PoissonRight():#■ Poisson右辺 global NumX,NumY, pbU, pbU, pbXG, pbYG, pbXA, pbYA, pbD,cndDt for j in range(1,NumY): for i in range(1,NumX): U1 = (pbU[i+1][j ]-pbU[i-1][j ])*0.5*pbXG[i] V2 = (pbV[i ][j+1]-pbV[i ][j-1])*0.5*pbYG[j] pbQ[i][j] = (U1 + V2) / cndDt for j in range(1,NumY): for i in range(1,NumX): pbD[i][j] = 0.5/(pbXA[i]+pbYA[j]) def compPress():#■圧力変化の計算 global NumX,NumY,boundLine, flowExist, boundType for j in range(1,NumY): for i in range(1,NumX): pbS[i][j]=(boundLine[i][j] * (pbP[i+1][j]+pbP[i-1][j]+pbP[i][j+1]+pbP[i][j-1])) for j in range(1,NumY):#障害物の考慮 for i in range(1,NumX): pbP[i][j]=(pbS[i][j]*(1-flowExist[i][j])*boundType[i][j]/12 +pbP[i][j]*flowExist[i][j]) def boundPress():#■圧力境界条件 global NumX,NumY, pbP I2 = NumX-1; I3 = NumX-2; J2 = NumY-1; J3 = NumY-2 for j in range(1,NumY): pbP[1][j] = pbP[2][j]; pbP[I2][j]= pbP[I3][j] for i in range(1,NumX): pbP[i][2] = pbP[i][3]; pbP[i][J2] = pbP[i][J3] def convCond():#■収束条件 global NumX, NumY, pbP, pbXC, pbYC, pbYC,pbYD global pbQ, pbUL, flowExist NY2 = NumY - 1; NX2 = NumX - 1 G1 = 0 for j in range(1, NumY): for i in range(1, NumX): pbUL[i] = (pbD[i][j]* (pbXC[i]*pbP[i+1][j]+pbYC[j]*pbP[i][j+1] + pbXD[i]*pbP[i-1][j]+pbYD[j]*pbP[i][j-1] - pbQ[i][j])-pbP[i][j]) for i in range(1, NumX): G1 += pbUL[i]*pbUL[i]*flowExist[i][j] for i in range(1, NumX): pbP[i][j] = pbUL[i] + pbP[i][j]*flowExist[i][j] return G1 def compPressCond(): #■圧力収束条件の計算 compPress(); boundPress() return convCond() def updateVel():#■仮の速度を圧力により変更 global NumX,NumY,cndDT, pbD, pbU, pbV, pbP,pbQ, pbXG, pbYG, cndDt global flowExist NX2 = NumX - 1; NY2 = NumY - 2; DT = cndDt for j in range(1,NumY): for i in range(1,NumX): pbD[i][j] = pbU[i][j]-DT*(pbP[i+1][j]-pbP[i-1][j])*0.5*pbXG[i] pbQ[i][j] = pbV[i][j]-DT*(pbP[i][j+1]-pbP[i][j-1])*0.5*pbYG[j] for j in range(1,NumY): for i in range(1,NumX): pbU[i][j] = pbD[i][j]*flowExist[i][j] pbV[i][j] = pbQ[i][j]*flowExist[i][j] GID=array(NumX+1,NumY+1) def initDsp():#■初期表示(各セルの四角形の図形識別番号を保存) global GID,NumX, NumY,canvas; CL="#555555" for i in range(NumX+1): ii=i*10 for j in range(NumY+1): jj=j*10; C="#FFFFFF" if i==0 or j==0 or i==NumX or j==NumY: C="#777777" GID[i][j]=canvas.create_rectangle(ii,jj, ii+20,jj+20,fill=C,outline=CL) def dspObstacle():#■障害物の表示(セルの色を変える) global GID,NumX, NumY,canvas for i in range(1, NumX): for j in range(1,NumY): #■流体と障害物の境界確認用 #if boundLine[i][j]== 1: # canvas.itemconfig(GID[i][j], fill="#FF0000") if flowExist[i][j]==0: C="#777777" #■障害物の圧力境界条件確認用 #if boundType[i][j]== 3: C="#FF0000" #elif boundType[i][j]== 4: C="#FFFF00" #elif boundType[i][j]== 6: C="#00FF00" #elif boundType[i][j]==12: C="#0000FF" #canvas.itemconfig(GID[i][j], fill=C) for j in range(1,NumY):#最左端と最右端は壁ではないので色を変える canvas.itemconfig(GID[0][j], fill="#AAAAFF") canvas.itemconfig(GID[NumX][j], fill="#AAAAFF") def chkC(R): #■R/G/B値のチェック if R>255: return 255 if R<0: return 0 return int(R) def dspCell(canvas, pbP, KX, KY): #■セルの色表示 global GID UMax = pbP[1][1]; UMin = pbP[1][1] #最大・最小値を求める for i in range(1, KX): for j in range(1,KY): if UMaxpbP[i][j]: UMin=pbP[i][j] UDX= UMax - UMin; B=0 if UDX<=0.0001: UDX=1 #最大最小の差がなければ一定値を基準とする for i in range(1, NumX): for j in range(1, NumY): if flowExist[i][j]!=0:#流体がある部分だけを色を変化させる V=int(((pbP[i][j]-UMin)/UDX)*765)#表示色の決定 if V<256: G=int(55+200*V/255); R=0 elif V<512: V-=256; R=int(V) ; G=255 else : V-=512; R=255 ; G=255-int(V) C="#%02X%02X%02X" % (chkC(R),chkC(G),B) canvas.itemconfig(GID[i][j], fill=C) #■矢印の描画 def drawArrow(canvas, X1,Y1,X2,Y2, color="black",width=1, tag="Arrow", r=10): canvas.create_line(X1,Y1,X2,Y2,fill=color,width=width,tag=tag) DX=X2-X1; DY=Y2-Y1; TH=math.pi/2-math.atan2(DY,DX) A=math.pi/8; TH2=TH-A;TH3=TH+A P=[[0,0], [0,0],[0,0]]; P[0][0]=X2; P[0][1]=Y2 P[1][0]=X2-r*math.sin(TH2); P[1][1]=Y2-r*math.cos(TH2) P[2][0]=X2-r*math.sin(TH3); P[2][1]=Y2-r*math.cos(TH3) canvas.create_line(P[0][0],P[0][1],P[1][0],P[1][1], fill=color,width=width,tag=tag) canvas.create_line(P[0][0],P[0][1],P[2][0],P[2][1], fill=color,width=width,tag=tag) def dspAbsVel(canvas,KX, KY):#■速度絶対値表示 global pbAbsV, pbU,pbV for i in range(1, KX): for j in range(1, KY): pbAbsV[i][j]=math.sqrt(pbU[i][j]*pbU[i][j]+pbV[i][j]*pbV[i][j]) dspCell(canvas, pbAbsV, KX, KY) def dspVel(canvas,KX, KY):#■速度を矢印で描画 canvas.delete("Arrow") for i in range(1, KX): X1=i*10+5 for j in range(1, KY): if flowExist[i][j]!=0: Y1=j*10+5;X2=X1+pbU[i][j]*10; Y2=Y1+pbV[i][j]*10 drawArrow(canvas, X1,Y1,X2,Y2, color="#0000FF", r=3) #■メイン処理 root=Tk(); root.title("障害物のある流体") canvas=Canvas(root, width=NumX*10+10,height=NumY*10+10, bg="white") canvas.pack() initDsp();dspObstacle();root.update() Ndsp=1 for L in range(cndLoop): boundVel() #　速度境界条件 tempVel() # 仮の速度を求める PoissonRight() # ポアソン方程式の右辺を計算する G1=compPressCond() # 圧力収束条件の計算 for kk in range(pbIter): #　収れん計算 G2=compPressCond() # 圧力収束条件の計算 if abs(1-G1/G2)<=0.00001:break G1=G2 #print(kk, G2,G1/G2) updateVel() # 仮の速度を圧力により変更 if L % Ndsp==0: root.title("障害物のある流体 計算回数 %d" % (L)) dspCheck=dspType & 3 if dspCheck==1: dspAbsVel(canvas,NumX, NumY) #速度絶対値表示 elif dspCheck==2: dspCell(canvas, pbP, NumX, NumY) #X方向圧力表示 elif dspCheck==3: dspCell(canvas, pbQ, NumX, NumY) #Y方向圧力表示 if dspType>=4 : dspVel(canvas,NumX, NumY) #速度(矢印)表示 root.update() time.sleep(0.0001) ５．５　円柱周りの粘性流体 【リスト5-8】円筒周りの粘性流体 from tkinter import * import random, math, time def array(N1, N2=0, N3=0): if N2 == 0 and N3==0: return [0 for i in range(N1)] elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)] else: return [[[0 for i in range(N3)] for j in range(N2)] for k in range(N1)] PI=math.pi MX=101; MY=101 VX =array(MX,MY); VY =array(MX,MY) P =array(MX,MY); Q =array(MX,MY); TMP=array(MX,MY) VXN=array(MX,MY); VYN=array(MX,MY) X =array(MX,MY); Y =array(MX,MY) NX0=40; NY0=40; NX=NX0+1; NY=NY0+1 RE = 5000; DT = 0.01; DY = 0.075 KK = 100; ACCPC=1 ; EPS=0.01; DX=2*PI/(NX-2) rDX = 1/DX; rDY=1/DY; rDT=1/DT; rRE=1/RE DX2 = rDX*rDX; DY2 = rDY*rDY rDX2= 1/(4*DX); rDY2 =1/(4*DY) FCT = 2/(2*DX2+2*DY2); for j in range(NY):#座標設定 for i in range(NX): THT=DX*(i-0.5); YY=DY*(j-0.5); THT=DX*(i-0.5) X[i][j]=math.exp(YY)*math.cos(THT) Y[i][j]=math.exp(YY)*math.sin(THT) for j in range(NY):#圧力と速度の初期設定 for i in range(NX): THT=DX*(i-1);P[i][j]=0.0; VX[i][j]=-math.cos(THT); VY[i][j]=math.sin(THT) def drawLine(cnvas, X1, Y1, X2, Y2, width=1, Color="black"):#■線分 return canvas.create_line(X1,Y1,X2,Y2,fill=Color, width=width) #多角形 def drawPolygon(canvas, poly, width=1, lColor="black", pColor="white"): return canvas.create_polygon(poly,fill=pColor, outline=lColor, width=width) def setX(X): return int(X*40+400) def setY(Y): return int(Y*40+250) def drawGrid(canvas): global NY0, NX0 for j in range(NY0): for i in range(NX0): drawLine(canvas, setX(X[i][j]), setY(Y[i][j]), setX(X[i+1][j]), setY(Y[i+1][j]),1,"#FF0000") drawLine(canvas, setX(X[i][j]), setY(Y[i][j]), setX(X[i][j+1]), setY(Y[i][j+1]),1,"#0000FF") for i in range(NX0): drawLine(canvas, setX(X[i][NY0]), setY(Y[i][NY0]), setX(X[i+1][NY0]), setY(Y[i+1][NY0]),1,"#FF0000") def average(A1, A2="", A3="", A4=""): if A2=="": return A1 if A3=="": return (A1+A2)*0.5 if A4=="": return (A1+A2+A3)/3 return (A1+A2+A3+A4)*0.25 def maxValue(P, nx, ny): MV=P[1][1] for i in range(1, nx-1): for j in range(1, ny-1): if P[i][j]>MV: MV=P[i][j] return MV def minValue(P, nx, ny): MV=P[1][1] for i in range(1, nx-1): for j in range(1, ny-1): if P[i][j] 255: return 255 if V < 0: return 0 return int(V) def color(V, minV, maxV):#■カラー設定 UDX=maxV-minV; B=0 if UDX<=0.0001: UDX=1 #最大最小の差がなければ一定値を基準とする V=int(((V-minV)/UDX)*765)#表示色の決定 if V<256: G=int(55+200*V/255); R=0 elif V<512: V-=256; R=V ; G=255 else : V-=512; R=255 ; G=255-V C="#%02X%02X%02X" % (chkC(R),chkC(G),B) return "#%02X%02X%02X" % (R, G, B) GID=array(NX,NY) LID=array(NX,NY) poly=array(8) def drawP(canvas): global P, NX, NY, GID,LID,poly MinV = minValue(P, NX, NY); MaxV=maxValue(P, NX,NY) for j in range(1,NY0): for i in range(1,NX0): poly[0]=setX(X[i][j]);poly[1]=setY(Y[i][j]) poly[2]=setX(X[i][j+1]);poly[3]=setY(Y[i][j+1]) poly[4]=setX(X[i+1][j+1]);poly[5]=setY(Y[i+1][j+1]) poly[6]=setX(X[i+1][j ]);poly[7]=setY(Y[i+1][j ]) C=color(P[i][j],MinV,MaxV) GID[i][j]=drawPolygon(canvas, poly,1, C, C) for j in range(1,NY0): for i in range(1,NX0): CX=average(X[i][j],X[i][j+1],X[i+1][j],X[i+1][j+1]) CY=average(Y[i][j],Y[i][j+1],Y[i+1][j],Y[i+1][j+1]) X1=setX(CX); Y1=setY(CY) X2=setX(CX+VXN[i][j]*1);Y2=setY(CY+VYN[i][j]*1) LID[i][j]=drawLine(canvas, X1,Y1, X2,Y2,1,"#0000CC") def updateP(canvas): global P, NX, NY, GID,LID,poly MinV = minValue(P, NX, NY); MaxV=maxValue(P, NX,NY) for j in range(1,NY0): for i in range(1,NX0): C=color(P[i][j],MinV,MaxV) canvas.itemconfig(GID[i][j],fill=C) for j in range(1,NY0): for i in range(1,NX0): CX=average(X[i][j],X[i][j+1],X[i+1][j],X[i+1][j+1]) CY=average(Y[i][j],Y[i][j+1],Y[i+1][j],Y[i+1][j+1]) X1=setX(CX); Y1=setY(CY) X2=setX(CX+VXN[i][j]*1);Y2=setY(CY+VYN[i][j]*1) canvas.coords(LID[i][j],X1,Y1, X2,Y2) def tempVelocity():#■速度設定 global NX,NY,NX0,NY0, VX,VY, VXN, VYN,TMP global rDX, rDY, DT for j in range(1,NY0):#　VY[]の計算 for i in range(1,NX): VYT=VY[i][j] VXR=average(VX[i-1][j],VX[i][j],VX[i-1][j+1],VX[i][j+1]) VYN[i][j]=(VYT*(VY[i+1][j]-VY[i-1][j])*rDX*0.5 + VXR*(VY[i][j+1]-VY[i][j-1])*rDY*0.5 + VY[i][j]*VXR - abs(VYT)*(VY[i+1][j]-2*VYT+VY[i-1][j])*rDX*0.5 - abs(VXR)*(VY[i][j+1]-2*VYT+VY[i][j-1])*rDY*0.5) for j in range(2,NY-2): for i in range(2,NX0): VYT=VY[i][j] VXR= average(VX[i-1][j],VX[i][j],VX[i-1][j+1],VX[i][j+1]) VYA=-VY[i+2][j]+8.0*(VY[i+1][j]-VY[i-1][j])+VY[i-2][j] VYB= VY[i+2][j]+VY[i-2][j]-4*(VY[i+1][j]+VY[i-1][j]) +6*VYT VYC=-VY[i][j+2]+8.0*(VY[i][j+1]-VY[i][j-1]) +VY[i][j-2] VYD= VY[i][j+2]+VY[i][j-2]-4*(VY[i][j+1]+VY[i][j-1]) +6*VYT VYN[i][j]=(VYT*VYA*rDX + abs(VYT)*VYB*rDX + VXR*VYC*rDY + abs(VXR)*VYD*rDY)/12.0 +VY[i][j]*VXR for j in range(1,NY0): for i in range(1,NX): YY=DY*(j-0.5); E1=math.exp(-YY); E2=math.exp(-2.0*YY) VYS=((VY[i+1][j]-2*VY[i][j]+VY[i-1][j])*DX2 + (VY[i][j+1]-2*VY[i][j]+VY[i][j-1])*DY2-VY[i][j] + (VX[i][j]+VX[i][j+1]-VX[i-1][j]-VX[i-1][j+1])*rDX) TMP[i][j]=VY[i][j]+DT*(-E1*VYN[i][j]+rRE*E2*VYS) for j in range(1,NY0): for i in range(1,NX): VY[i][j]=TMP[i][j] for j in range(1,NY):#VX[][]の計算 for i in range(1,NX0): VYT=average(VY[i][j-1],VY[i][j],VY[i+1][j-1],VY[i+1][j]) VXR=VX[i][j] VXN[i][j]=(VYT*(VX[i+1][j]-VX[i-1][j])*rDX*0.5 + VXR*(VX[i][j+1]-VX[i][j-1])*rDY*0.5 - VYT*VYT- abs(VYT)*(VX[i+1][j]-2*VXR+VX[i-1][j])*rDX*0.5 - abs(VXR)*(VX[i][j+1]-2*VXR+VX[i][j-1])*rDY*0.5) for j in range(2,NY0): for i in range(2,NX-2): VYT= average(VY[i][j-1],VY[i][j],VY[i+1][j-1],VY[i+1][j]) VXR= VX[i][j] VXA=-VX[i+2][j]+8*(VX[i+1][j]-VX[i-1][j])+VX[i-2][j] VXB= VX[i+2][j]+VX[i-2][j]-4*(VX[i+1][j]+VX[i-1][j]) +6*VXR VXC=-VX[i][j+2]+8*(VX[i][j+1]-VX[i][j-1])+VX[i][j-2] VXD= VX[i][j+2]+VX[i][j-2]-4*(VX[i][j+1]+VX[i][j-1]) +6*VXR VXN[i][j]=(VYT*VXA*rDX+abs(VYT)*VXB*rDX+ VXR*VXC*rDY+abs(VXR)*VXD*rDY)/12.0-VYT*VYT for j in range(1,NY): for i in range(1,NX0): YY=DY*(j-1); E1=math.exp(-YY); E2=math.exp(-2.0*YY) VXS=((VX[i+1][j]-2*VX[i][j]+VX[i-1][j])*DX2 + (VX[i][j+1]-2*VX[i][j]+VX[i][j-1])*DY2 -VX[i][j] - (VY[i+1][j]+VY[i+1][j-1]-VY[i][j]-VY[i][j-1])*rDX) TMP[i][j]=VX[i][j]+DT*(-E1*VXN[i][j]+rRE*E2*VXS) for j in range(1,NY0): for i in range(1,NX): VX[i][j]=TMP[i][j] def solvePoisson():#■ポアソン方程式の解 global NX,NY,NX0,NY0, VX,VY, VXN, VYN,TMP,Q global rDT,rDX,rDY,DY, DX2, DY2, KK, ACCPC for j in range(1,NY0): for i in range(1,NX0): VXR=average(VX[i][j+1],VX[i][j]) DVYDT=(VY[i+1][j]-VY[i][j])*rDX DVXDY=(VX[i][j+1]-VX[i][j])*rDY Q[i][j]=math.exp(DY*(j-0.5))*rDT*(DVXDY+VXR+DVYDT) FCT=1.0/(2*DX2+2*DY2) ERR0=0 for k in range(KK): ERR=0; j=0 for i in range(NX): P[i][j]=P[i][j+1]-2*rRE*math.exp(-(DY*j))*rDY*VX[i][j+2] P[i][NY0]=0 for j in range(NY): P[0][j]=P[NX-2][j]; P[NX0][j]=P[1][j]; for j in range(1,NY0): for i in range(1,NX0): RHS=((P[i+1][j]+P[i-1][j])*DX2+ (P[i][j+1]+P[i][j-1])*DY2-Q[i][j])*FCT dRHS = RHS-P[i][j] ERR += dRHS*dRHS P[i][j]=P[i][j]*(1-ACCPC)+RHS*ACCPC #print(k, ERR) if abs(1-ERR0/ERR)<0.000001: break ERR0=ERR return ERR; def newVelocity():#■圧力による速度変更 global NX,NY,NX0,NY0, VX,VY, VXN, VYN,TMP,P global rDT,rDX,rDY,DY, DX2, DY2, KK, ACCPC for j in range(1,NY0): for i in range(1,NX): YY=DY*(j - 0.5); E1=math.exp(-YY) VY[i][j]+= DT*(-E1*(P[i][j]-P[i-1][j])*rDX) for j in range(1,NY): for i in range(1,NX0): YY=DY*(j - 1); E1=math.exp(-YY) VX[i][j]+= DT*(-E1*(P[i][j]-P[i][j-1])*rDY) def boundaryCondition():#■境界条件 global NX,NY,NX0,NY0, VX,VY, VXN, VYN,TMP,P global rDT,rDX,rDY,DY, DX2, DY2, KK, ACCPC for i in range(NX): VY[i][0]=-VY[i][1];VX[i][0]=VX[i][2]; VX[i][1]=0 for i in range(NX): THT=DX*(i-1) VY[i][NY0]=math.sin(THT); VX[i][NY0]=-math.cos(THT); VX[i][NY]=-math.cos(THT) for j in range(NY): VX[0 ][j] = VX[NX-2][j]; VY[0 ][j] = VY[NX-2][j] VX[NX0][j] = VX[1 ][j]; VY[NX0][j] = VY[1 ][j] VX[NX ][j] = VX[2 ][j]; VY[NX ][j] = VY[2 ][j] def velocityVector():#■速度ベクトルの設定 global NX,NY,NX0,NY0, VX,VY, VXN, VYN,TMP,P global rDT,rDX,rDY,DY, DX2, DY2, KK, ACCPC for j in range(NY): for i in range(NX): THT=DX*(i-0.5) VXA=average(VX[i][j],VX[i+1][j]) VYA=average(VY[i][j],VY[i][j+1]) VXN[i][j]=VXA*math.cos(THT)-VYA*math.sin(THT) VYN[i][j]=VXA*math.sin(THT)+VYA*math.cos(THT) def execSim(N):#シミュレーション実行 for k in range(N): tempVelocity(); solvePoisson() newVelocity() ; boundaryCondition(); velocityVector(); #■メイン処理 root=Tk(); root.title("円柱回りの粘性流体") canvas=Canvas(root, width=500,height=500, bg="white");canvas.pack() drawGrid(canvas) drawP(canvas) for k in range(1000): root.title("円柱回りの粘性流体 計算回数 %d" % k) execSim(1); updateP(canvas); root.update() time.sleep(0.001)