QA@IT

分子動力学(速度ベルレ法)について

397 PV

下記のコードについてです。メインループは速度ベルレ法を使用してます。この際q_in.q_outの値が永遠に0となってしまいます。dqの値は出力されていますが次のメインループに行く前にq_in,q_outが、0クリアされてしまいます。どうかご享受お願いします。

コード

q_in=0.0

q_out=0.0



########x:周期境界#######

for i in range(N):

    if rx[i]>lx:

        rx[i]-=lx

    elif rx[i]<=0.0:

        rx[i]+=lx

    else:

        rx[i]=rx[i]



#########熱壁(y)###########



    if ry[i]<=0.5:                       

        vy_l=vy[i]

        vy[i]=canon_b(temp_l)

        dq=0.5*(vy[i]**2-vy_l**2)

if dq_1>=0.0:

q_in+=dq_1

else:

q_out+=dq_1

        q_in+=dq



        hit_l+=1









    elif ry[i]>=(ly-0.5):

        vy_l=vy[i]

        vy[i]=-canon_b(temp_h)

        dq=0.5*(vy[i]**2-vy_l**2)

if dq_2>=0.0:

q_in+=dq_2

else:

q_out+=dq_2

        q_in+=dq



        hit_h+=1

##########圧力、温度###########

temp_d=0.0

temp_u=0.0

N_d=0

N_u=0

for i in range(N):

if ry[i]<py:

temp_d+=(vx[i]2+vy[i]2)

N_d+=1

else:

temp_u+=(vx[i]2+vy[i]2)

N_u+=1

temp_d=temp_d/2

temp_u=temp_u/2

T_sum=0.0

for i in range(N):

    T=0.5*(vx[i]**2+vy[i]**2)

    T_sum+=T

T_res=T_sum*2/N
ウォッチ

この質問への回答やコメントをメールでお知らせします。