Relationship between critical depth and discharge coefficient for free overflow on dam spillway
In general, the discharge for free overflow on dam spillway can be expressed as follow:
On the other hand, a critical depth of rectangular cross section channel can be expressed as follow:
When it is assumed the water depth at the overflow crest of the dam is equal to the critical depth, the design head including approach velocity head can be expressed as follow:
Considering a unit length of overflow crest which means ,
From above, can be expressed as follow:
Therefore, discharge coefficient becomes a constant and it can be expressed as follow:
That is all.
Python:計算結果をエクセル出力し加工する(月平均値の算出とエクセル化)
能書き
報告書を作成しているとき、Pythonで計算した結果をWordに表形式で貼り付けたい場合がある。 これが、TeXやhtmlであれば、プロブラム内でタグを埋め込んだ形で標準出力し、コピー&ペーストでTeXもしくはhtml文書に貼り付けるのだが、相手がWordとなるとなかなか面倒である。PythonからWordを制御することができることも知っている(使ったことはない)が、そのためのみのコードを書くのも大変だし、後処理計算での活用などを考えるとエクセル出力しておき、それをwordに貼り付けるほうが、実用的である。
私がよく使う数表をwordに貼り付ける方法は以下の3つである。
(1) 計算結果をコンマ区切りで標準出力し、Wordにコピー&ペースト。その後はWordで[Insert]=>[table]=>[Convert Text to Table]=>[Separate test at Commas]=.[ok]の手順で、Word上で表に加工する方法。Word上で加工できるので仕上がりは結構綺麗にいく。
(2) 計算結果をエクセルに出力しフォーマットを整える。その後はエクセル上で範囲指定しクリップボードコピー、Word上で[Paste Special]=>[As Picture(PNG)]=>[ok]の手順で、画像としてWordに貼り付ける方法。画像としてwordに貼り付けるためWord上で加工できず、結果がぼやける場合もあり、最近は使わない。
(3)計算結果の表をhtml出力してブラウザで表示し、それをコピー&ペーストでwordに貼り付け、仕上げはword上で行う。ブラウザはsafariで行うのがやりやすいようである、
方法(3)は最近使い出したのだが、面倒なようでなかなか良い。 というのも計算結果は一回出力するのみでなく、次の処理で使う場合も多く、そうなると1ファイルに複数シートを配置できるエクセルで保存しておくのも後工程での使い勝手がよい。
エクセルからhtmlへの変換
エクセルからhtmlへの出力は、ネット上で便利なものが有り、以下のものを使わせてもらっている。
エクセルから変換したhtmlの表を貼り付けたものが下の表。二重線が消えているがhtml出力が最終成果ではないので気にしない。 必要あれば装飾をかければ良い。
Year | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ave |
2001 | 19.6 | 17.8 | 16.1 | 14.5 | 21.3 | 36.1 | 35.9 | 84.0 | 36.0 | 47.7 | 39.4 | 26.7 | 33.1 |
2002 | 22.1 | 20.4 | 18.5 | 16.7 | 24.1 | 19.9 | 43.1 | 36.7 | 44.8 | 32.9 | 29.6 | 23.8 | 27.8 |
2003 | 21.4 | 19.7 | 18.0 | 16.2 | 14.6 | 13.9 | 20.6 | 22.9 | 30.7 | 29.9 | 28.1 | 21.5 | 21.5 |
2004 | 19.2 | 17.6 | 15.9 | 14.3 | 21.1 | 37.7 | 47.0 | 48.7 | 110.6 | 47.1 | 32.0 | 24.2 | 36.2 |
2005 | 22.0 | 20.0 | 18.1 | 16.3 | 14.9 | 15.0 | 44.5 | 49.7 | 92.6 | 45.4 | 32.2 | 25.6 | 33.0 |
2006 | 22.3 | 20.5 | 18.6 | 16.7 | 15.2 | 16.1 | 32.7 | 65.6 | 104.4 | 62.3 | 34.7 | 25.6 | 36.3 |
2007 | 22.8 | 20.8 | 18.8 | 16.9 | 17.7 | 20.3 | 24.5 | 35.7 | 29.5 | 33.4 | 26.8 | 21.8 | 24.1 |
2008 | 20.0 | 18.3 | 16.6 | 15.0 | 13.7 | 25.0 | 56.6 | 83.9 | 74.4 | 41.0 | 36.4 | 28.9 | 35.9 |
2009 | 24.1 | 22.2 | 20.1 | 18.1 | 16.4 | 15.3 | 18.0 | 37.2 | 31.9 | 29.1 | 25.4 | 20.2 | 23.2 |
2010 | 18.6 | 17.0 | 15.3 | 13.8 | 12.4 | 11.8 | 34.5 | 114.9 | 89.1 | 45.6 | 32.9 | 24.4 | 36.0 |
2011 | 21.8 | 19.9 | 18.0 | 16.4 | 17.9 | 28.9 | 42.7 | 129.4 | 111.2 | 55.2 | 36.3 | 27.8 | 43.9 |
2012 | 24.9 | 22.7 | 20.4 | 18.4 | 17.8 | 25.9 | 41.6 | 60.0 | 72.2 | 42.4 | 33.8 | 26.5 | 33.9 |
2013 | 23.6 | 21.6 | 19.6 | 17.6 | 15.9 | 16.2 | 28.0 | 103.1 | 82.5 | 52.8 | 38.4 | 28.2 | 37.4 |
2014 | 23.7 | 21.8 | 19.8 | 17.8 | 16.0 | 18.3 | 57.8 | 56.8 | 58.2 | 34.4 | 28.8 | 23.1 | 31.5 |
2015 | 21.2 | 19.3 | 17.5 | 15.7 | 14.6 | 18.5 | 20.7 | 43.4 | 29.2 | 28.2 | 24.8 | 20.3 | 22.8 |
2016 | 18.8 | 17.1 | 15.5 | 14.0 | 12.6 | 13.6 | 39.1 | 48.7 | 58.3 | 35.8 | 29.9 | 23.9 | 27.3 |
2017 | 21.0 | 19.5 | 17.7 | 16.0 | 14.5 | 14.5 | 29.3 | 36.6 | 51.0 | 33.2 | 29.9 | 23.8 | 25.6 |
2018 | 21.1 | 19.6 | 17.8 | 16.1 | 14.8 | 33.8 | 40.6 | 58.9 | 41.9 | 64.7 | 35.2 | 26.6 | 32.7 |
2019 | 23.2 | 21.3 | 19.4 | 17.4 | 15.6 | ||||||||
2020 | |||||||||||||
Ave | 21.7 | 19.8 | 18.0 | 16.2 | 16.4 | 21.2 | 36.5 | 62.0 | 63.8 | 42.3 | 31.9 | 24.6 | 31.2 |
RGS1 (CA=3062km2) tank model monthly and annual discharge, unit: m3/s |
プログラム上の処理
Pythonプログラム上での出力したい表の処理は以下の手順で行う。
- 1.pandasのデータフレーム作成
- 2.一度データフレームをそのままエクセルに保存
- 3.保存したエクセルファイルを再度pythonで呼び出し、エクセうの書式指定、罫線などの可能を行い、再保存する。
データフレーム作成
df = pd.DataFrame({ 'Year' : lyy, 'Jan' : qqm[0:n,0], 'Feb' : qqm[0:n,1], 'Mar' : qqm[0:n,2], 'Apr' : qqm[0:n,3], 'May' : qqm[0:n,4], 'Jun' : qqm[0:n,5], 'Jul' : qqm[0:n,6], 'Aug' : qqm[0:n,7], 'Sep' : qqm[0:n,8], 'Oct' : qqm[0:n,9], 'Nov' : qqm[0:n,10], 'Dec' : qqm[0:n,11], 'Ave' : qqm[0:n,12] }) df = df.set_index('Year')
データフレーム保存
関連する内容を、1ファイルに、複数シートに分割して保存する。
fnameW='out_dis_month.xlsx' with pd.ExcelWriter(fnameW) as writer: df1_tank.to_excel(writer, sheet_name='RGS1_tank') df2_tank.to_excel(writer, sheet_name='RGS2_tank') df1_mlp.to_excel(writer, sheet_name='RGS1_mlp') df2_mlp.to_excel(writer, sheet_name='RGS2_mlp')
エクセルファイルの再読み込みと保存
この間に処理を記載する。
wb =openpyxl.load_workbook(fnameW) # .... # wb.save(fnameW)
フォーマット指定
# format for i in range(2, n+1): for j in range(2, m+1): ws.cell(row=i, column=j).number_format = '0.0'
罫線(前半は1重線、後半は2重線を指定)
from openpyxl.styles.borders import Border, Side # border line bcc = Border(top=Side(style='thin', color='000000'), bottom=Side(style='thin', color='000000'), left=Side(style='thin', color='000000'), right=Side(style='thin', color='000000') ) for i in range(1, n+1): for j in range(1, m+1): ws.cell(row=i, column=j).border = bcc bcc = Border(top=Side(style='double', color='000000'), bottom=Side(style='double', color='000000'), left=Side(style='double', color='000000'), right=Side(style='double', color='000000') ) for j in range(1, m+1): ws.cell(row=n, column=j).border = bcc for i in range(1, n+1): ws.cell(row=i, column=m).border = bcc
セルの結合(長い文字列を入力しセルを結合する)
ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) tank model monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)
セルの塗りつぶし
このプログラムには含まれていないが、塗りつぶしを行う場合は以下の要領で行う。
from openpyxl.styles import PatternFill # cell color fill = PatternFill(patternType='solid', fgColor='00ffff') for i in [6,11,16,21,26,31]: for j in range(3,17): ws.cell(row=i,column=j).fill = fill
主要計算部分の説明
例外処理を使っているので、その部分を解説。 このプログラムでは、1日でも欠測があれば、その月および年の平均は欠測としている。
main
で数表にする年と月を、リスト lyy
と lmm
で定義している。
lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009','2010', '2011','2012','2013','2014','2015','2016','2017','2018','2019','2020'] lmm=['01','02','03','04','05','06','07','08','09','10','11','12']
以下の関数で各月および各年の平均値を算出している。 引数の意味は以下の通り。
rf
: 日データを示すseries.lyyy
: リストlyy
と同じlmmm
: リストlmm
と同じ。
def cal_m(rf,lyyy,lmmm): qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) nac=0 for i,yy in enumerate(lyyy): for j,mm in enumerate(lmmm): ss=yy+'/'+mm # 処理する年/月を指定 try: # Series: rf[ss]に対する処理(合計およびデータ数カウント)を行う qqm[i,j]=rf[ss].sum() # sum in a month kda[i,j]=rf[ss].count() # availavle days nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days except KeyError: # もし[ss]に合致する年/月のデータがなければ配列にnp.nanを代入 qqm[i,j]=np.nan kda[i,j]=np.nan if 0<nac: # 仮に指定した年/月のデータがあってもnp.nanがあれば、すなわち月の全日のデータが揃っていなければnp.nanを代入 qqm[i,j]=np.nan kda[i,j]=np.nan qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year 年間合計値 kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days 年間合計値 for j in range(len(lmmm)+1): qqm[-1,j]=np.nansum(qqm[:,j]) # nanを除外したデータの合計(各月合計値) kda[-1,j]=np.nansum(kda[:,j]) # nanを除外したデータ数の合計(各月合計値) qmean=qqm/kda # 月および年のデータ合計をデータ数で除すことにより平均値を算出 return qmean
プログラム例(全文)
プログラムは、ファイルから値を読み取り、上に示した表を作るもの。 実際のプログラムの出力はエクセル。 以下にPythonによるプログラム全文を示す。
import pandas as pd import numpy as np import datetime import openpyxl from openpyxl.styles.borders import Border, Side def formx(ws,n,m): # format for i in range(2, n+1): for j in range(2, m+1): ws.cell(row=i, column=j).number_format = '0.0' # border line bcc = Border(top=Side(style='thin', color='000000'), bottom=Side(style='thin', color='000000'), left=Side(style='thin', color='000000'), right=Side(style='thin', color='000000') ) for i in range(1, n+1): for j in range(1, m+1): ws.cell(row=i, column=j).border = bcc bcc = Border(top=Side(style='double', color='000000'), bottom=Side(style='double', color='000000'), left=Side(style='double', color='000000'), right=Side(style='double', color='000000') ) for j in range(1, m+1): ws.cell(row=n, column=j).border = bcc for i in range(1, n+1): ws.cell(row=i, column=m).border = bcc def wexcel(df1_tank,df2_tank,df1_mlp,df2_mlp): fnameW='out_dis_month.xlsx' with pd.ExcelWriter(fnameW) as writer: df1_tank.to_excel(writer, sheet_name='RGS1_tank') df2_tank.to_excel(writer, sheet_name='RGS2_tank') df1_mlp.to_excel(writer, sheet_name='RGS1_mlp') df2_mlp.to_excel(writer, sheet_name='RGS2_mlp') n=len(df1_tank.index)+1 m=len(df1_tank.columns)+1 wb =openpyxl.load_workbook(fnameW) # ws = wb.get_sheet_by_name('RGS1_tank') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS1 (CA=3062km2) tank model monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # ws = wb.get_sheet_by_name('RGS2_tank') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) tank model monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # ws = wb.get_sheet_by_name('RGS1_mlp') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS1 (CA=3062km2) MLP monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # ws = wb.get_sheet_by_name('RGS2_mlp') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) MLP monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # wb.save(fnameW) def rdata(fnameR): df=pd.read_csv(fnameR, header=0, index_col=0) # read excel data df.index = pd.to_datetime(df.index, format='%Y/%m/%d') return df def cal_m(rf,lyyy,lmmm): qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) nac=0 for i,yy in enumerate(lyyy): for j,mm in enumerate(lmmm): ss=yy+'/'+mm try: qqm[i,j]=rf[ss].sum() # sum in a month kda[i,j]=rf[ss].count() # availavle days nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days except KeyError: qqm[i,j]=np.nan kda[i,j]=np.nan if 0<nac: qqm[i,j]=np.nan kda[i,j]=np.nan qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days for j in range(len(lmmm)+1): qqm[-1,j]=np.nansum(qqm[:,j]) kda[-1,j]=np.nansum(kda[:,j]) qmean=qqm/kda return qmean def qq_mon(lyy,lmm,qq): qqm=cal_m(qq,lyy,lmm) lyy=lyy+['Ave'] n=len(lyy) df = pd.DataFrame({ 'Year' : lyy, 'Jan' : qqm[0:n,0], 'Feb' : qqm[0:n,1], 'Mar' : qqm[0:n,2], 'Apr' : qqm[0:n,3], 'May' : qqm[0:n,4], 'Jun' : qqm[0:n,5], 'Jul' : qqm[0:n,6], 'Aug' : qqm[0:n,7], 'Sep' : qqm[0:n,8], 'Oct' : qqm[0:n,9], 'Nov' : qqm[0:n,10], 'Dec' : qqm[0:n,11], 'Ave' : qqm[0:n,12] }) df = df.set_index('Year') return df def main(): lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009','2010', '2011','2012','2013','2014','2015','2016','2017','2018','2019','2020'] lmm=['01','02','03','04','05','06','07','08','09','10','11','12'] fname_list=[ 'df_rgs1_tank_result.csv', 'df_rgs2_tank_result.csv', 'df_rgs1_mlp_result.csv', 'df_rgs2_mlp_result.csv' ] for iii,fnameR in enumerate(fname_list): df0=rdata(fnameR) # daily discharge df0=df0['2001/01/01':'2019/05/31'] qq=pd.Series(df0['Q_pred'], index=df0.index) df=qq_mon(lyy,lmm,qq) print(fnameR) if iii==0: df1_tank=df if iii==1: df2_tank=df if iii==2: df1_mlp=df if iii==3: df2_mlp=df wexcel(df1_tank,df2_tank,df1_mlp,df2_mlp) #============== # Execution #============== if __name__ == '__main__': main()
以 上
GraphvizでMLPと通常の重回帰分析のイメージ図を作る
作例
GraphvizでMLPと通常重回帰分析のイメージ図(conceptual diagram あるいは schematic diagram かな)を作ってみた。 作例は以下の通り。
プログラミング環境は以下の通り。
Graphvizのインストール
まずはGraphvizのインストール。 brewとpipで2回行う。
brew install graphviz pip install graphviz
MLPのイメージ図作成プログラム
node の表示順を制御するため、gi.node('1','x[0]' ,pos='0,0!')
というように pos
を追記している。
from graphviz import Digraph g = Digraph(format='png') g.attr(rankdir='LR') g.attr(splines='false') g.attr(dpi='300') # Input subgraph gi=Digraph(name='cluster_i') gi.attr(label='Input') gi.attr(penwidth='0') gi.node('1','x[0]' ,pos='0,0!') gi.node('2','x[1]' ,pos='0,1!') gi.node('3','x[..]' ,pos='0,2!') gi.node('4','x[n]',pos='0,3!') # Hidden subgraph 1 gh1=Digraph(name='cluster_h1') gh1.attr(label='Hidden layer 1') gh1.attr(penwidth='0') gh1.node('5', 'h1[0]' ,pos='0,0!') gh1.node('6', 'h1[1]' ,pos='0,1!') gh1.node('7', 'h1[..]' ,pos='0,2!') gh1.node('8', 'h1[n]',pos='0,3!') # Hidden subgraph 2 gh2=Digraph(name='cluster_h2') gh2.attr(label='Hidden layer 2') gh2.attr(penwidth='0') gh2.node('9', 'h2[0]' ,pos='0,0!') gh2.node('10','h2[1]' ,pos='0,1!') gh2.node('11','h2[..]' ,pos='0,2!') gh2.node('12','h2[n]',pos='0,3!') # Output subgraph go=Digraph(name='cluster_o') go.attr(label='Output') go.attr(penwidth='0') go.node('13','y',pos='0,1.5!') g.subgraph(gi) g.subgraph(gh1) g.subgraph(gh2) g.subgraph(go) g.edge('1', '5') g.edge('1', '6') g.edge('1', '7') g.edge('1', '8') g.edge('2', '5') g.edge('2', '6') g.edge('2', '7') g.edge('2', '8') g.edge('3', '5') g.edge('3', '6') g.edge('3', '7') g.edge('3', '8') g.edge('4', '5') g.edge('4', '6') g.edge('4', '7') g.edge('4', '8') g.edge('5', '9') g.edge('5', '10') g.edge('5', '11') g.edge('5', '12') g.edge('6', '9') g.edge('6', '10') g.edge('6', '11') g.edge('6', '12') g.edge('7', '9') g.edge('7', '10') g.edge('7', '11') g.edge('7', '12') g.edge('8', '9') g.edge('8', '10') g.edge('8', '11') g.edge('8', '12') g.edge('9', '13') g.edge('10', '13') g.edge('11', '13') g.edge('12', '13') g.render('fig_3_mlp', view=True)
重回帰のイメージ図作成プログラム
edge(矢印線)に回帰係数を示す w[..]
を表示するため、g.edge('3', '5', label='w[..]')
というように label
を追記している。
from graphviz import Digraph g = Digraph(format='png') g.attr(rankdir='LR') g.attr(splines='false') g.attr(dpi='300') # Input subgraph gi=Digraph(name='cluster_i') gi.attr(label='Input') gi.attr(penwidth='0') gi.node('1','x[0]' ,pos='0,0!') gi.node('2','x[1]' ,pos='0,1!') gi.node('3','x[..]' ,pos='0,2!') gi.node('4','x[n]',pos='0,3!') # Output subgraph go=Digraph(name='cluster_o') go.attr(label='Output') go.attr(penwidth='0') go.node('5','y',pos='0,0!') g.subgraph(gi) g.subgraph(go) g.edge('1', '5', label='w[0]') g.edge('2', '5', label='w[1]') g.edge('3', '5', label='w[..]') g.edge('4', '5', label='w[n]') g.render('fig_3_reg', view=True)
以 上
Structural design formulas for penstocks embedded in rock
The structural design formulas for penstocks embedded in rock which are shown in 'Technical standards for gates and penstock' in Japan are described.
Allowable stresss of steel material
Welded joint efficiency
Location of welding | RT or UT carrying out for more than 5% of welded joint length |
RT or UT carrying out for less than 5% of welded joint length |
---|---|---|
Factory welding | 95% | 85% |
Site welding | 90% | 80% |
Required safety factor for buckling due to external pressure
Pipe shell
Stiffener
Design internal pressure
Definition of dimensions
Corrosion allowance
An allowance of more than or equal to 1.5mm shall be provided against corrossion and wear.
Minimum shell thickness
The minimum shell thickness used for the pressure lining part shall be more than or equal to those determined from following formula, if stiffeners are not used.
The minimum shell thickness shall not be less than 6mm even if the pipe diameter is small and stiffeners are used.
Stress calculation
Tensile stress in circumferential direction due to internal pressure
General
Internal pressure shared design by bedrock
Tensile stress in penstock axis direction
Stress due to restraint by stiffener
Stress due to temperature change
Stress due to Poisson's effect
Calculation for external presure
Without stiffener (E. Amstutz's formula)
With stiffener
Pipe shell proper (S. Timoshenko's formula)
Stiffener (E. Amstutz's formula)
That's all.
Python 岩盤内埋設式水圧鉄管設計プログラム(新版)
はじめに
岩盤内埋設式水圧鉄管プログラムは、過去に一度掲載している(以下参照)。今回のプログラミング内容は工学的には過去のものとダブルが、過去投稿も、エクセルの罫線の書き方やセルの色付などの参考になるため、削除せずにおく。
今回は、以前のプログラムに対し、更に汎用性をもたせるよう若干の変更を加えるとともに、エクセル書出部があまりにも冗長なので、プログラムを「設計計算部」と「エクセル書出部」に分離した。
「設計計算部」ではデータを読み込み、必要諸元を定めたらデータフレームを作成し、それをエクセル形式で出力するまでとした。具体的には、1つの出力用temporaryファイル(Excel)に、3つのシートを作成し、入力データ、耐内圧設計結果、耐外圧設計結果を格納したデータフレームを保存している。
with pd.ExcelWriter('_df_ps.xlsx') as writer: # temporary file of outpit df0.to_excel(writer, sheet_name='Load', index=False) # dataframe for input data df1.to_excel(writer, sheet_name='Pin', index=False) # dataframe for design against internal pressure df2.to_excel(writer, sheet_name='Pex', index=False) # dataframe for design against external pressure
この段階では出力数値のフォーマットが統一されておらず、きれいとは言えない。
このため、「エクセル書出部」で、「設計計算部」で吐き出したエクセル形式のデータフレームを読み込み、数値の書式を揃えてエクセルに書き出すようにした。この段階で、等幅フォントの使用、セル幅の調整、各セルを破線で囲む、というところまで行っている。これ以上の装飾はプログラムが煩雑になるとともに汎用性をもたせるのが難しいため、以降はエクセルから装飾するという考え方にしている。
設計計算式は、日本の水門鉄管技術基準に基づいている。なお、本プログラムでは、管軸方向応力の検討は行っていない。
出力データ事例(file: out_penstock.xlsx)
入力データ出力(sheet: Load)
耐内圧設計結果(sheet: Pin)
耐外圧設計結果(sheet: Pout)
プログラミング
入出力データファイル
プログラム内で指定している入出力ファイルは以下の通り。 設計計算部の入力データファイル「inp_load.xlsx」の内容は、エクセル書出部出力の入力データ出力と同一(ただし記号の凡例は含まず)である。
プログラム | 入力データファイル | 出力データファイル |
---|---|---|
設計計算部 | inp_load.xlsx | _df_ps.xlsx |
エクセル書出部 | _df_ps.xlsx | out_penstock.xlsx |
エクセル書き出し部の出力は3つのシートからなっており、それぞれ以下の内容が出力されている。
Sheet | 内容 |
---|---|
Load | 入力データ(板割り・荷重計算・岩盤物性) |
Pin | 耐内圧設計結果 |
Pex | 耐外圧設計結果 |
入出力項目
以下で示す入力データおよび出力データは、考えるべき「ある長さを有する鉄管」の「下流端」を示すことに注意する。
列名 | 説明 |
---|---|
No | 通し番号 |
L | 区間鉄管長さ |
Sum(L) | 鉄管累計長さ(サージタンクからの長さ) |
EL | 鉄管中心標高 |
D0 | 鉄管設計内径 |
GWL | 鉄管中心に対する地表標高 |
Hst | 静水頭 |
Hsg | サージング水頭(一定値) |
Hwh | 水撃圧水頭(水車中心で最大、サージタンクでゼロ) |
Hin | 設計内水頭(Hst + Hsg + Hwh) |
Hex | 設計外水頭(GWL - EL) |
Dr | トンネル掘削径(岩盤負担設計用) |
Eg | 岩盤弾性係数(岩盤負担設計用) |
Remarks | 部位説明 |
列名 | 説明 |
---|---|
No | 通し番号 |
Pi | 設計内水圧(設計内水頭 x 0.01) |
L | 区間鉄管長さ |
D0 | 鉄管設計内径 |
t0 | 鉄管設計板厚 |
steel | 鋼種 |
Eg | 岩盤弾性係数 |
lam | 岩盤負担率 |
sig | 鉄管発生応力 |
siga | 鉄管許容応力 |
weight | 区間鉄管管胴重量 |
Remarks | 部位説明 |
列名 | 説明 |
---|---|
No | 通し番号 |
Pe | 設計外水圧(設計外水頭 x 0.01) |
L | 区間鉄管長さ |
D0 | 鉄管設計内径 |
t0 | 鉄管設計板厚 |
steel | 鋼種 |
pk_0 | 限界座屈圧力(補剛材無し) |
SF_0 | 外水圧に対する安全率(補剛材無し) |
pk_s | 限界座屈圧力(補剛材有り) |
SF_s | 外水圧に対する安全率(補剛材有り) |
sc_r | 補剛材の限界座屈応力 |
sc_c | 補剛材の合成圧縮応力 |
SF_c | 補剛材の合成圧縮応力に対する安全率 |
stiffener | 補剛材ピッチ |
プログラム内で指定している変数
設計計算に必要な定数は、入力データの他、以下のものをプログラム内で定義している。もちろん値の変更は可能である。
変数 | 説明 |
---|---|
E_steel = 206000 | 鉄管弾性係数(MPa) |
po_steel = 0.3 | 鉄管ポアソン比 |
t_eps = 2.0 | 鉄管余裕厚(mm) |
ef = 0.85 | 溶接継手効率 |
t0_rec = 25 | 推奨最小板厚 |
(note) 推奨最小板厚 t0_rec
は、水門鉄管技術基準で定める最小板厚とは別に、施工性・耐外圧設計のために独自に定める最小板厚である。計算される板厚はこの値以上となる。計算された補剛材ピッチが小さすぎる場合など耐外圧設計において管胴本体の耐外圧性能を高めたい場合、溶接施工用台車吊り下げのための必要最小板厚を指定したい場合に用いる。
変数 | 説明 |
---|---|
alpha = 1.2e-5 | 鉄管の熱膨張係数(1/度) |
deltaT = 20.0 | 鉄管の温度降下量(度) |
Ec = 20600.0 | 充填コンクリートの弾性係数(MPa) |
Bc = 0.0 | 充填コンクリートの塑性変形係数 |
Bg = 1.0 | 岩盤の塑性変形係数 |
mg = 4.0 | 岩盤のポアソン数(岩盤のポアソン比 = 0.25) |
(note) 鉄管と重点コンクリート間の空隙は、 で計算している。
鋼種 | 適用板厚 | 許容応力 | 降伏点 |
---|---|---|---|
JIS: SM400 | t0 < 40mm | 130 MPa | 235 MPa |
JIS: SM490 | t0 < 40mm | 175 MPa | 315 MPa |
JIS: SM570 | t0 < 40mm | 240 MPa | 450 MPa |
JIS: SM570 | 40mm < t0 | 235 MPa | 430 MPa |
(note) 設計計算では耐内圧設計において上記鋼種を自動的に選定する。SM400、SM490では許容応力が40mmで切り替わることを考慮し板厚は40mm以下としている。このプログラム内では使用鋼種はJIS: SM570までとしているが、JIS: SHY685を加えることにより、更に高落差への対応が可能である。
変数 | 説明 |
---|---|
hr=75 | 補剛材高さ(mm) |
tr=20 | 補剛材板厚(mm) |
el=3000, 1500, 1000 | 補剛材ピッチ(mm) |
(note) 補剛材設計では、補剛材寸法を、hr x tr = 75 x 20 と設定し、補剛材ピッチを el=3000, 1500, 1000と変化させて、外水圧に対する安全率が、管胴・補剛材とも1.5以上となるピッチを選定している。所要安全率が満たされない場合、補剛材寸法・ピッチ共に変更可能であるが、Global変数として定義している推奨最小板厚(変数名:t0_rec)を増加させて管胴の外水圧に対する安全率を高める方法もある。なお、補剛材ピッチは3m単位管に設置することを意識しており、また補剛材高さはトンネル掘削径との兼ね合いを考慮して設定することに注意する。
プログラムコード
設計計算プログラム
# ================================ # Embedded penstock Design # ================================ import numpy as np import pandas as pd from scipy import optimize import openpyxl # declaration of global variables E_steel = 206000 # elastic mpdulus of steel (MPa) po_steel = 0.3 # Poisson's ratio of steel t_eps = 2.0 # margin thickness (mm) ef = 0.85 # welding joint efficiency t0_rec = 25 # recommended minimum plate thickness def stcr(dd0, t0, hr, tr, el, pp, sigF): # Calculation of critical stress and compressive stress of stiffener # dd0 : design internal diameter of pipe # t0 : design plate thickness # hr : height of stiffener # tr : plate thickness of stiffener # el : stiffener interval # pp : external pressure # sigF : yield point of steel plate material def calrg(dd0, t0, eps, hr, tr): # calculation of section properties rm = 0.5*(dd0+t0) tt = t0-eps # change symbols b = 1.56*np.sqrt(rm*tt)+tr d = tt+hr s = tt t = tr # calculation of radius of gyration of area e2 = (d**2*t+s**2*(b-t))/2/(b*s+t*(d-s)) e1 = d-e2 ii = 1/3*(t*e1**3+b*e2**3-(b-t)*(e2-s)**3) aa = s*b+t*(d-s) rg = np.sqrt(ii/aa) # radius of gyration of area ee = np.abs(e2-s) return rg, ee def func(x, params): # function for Brent-method for obtaining sigN # x=sigN k0 = params[0] rm = params[1] t = params[2] Es = params[3] sigF = params[4] rg = params[5] ee = params[6] aa = k0/rm+x/Es bb = 1+rm**2/rg**2*x/Es cc = 1.68*rm/ee*(sigF-x)/Es*(1-1/4*rm/ee*(sigF-x)/Es) f = aa*bb**1.5-cc return f # main routine Es = E_steel # elastic modulus of steel plate pos = po_steel # Poisson's ratio of steel plate eps = t_eps # margin of plate thickness t = t0-eps rm = (dd0+t0)/2 r0d = (dd0+2*t0)/2 k0 = 0.4e-3*rm # calculation of critical stress rg, ee = calrg(dd0, t0, eps, hr, tr) params = np.array([k0, rm, t, Es, sigF, rg, ee]) sigN = optimize.brentq(func, 0.0, 5*sigF, args=(params)) scr = sigN*(1-r0d/ee*(sigF-sigN)/(1+3/2*np.pi)/Es) # calculation of compressive stress s0 = tr*(t+hr) iis = tr/12*(hr+t)**3 beta = (3*(1-pos**2))**0.25/np.sqrt(rm*t) c1 = r0d**2/t-(tr+1.56*np.sqrt(rm*t))*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) c2 = 3/(3*(1-pos**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sin(beta*el)) / \ (np.cosh(beta*el)-np.cos(beta*el))+2*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) pd = pp/(tr+1.56*np.sqrt(rm*t))*((tr+1.56*np.sqrt(rm*t))+2*c1/c2) sc = pd*r0d*(tr+1.56*np.sqrt(rm*t))/(s0+1.56*t*np.sqrt(rm*t)) return scr, sc def timo(el, hr, tr, dd0, t0): # Calculation of criticsl buckling pressure # of steel pipe with stiffener by Timoshenko's formula # el : interval of stiffener # hr : height of stiffener # tr : plate thickness of stiffener # dd0 : design internal diameter of steel pipe # t0 : design plate thickness of steel pipe pi = np.pi Es = E_steel # elastic modulus of steel pos = po_steel # Poisson's ratio of steel eps = t_eps # margin of plate thickness rm = 0.5*(dd0+t0) r0d = 0.5*(dd0+2*t0) t = t0-eps s0 = tr*(t+hr) iis = tr/12*(hr+t)**3 beta = (3*(1-pos**2))**0.25/np.sqrt(rm*t) c1 = r0d**2/t-(tr+1.56*np.sqrt(rm*t))*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) c2 = 3/(3*(1-pos**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sin(beta*el)) / \ (np.cosh(beta*el)-np.cos(beta*el))+2*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) lt = 2/(tr+1.56*np.sqrt(rm*t))*c1/c2 lam = 1-(1+lt)*(1+tr/(1.56*np.sqrt(rm*t)))/(1+s0/(1.56*t*np.sqrt(rm*t))) ell = (el+1.56*np.sqrt(rm*t)*np.arccos(lam))*(1+0.037 * np.sqrt(rm*t)/(el+1.56*np.sqrt(rm*t)*np.arccos(lam))*t**3/iis) apk = [] for n in range(2, 30): c1 = (1-pos**2)/(n**2-1)/(1+n**2*ell**2/pi**2/r0d**2)**2 c2 = t**2/12/r0d**2*((n**2-1)+(2*n**2-1-pos) / (1+n**2*ell**2/pi**2/r0d**2)) _pk = (c1+c2)*Es*t/(1-pos**2)/r0d apk = apk+[_pk] pk = min(apk) return pk def ams(dd0, t0, sigF): # Calculate criticsl buckling pressure # of steel pipe without stiffener by Amstutz's formula # dd0 : design internal diameter of steel pipe # t0 : design plate thickness of steel pipe # sigF : yield point of steel material def func(x, params): # function for Brent-method for obtaining sigN # x=sigN k0 = params[0] rm = params[1] t = params[2] Ess = params[3] sigFs = params[4] aa = k0/rm+x/Ess bb = 1+12*rm**2/t**2*x/Ess cc = 3.36*rm/t*(sigFs-x)/Ess*(1-1/2*rm/t*(sigFs-x)/Ess) f = aa*bb**1.5-cc return f def cal_pk(sigN, sigFs, Ess, rm, t): # calculation of critical buckling pressure aa = rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess) pk = sigN/aa return pk # main routine eps = t_eps # margin of plate thickness Es = E_steel # elastic modulus of steel plate pos = po_steel # Poisson's ratio of steel plate Ess = Es/(1-pos**2) mu = 1.5-0.5*1/(1+0.002*Es/sigF)**2 sigFs = mu*sigF/(1-pos+pos**2)**0.5 t = t0-eps rm = (dd0+t0)/2 k0 = 0.4e-3*rm params = np.array([k0, rm, t, Ess, sigFs]) sigN = optimize.brentq(func, 0.0, sigFs, args=(params)) pk = cal_pk(sigN, sigFs, Ess, rm, t) return pk def calc1(pp, dd0, sa, Eg, dr): # Design for internal pressure # (Calculate required plate thickness) # pp : internal pressure # dd0 : design internal diameter # sa : allowable stress of steel pipe # Eg : elastic modulus of bedrock # dr : excavation diameter of tunnel Es = E_steel # elastic modulus of steel eps = t_eps # margin thickness alpha = 1.2e-5 # thermal expansion coefficient of steel deltaT = 20.0 # temperatire change of steel Ec = 20600.0 # elastic modulus of backfill concrete Bc = 0.0 # plastic deformation modulus of cbackfill concrete Bg = 1.0 # plastoc deformation modulus of bedrock mg = 4.0 # Poisson number of bedrock t0t = t0_rec # recommended minimum thickness t0min = np.ceil((dd0+800)/400) dd = dd0-eps lam = 0 tt = pp*dd/2/sa t0 = np.ceil(tt+eps) if t0 < t0min: t0 = t0min if t0 < t0t: t0 = t0t # Internal pressure shared design by bedrock if 1 <= Eg: c1 = sa-Es*alpha*deltaT c2 = (1+Bc)*Es/Ec*2/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2/dd tt = pp*dd/2/sa-c1/c2/sa lam1 = 1-Es/pp*alpha*deltaT*2*tt/dd lam2 = 1+(1+Bc)*Es/Ec*2*tt/dd*np.log(dr/dd) + \ (1+Bg)*Es/Eg*(mg+1)/mg*2*tt/dd lam = lam1/lam2 t0 = np.ceil(tt+eps) if t0 < t0min: t0 = t0min if t0 < t0t: t0 = t0t tt = t0-eps lam1 = 1-Es/pp*alpha*deltaT*2*tt/dd lam2 = 1+(1+Bc)*Es/Ec*2*tt/dd*np.log(dr/dd) + \ (1+Bg)*Es/Eg*(mg+1)/mg*2*tt/dd lam = lam1/lam2 sig = pp*dd/2/(t0-eps)*(1-lam) return t0, sig, lam def datainp(): df0 = pd.read_excel('inp_load.xlsx') return df0 def main(): df0 = datainp() no = np.array(df0['No'], dtype=np.int) # pipe number al = np.array(df0['L(m)'], dtype=np.float64) # pipe length (m) d0 = np.array(df0['D0(m)'], dtype=np.float64) * \ 1000 # internal diameter (mm) pi = np.array(df0['Hin(m)'], dtype=np.float64) * \ 0.01 # internal pressure (MPa) pe = np.array(df0['Hex(m)'], dtype=np.float64) * \ 0.01 # external pressure (MPa) dte = np.array(df0['Dr(m)'], dtype=np.float64) * \ 1000 # tunnel excavation diameter (mm) # elastic modulus of bedrock (MPa) egg = np.array(df0['Eg(MPa)'], dtype=np.float64) a_t0 = np.zeros(len(no), dtype=np.float64) a_sig = np.zeros(len(no), dtype=np.float64) a_lam = np.zeros(len(no), dtype=np.float64) a_sa = np.zeros(len(no), dtype=np.float64) a_sname = [] a_pk0 = np.zeros(len(no), dtype=np.float64) a_sf0 = np.zeros(len(no), dtype=np.float64) a_pks = np.zeros(len(no), dtype=np.float64) a_sfs = np.zeros(len(no), dtype=np.float64) a_scr = np.zeros(len(no), dtype=np.float64) a_scc = np.zeros(len(no), dtype=np.float64) a_sfc = np.zeros(len(no), dtype=np.float64) a_stiff = [] for num, dd0, ppi, ppe, dr, Eg in zip(no, d0, pi, pe, dte, egg): sname = 'SM400' sa = 130*ef t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr) if sname == 'SM400' and 40 < t0: sname = 'SM490' sa = 175*ef t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr) if sname == 'SM490' and 40 < t0: sname = 'SM570' sa = 240*ef t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr) if sname == 'SM570' and 40 < t0: sname = 'SM570' sa = 235*ef t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr) a_t0[num-1] = t0 a_sig[num-1] = sig/ef a_lam[num-1] = lam a_sa[num-1] = sa/ef a_sname = a_sname+[sname] if sname == 'SM400': sigF = 235 if sname == 'SM490': sigF = 315 if sname == 'SM570' and t0 <= 40: sigF = 450 if sname == 'SM400' and 40 < t0: sigF = 430 # without stiffener pk0 = ams(dd0, t0, sigF) sf0 = pk0/ppe if sf0 < 1.5: # stiffener size hr = 75 # height of stiffener plate tr = 20 # thickness of stiffener plate # stiffener interval=2000mm el = 3000.0 pk1 = timo(el, hr, tr, dd0, t0) sf1 = pk1/ppe scr1, sc1 = stcr(dd0, t0, hr, tr, el, ppe, sigF) sfc1 = scr1/sc1 # stiffener interval=1500mm el = 1500.0 pk2 = timo(el, hr, tr, dd0, t0) sf2 = pk2/ppe scr2, sc2 = stcr(dd0, t0, hr, tr, el, ppe, sigF) sfc2 = scr2/sc2 # stiffener interval=1000mm el = 1000.0 pk3 = timo(el, hr, tr, dd0, t0) sf3 = pk3/ppe scr3, sc3 = stcr(dd0, t0, hr, tr, el, ppe, sigF) sfc3 = scr3/sc3 # store results a_pk0[num-1] = pk0 a_sf0[num-1] = sf0 if 1.5 <= sf0: a_pks[num-1] = 0 a_sfs[num-1] = 0 a_stiff = a_stiff+['no-need'] a_scr[num-1] = 0 a_scc[num-1] = 0 a_sfc[num-1] = 0 elif sf0 < 1.5 and 1.5 <= sf1: a_pks[num-1] = pk1 a_sfs[num-1] = sf1 a_stiff = a_stiff+['interval=3000mm'] a_scr[num-1] = scr1 a_scc[num-1] = sc1 a_sfc[num-1] = sfc1 elif sf1 < 1.5 and 1.5 <= sf2: a_pks[num-1] = pk2 a_sfs[num-1] = sf2 a_stiff = a_stiff+['interval=1500mm'] a_scr[num-1] = scr2 a_scc[num-1] = sc2 a_sfc[num-1] = sfc2 elif sf2 < 1.5 and 1.5 <= sf3: a_pks[num-1] = pk3 a_sfs[num-1] = sf3 a_stiff = a_stiff+['interval=1000mm'] a_scr[num-1] = scr3 a_scc[num-1] = sc3 a_sfc[num-1] = sfc3 ww = 0.25*np.pi*((d0+2*a_t0)**2-d0**2)/1000/1000*al*7.85 # weight print('Weight:{0:10.3f} ton'.format(np.sum(ww))) print('Average thickness: {0:6.1f} mm'.format(np.sum(a_t0*al)/np.sum(al))) df1 = pd.DataFrame({ 'No': no, 'Pi(MPa)': pi, 'L(m)': al, 'D0(mm)': d0, 't0(mm)': a_t0, 'steel': a_sname, 'Eg(MPa)': egg, 'lam': a_lam, 'sig(MPa)': a_sig, 'siga(MPa)': a_sa, 'weight(t)': ww, 'Remarks': df0['Remarks'] }) df2 = pd.DataFrame({ 'No': no, 'Pe(MPa)': pe, 'L(m)': al, 'D0(mm)': d0, 't0(mm)': a_t0, 'steel': a_sname, 'pk_0(MPa)': a_pk0, 'SF_0': a_sf0, 'pk_s(MPa)': a_pks, 'SF_s': a_sfs, 'sc_r(MPa)': a_scr, 'sc_c(MPa)': a_scc, 'SF_c': a_sfc, 'stiffener': a_stiff }) with pd.ExcelWriter('_df_ps.xlsx') as writer: df0.to_excel(writer, sheet_name='Load', index=False) df1.to_excel(writer, sheet_name='Pin', index=False) df2.to_excel(writer, sheet_name='Pex', index=False) if __name__ == '__main__': main()
エクセル書出プログラム
# ================================ # Makeing formated excel file # ================================ import numpy as np import pandas as pd import openpyxl from openpyxl.styles.borders import Border, Side from openpyxl.styles import PatternFill from openpyxl.styles.fonts import Font def xlline(ws, row_m, col_m): # border line bcc = Border(top=Side(style='dashed', color='000000'), bottom=Side(style='dashed', color='000000'), left=Side(style='dashed', color='000000'), right=Side(style='dashed', color='000000') ) for i in range(0, row_m): for j in range(0, col_m): ws.cell(row=i+1, column=j+1).border = bcc def xlswrite(fnameW, df0, df1, df2): # Write results to Excel wb = openpyxl.Workbook() row_height = 20 font = Font(name='Courier New', size=12, bold=False, italic=False, underline='none', strike=False, color='000000') # # ws = wb.create_sheet(index=1, title='Load') df = df0 pcol = [ ['No', '0', 'A', 12], ['L(m)', '0.000', 'B', 12], ['Sum(L)', '0.000', 'C', 12], ['EL(m)', '0.000', 'D', 12], ['D0(m)', '0.000', 'E', 12], ['GWL(m)', '0.000', 'F', 12], ['Hst(m)', '0.000', 'G', 12], ['Hsg(m)', '0.000', 'H', 12], ['Hwh(m)', '0.000', 'I', 12], ['Hin(m)', '0.000', 'J', 12], ['Hex(m)', '0.000', 'K', 12], ['Dr(m)', '0.000', 'L', 12], ['Eg(MPa)', '0', 'M', 12], ['Remarks', '', 'N', 30] ] for j in range(0, len(pcol)): ws.cell(row=1, column=j+1).value = pcol[j][0] ws.cell(row=1, column=j+1).font = font for j in range(0, len(pcol)): temp = np.array(df[pcol[j][0]]) for i in range(len(temp)): ws.cell(row=i+2, column=j+1).number_format = pcol[j][1] ws.cell(row=i+2, column=j+1).value = temp[i] ws.cell(row=i+2, column=j+1).font = font # legend snote = [ 'No: section number', 'L: length of pipe section', 'Sum(L): cumulative length', 'EL: elevation of pipe section', 'D0: internal diameter of pipe section', 'GWL: ground water level', 'Hst: hydrostatic head', 'Hsg: pressure head rise by surging', 'Hwh: pressure head rise by water hammering', 'Hin: design internal pressure head', 'Hex: design external pressure head', 'Dr: excavation diameter of tunnel', 'Eg: elastic modulus of bedrock' ] k = 1+len(temp)+1 for i, ss in enumerate(snote): ws.cell(row=k+i, column=1).value = ss ws.cell(row=k+i, column=1).font = font # cell height and cell width row_max = k+len(snote) col_max = len(pcol) for i in range(0, row_max): ws.row_dimensions[i+1].height = row_height for j in range(0, col_max): ws.column_dimensions[pcol[j][2]].width = pcol[j][3] # line and font xlline(ws, row_max-len(snote)-1, col_max) # # ws = wb.create_sheet(index=1, title='Pin') df = df1 pcol = [ ['No', '0', 'A', 12], ['Pi(MPa)', '0.000', 'B', 12], ['L(m)', '0.000', 'C', 12], ['D0(mm)', '0', 'D', 12], ['t0(mm)', '0', 'E', 12], ['steel', '', 'F', 12], ['Eg(MPa)', '0', 'G', 12], ['lam', '0.000', 'H', 12], ['sig(MPa)', '0.000', 'I', 12], ['siga(MPa)', '0', 'J', 12], ['weight(t)', '0.000', 'K', 12], ['Remarks', '', 'L', 30] ] for j in range(0, len(pcol)): ws.cell(row=1, column=j+1).value = pcol[j][0] ws.cell(row=1, column=j+1).font = font for j in range(0, len(pcol)): temp = np.array(df[pcol[j][0]]) for i in range(len(temp)): ws.cell(row=i+2, column=j+1).number_format = pcol[j][1] ws.cell(row=i+2, column=j+1).value = temp[i] ws.cell(row=i+2, column=j+1).font = font k = 1+len(temp)+1 totalw = np.sum(np.array(df['weight(t)'])) ws.cell(row=k, column=1).number_format = '' ws.cell(row=k, column=1).value = 'Sum' ws.cell(row=k, column=1).font = font ws.cell(row=k, column=11).number_format = '0.000' ws.cell(row=k, column=11).value = totalw ws.cell(row=k, column=11).font = font # legend snote = [ 'No: section number', 'Pi: design internal pressure', 'L: length of pipe section', 'D0: internal diameter of pipe section', 't0: plate thickness of pipe section', 'steel: kind of steel', 'Eg: elastic modulus of bedrock', 'lam: internal pressure shared ratio be bedrock', 'sig: stress of pipe shell', 'siga: allowable stress', 'Weight: weight of pipe section' ] k = k+1 for i, ss in enumerate(snote): ws.cell(row=k+i, column=1).value = ss ws.cell(row=k+i, column=1).font = font # cell height and cell width row_max = k+len(snote) col_max = len(pcol) for i in range(0, row_max): ws.row_dimensions[i+1].height = row_height for j in range(0, col_max): ws.column_dimensions[pcol[j][2]].width = pcol[j][3] # line xlline(ws, row_max-len(snote)-1, col_max) # # ws = wb.create_sheet(index=1, title='Pex') df = df2 pcol = [ ['No', '0', 'A', 12], ['Pe(MPa)', '0.000', 'AB', 12], ['L(m)', '0.000', 'C', 12], ['D0(mm)', '0', 'D', 12], ['t0(mm)', '0', 'E', 12], ['steel', '', 'F', 12], ['pk_0(MPa)', '0.000', 'G', 12], ['SF_0', '0.000', 'H', 12], ['pk_s(MPa)', '0.000', 'I', 12], ['SF_s', '0.000', 'J', 12], ['sc_r(MPa)', '0.000', 'K', 12], ['sc_c(MPa)', '0.000', 'L', 12], ['SF_c', '0.000', 'M', 12], ['stiffener', '', 'N', 30] ] for j in range(0, len(pcol)): ws.cell(row=1, column=j+1).value = pcol[j][0] ws.cell(row=1, column=j+1).font = font for j in range(0, len(pcol)): temp = np.array(df[pcol[j][0]]) for i in range(len(temp)): ws.cell(row=i+2, column=j+1).number_format = pcol[j][1] ws.cell(row=i+2, column=j+1).value = temp[i] ws.cell(row=i+2, column=j+1).font = font # legend snote = [ 'No: section number', 'Pe: design external pressure', 'L: length of pipe section', 'D0: internal diameter of pipe section', 't0: plate thickness of pipe section', 'steel: kind of steel', 'pk_0: critical buckling pressure of pipe section without stiffener', 'SF_0: safety factor against external pressure without pressure', 'pk_s: critical buckling pressure of pipe section with stiffener', 'SF_s: safety factor against external pressure with stiffener', 'sc_r: critical buckling stress of stiffener', 'sc_c: compressive stress of stiffener', 'SF_c: safety factor against compressive stress in stiffener', 'size of stiffener plate: 75mm height x 20mm thickness' ] k = 1+len(temp)+1 for i, ss in enumerate(snote): ws.cell(row=k+i, column=1).value = ss ws.cell(row=k+i, column=1).font = font # cell height and cell width row_max = k+len(snote) col_max = len(pcol) for i in range(0, row_max): ws.row_dimensions[i+1].height = row_height for j in range(0, col_max): ws.column_dimensions[pcol[j][2]].width = pcol[j][3] # line xlline(ws, row_max-len(snote)-1, col_max) # # wb.save(fnameW) def main(): fnameR = '_df_ps.xlsx' fnameW = 'out_penstock.xlsx' sn0 = 'Load' sn1 = 'Pin' sn2 = 'Pex' df0 = pd.read_excel(fnameR, sheet_name=sn0, header=0) df1 = pd.read_excel(fnameR, sheet_name=sn1, header=0) df2 = pd.read_excel(fnameR, sheet_name=sn2, header=0) xlswrite(fnameW, df0, df1, df2) if __name__ == '__main__': main()
以 上
Python 軸対称応力解析プログラム
はじめに
軸対称モデルとして扱う構造物として、水平トンネルなど水平方向に回転軸を持つ構造物と、調圧水槽など鉛直方向に回転軸を持つ構造物がある。 これまで、四角形要素の2次元応力解析プログラムを軸対称解析用に書き換えたプログラムを使ってきたが、座標系(座標軸の方向)を意識せず回転軸を鉛直にして計算を行ったところ、荷重に対して変位及び応力が本来あるべき方向と逆方向の解が計算されてきた。 このため、回転軸(z軸)が水平の場合と、回転軸が鉛直の場合で、場合分けして解析するようプログラムの見直しを行った。
はじめから下図Case-2の座標系を反時計回りに90度回転させてz軸を上向きとした座標を与えればいいのだが、実際のメッシュ切では、2次元応力解析用のメッシャーで作成した要素座標を使うのが便利なため、今回の対応策を導入することにしたものである。
プログラムの改良
以下に座標系と要素のイメージ図を示す。
Case-1およびCase-2では、回転軸をz軸、半径方向軸をr軸としている。
Reference
Reference は通常の2次元応力解析用の座標である。x軸を水平右向き、y軸を鉛直上向きに設定している。 この場合、通常、要素番号は、節点番号を反時計回りに指定して要素を定義する。
Case-1
Case-1は、z軸(回転軸)を水平右向き、r軸(半径方向)を鉛直上向きに設定している。 この場合も、Referenceの2次元応力解析と同様に、要素番号は、節点番号を反時計回りに指定して要素を定義する。(そうするようにプログラムが組まれている)
Case-2
Case-2は、z軸(回転軸)を鉛直上向き、r軸(半径方向)を水平右向きに設定している。 この場合、Reference の2次元応力解析と同様に反時計回りの節点番号で要素を定義するすると、Case-1と比べてz軸とr軸が入れ替わっているため、要素の節点番号が逆周りで定義されたことになり、剛性マトリックスや物体力ベクトル、温度荷重ベクトルなど要素座標を用いて計算する行列やベクトルが通常と逆符号で計算されることとなる。
対応策
このため、今回整備した軸対称応力解析プログラムでは、変数 nzdir
を導入し、座標の方向に応じて、剛性マトリックス、物体力ベクトル、温度荷重ベクトルの方向を逆転させることにした。
すなわち、上図Case-1では、nzdir=1
、上図Case-2ではnzdir=-1
とし、メインプログラム上で、求まった剛性マトリックス、物体力ベクトル、温度荷重ベクトルの符号を調整するようにしている。
sm=sm*float(nzdir) # 剛性マトリックス tfe=tfe*float(nzdir) # 温度荷重ベクトル bfe=bfe*float(nzdir) # 物体力(慣性力)ベクトル
入力データ・フォーマット
npoin nele nsec npfix nlod nzdir # Basic values for analysis E po alpha gamma gkz # Material properties ..... (1 to nsec) ..... # node1 node2 node3 node4 isec # Element connectivity, material set number ..... (1 to nele) ..... # (counter-clockwise order of node number) z r deltaT # Coordinate, Temperature change of node ..... (1 to npoin) ..... # node koz kor rdisz rdisr # Restricted node and restrict conditions ..... (1 to npfix) ..... # node fz fr # Loaded node and loading conditions ..... (1 to nlod) ..... # (omit data input if nlod=0)
変数 | 説明 |
---|---|
npoin, nele, nsec | 節点数,要素数,材料特性数 |
npfix, nlod | 拘束節点数,載荷節点数 |
nzdir | z軸方向(水平右向き: 1,鉛直上向き: −1) |
E, po, alpha | 弾性係数,ポアソン比,線膨張係数 |
gamma, gkz | 単位体積重量,z方向加速度(gの比) |
node1, node2, node3, node4, isec | 節点1,節点2,節点3,節点4,材料特性番号 |
z, r, deltaT | 節点z座標,節点r座標,節点温度変化 |
node, koz, kor | 拘束節点番号,z及びr方向拘束の有無(拘束: 1,自由: 0) |
rdisz, rdisr | z及びr方向変位(無拘束でも0を入力) |
node, fz, fr | 載荷重節点番号,z方向荷重,r方向荷重 |
(注)回転軸zの方向が水平(nzdir=1)であっても鉛直(nzdir=-1)であっても、座標入力は (z , r) の順に行う。
プログラム全文
################################ # Axisymmetric Stress Analysis ################################ import numpy as np from scipy.sparse.linalg import spsolve from scipy.sparse import csr_matrix import sys import time def inpdata_ax4(fnameR,nod,nfree): f=open(fnameR,'r') text=f.readline() text=text.strip() text=text.split() npoin =int(text[0]) # Number of nodes nele =int(text[1]) # Number of elements nsec =int(text[2]) # Number of sections npfix =int(text[3]) # Number of restricted nodes nlod =int(text[4]) # Number of loaded nodes nzdir =int(text[5]) # direction of z-axis (=1 or =-1) # array declanation ae =np.zeros([5,nsec],dtype=np.float64) # Section characteristics node =np.zeros([nod+1,nele],dtype=np.int) # Node-element relationship x =np.zeros([2,npoin],dtype=np.float64) # Coordinates of nodes deltaT=np.zeros(npoin,dtype=np.float64) # Temperature change of node mpfix =np.zeros([nfree,npoin],dtype=np.int) # Ristrict conditions rdis =np.zeros([nfree,npoin],dtype=np.float64) # Ristricted displacement fp =np.zeros(nfree*npoin,dtype=np.float64) # External force vector # section characteristics for i in range(0,nsec): text=f.readline() text=text.strip() text=text.split() ae[0,i]=float(text[0]) #E : Elastic modulus ae[1,i]=float(text[1]) #po : Poisson's ratio ae[2,i]=float(text[2]) #alpha: Thermal expansion coefficient ae[3,i]=float(text[3]) #gamma: Unit weight of material ae[4,i]=float(text[4]) #gkz : Acceleration in z-direction # element-node for i in range(0,nele): text=f.readline() text=text.strip() text=text.split() node[0,i]=int(text[0]) #node_1 node[1,i]=int(text[1]) #node_2 node[2,i]=int(text[2]) #node_3 node[3,i]=int(text[3]) #node_4 node[4,i]=int(text[4]) #section characteristic number # node coordinates for i in range(0,npoin): text=f.readline() text=text.strip() text=text.split() x[0,i] =float(text[0]) # z-coordinate x[1,i] =float(text[1]) # r-coordinate deltaT[i]=float(text[2]) # Temperature change of node # boundary conditions (0:free, 1:restricted) if 0<npfix: for i in range(0,npfix): text=f.readline() text=text.strip() text=text.split() lp=int(text[0]) #fixed node mpfix[0,lp-1]=int(text[1]) #fixed in z-direction mpfix[1,lp-1]=int(text[2]) #fixed in r-direction rdis[0,lp-1]=float(text[3]) #fixed displacement in z-direction rdis[1,lp-1]=float(text[4]) #fixed displacement in r-direction # load if 0<nlod: for i in range(0,nlod): text=f.readline() text=text.strip() text=text.split() lp=int(text[0]) #loaded node fp[2*lp-2]=float(text[1]) #load in z-direction fp[2*lp-1]=float(text[2]) #load in r-direction f.close() return npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp def prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node): fout=open(fnameW,'w') # print out of input data print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}' .format('npoin','nele','nsec','npfix','nlod','nzdir'),file=fout) print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}' .format(npoin,nele,nsec,npfix,nlod,nzdir),file=fout) print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s}' .format('sec','E','po','alpha','gamma','gkz'),file=fout) for i in range(0,nsec): print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}' .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i]),file=fout) print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}' .format('node','z','r','fz','fr','deltaT','koz','kor'),file=fout) for i in range(0,npoin): lp=i+1 print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d}' .format(lp,x[0,i],x[1,i],fp[2*i],fp[2*i+1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout) if 0<npfix: print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}' .format('node','koz','kor','rdis_z','rdis_r'),file=fout) for i in range(0,npoin): if 0<mpfix[0,i]+mpfix[1,i]: lp=i+1 print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}' .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout) print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}' .format('elem','i','j','k','l','sec'),file=fout) for ne in range(0,nele): print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}' .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne],node[4,ne]),file=fout) fout.close() def prout_ax4(fnameW,npoin,nele,disg,strs): fout=open(fnameW,'a') # displacement print('{0:>5s} {1:>15s} {2:>15s}'.format('node','dis-z','dis-r'),file=fout) for i in range(0,npoin): lp=i+1 print('{0:5d} {1:15.7e} {2:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1]),file=fout) # section force print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}' .format('elem','sig_z','sig_r','sig_t','tau_zr','p1','p2','ang'),file=fout) for ne in range(0,nele): sigz =strs[0,ne] sigr =strs[1,ne] sigt =strs[2,ne] tauzr=strs[3,ne] ps1,ps2,ang=pst_ax4(sigz,sigr,tauzr) print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}' .format(ne+1,sigz,sigr,sigt,tauzr,ps1,ps2,ang),file=fout) fout.close() def dmat_ax4(E,po): d=np.zeros([4,4],dtype=np.float64) d[0,0]=1.0-po; d[0,1]=po ; d[0,2]=po ; d[0,3]=0.0 d[1,0]=po ; d[1,1]=1.0-po; d[1,2]=po ; d[1,3]=0.0 d[2,0]=po ; d[2,1]=po ; d[2,2]=1.0-po; d[2,3]=0.0 d[3,0]=0.0 ; d[3,1]=0.0 ; d[3,2]=0.0 ; d[3,3]=0.5*(1.0-2.0*po) d=d*E/(1.0+po)/(1.0-2.0*po) return d def bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4): bm=np.zeros([4,8],dtype=np.float64) #[N] nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 #dN/da,dN/db dn1a=-0.25*(1.0-b) dn2a= 0.25*(1.0-b) dn3a= 0.25*(1.0+b) dn4a=-0.25*(1.0+b) dn1b=-0.25*(1.0-a) dn2b=-0.25*(1.0+a) dn3b= 0.25*(1.0+a) dn4b= 0.25*(1.0-a) #Jacobi matrix and det(J) J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4 J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4 J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4 J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4 detJ=J11*J22-J12*J21 #[B]=[dN/dx][dN/dy] bm[0,0]= J22*dn1a-J12*dn1b bm[0,2]= J22*dn2a-J12*dn2b bm[0,4]= J22*dn3a-J12*dn3b bm[0,6]= J22*dn4a-J12*dn4b bm[1,1]=-J21*dn1a+J11*dn1b bm[1,3]=-J21*dn2a+J11*dn2b bm[1,5]=-J21*dn3a+J11*dn3b bm[1,7]=-J21*dn4a+J11*dn4b bm[2,1]=nm1/rr*detJ bm[2,3]=nm2/rr*detJ bm[2,5]=nm3/rr*detJ bm[2,7]=nm4/rr*detJ bm[3,0]=-J21*dn1a+J11*dn1b bm[3,1]= J22*dn1a-J12*dn1b bm[3,2]=-J21*dn2a+J11*dn2b bm[3,3]= J22*dn2a-J12*dn2b bm[3,4]=-J21*dn3a+J11*dn3b bm[3,5]= J22*dn3a-J12*dn3b bm[3,6]=-J21*dn4a+J11*dn4b bm[3,7]= J22*dn4a-J12*dn4b bm=bm/detJ return bm,detJ def ab_ax4(kk): if kk==0: a=-0.5773502692 b=-0.5773502692 if kk==1: a= 0.5773502692 b=-0.5773502692 if kk==2: a= 0.5773502692 b= 0.5773502692 if kk==3: a=-0.5773502692 b= 0.5773502692 return a,b def sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4): sm=np.zeros([8,8],dtype=np.float64) #Stiffness matrix [sm]=[B]T[D][B]*t*det(J) d=dmat_ax4(E,po) for kk in range(0,4): a,b=ab_ax4(kk) bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4) nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 sm=sm+np.dot(bm.T,np.dot(d,bm))*rr*detJ return sm def calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4): eps0 =np.zeros(4,dtype=np.float64) stress=np.zeros(4,dtype=np.float64) #stress vector {stress}=[D][B]{u} d=dmat_ax4(E,po) for kk in range(0,4): a,b=ab_ax4(kk) bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4) eps=np.dot(bm,wd) # Strain due to temperature change nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4 #Thermal strain eps0[0]=alpha*tem eps0[1]=alpha*tem eps0[2]=alpha*tem eps0[3]=0.0 stress=stress+np.dot(d,(eps-eps0)) stress=0.25*stress return stress def pst_ax4(sigz,sigr,tauzr): ps1=0.5*(sigz+sigr)+np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr) ps2=0.5*(sigz+sigr)-np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr) if sigz==sigr: if tauzr >0.0: ang= 45.0 if tauzr <0.0: ang=135.0 if tauzr==0.0: ang= 0.0 else: ang=0.5*np.arctan(2.0*tauzr/(sigz-sigr)) ang=180.0/np.pi*ang if sigz>sigr and tauzr>=0.0: ang=ang if sigz>sigr and tauzr <0.0: ang=ang+180.0 if sigz<sigr: ang=ang+90.0 return ps1,ps2,ang def tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4): eps0=np.zeros(4,dtype=np.float64) tfe =np.zeros(8,dtype=np.float64) # {tfe=[B]T[D]{eps0} d=dmat_ax4(E,po) for kk in range(0,4): a,b=ab_ax4(kk) bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4) # Strain due to temperature change nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4 #Thermal strain eps0[0]=alpha*tem eps0[1]=alpha*tem eps0[2]=alpha*tem eps0[3]=0.0 rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 tfe=tfe+np.dot(bm.T,np.dot(d,eps0))*rr*detJ return tfe def bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4): ntn=np.zeros([8,8],dtype=np.float64) nm =np.zeros([2,8],dtype=np.float64) bfe=np.zeros(8,dtype=np.float64) for kk in range(0,4): a,b=ab_ax4(kk) nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 dn1a=-0.25*(1.0-b) dn2a= 0.25*(1.0-b) dn3a= 0.25*(1.0+b) dn4a=-0.25*(1.0+b) dn1b=-0.25*(1.0-a) dn2b=-0.25*(1.0+a) dn3b= 0.25*(1.0+a) dn4b= 0.25*(1.0-a) #Jacobi matrix and det(J) J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4 J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4 J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4 J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4 detJ=J11*J22-J12*J21 nm[0,0]=nm1 nm[0,2]=nm2 nm[0,4]=nm3 nm[0,6]=nm4 nm[1,1]=nm[0,0] nm[1,3]=nm[0,2] nm[1,5]=nm[0,4] nm[1,7]=nm[0,6] ntn=ntn+np.dot(nm.T,nm)*gamma*rr*detJ w=np.array([gkz,0,gkz,0,gkz,0,gkz,0],dtype=np.float64) bfe=np.dot(ntn,w) return bfe def main_ax4(): start=time.time() args = sys.argv fnameR=args[1] # input data file fnameW=args[2] # output data file nod=4 # Number of nodes per element nfree=2 # Degree of freedom per node # data input npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_ax4(fnameR,nod,nfree) # print out of input data prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node) # array declanation ir =np.zeros(nod*nfree,dtype=np.int) # Work vector for matrix assembly gk =np.zeros([nfree*npoin,nfree*npoin],dtype=np.float64) # Global stiffness matrix # assembly of stiffness matrix & load vector for ne in range(0,nele): i=node[0,ne]-1 j=node[1,ne]-1 k=node[2,ne]-1 l=node[3,ne]-1 m=node[4,ne]-1 z1=x[0,i]; r1=x[1,i] z2=x[0,j]; r2=x[1,j] z3=x[0,k]; r3=x[1,k] z4=x[0,l]; r4=x[1,l] E =ae[0,m] #E : Elastic modulus po =ae[1,m] #po : Poisson's ratio alpha =ae[2,m] #alpha: Thermal expansion coefficient gamma =ae[3,m] #gamma: Unit weight of material gkz =ae[4,m] #gkz : Acceleration in z-direction dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] # temperature change sm=sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4) # element stiffness matrix tfe=tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4) # termal load vector bfe=bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4) # body force vector sm=sm*float(nzdir) tfe=tfe*float(nzdir) bfe=bfe*float(nzdir) ir[7]=2*l+1; ir[6]=ir[7]-1 ir[5]=2*k+1; ir[4]=ir[5]-1 ir[3]=2*j+1; ir[2]=ir[3]-1 ir[1]=2*i+1; ir[0]=ir[1]-1 for i in range(0,nod*nfree): it=ir[i] fp[it]=fp[it]+tfe[i]+bfe[i] for j in range(0,nod*nfree): jt=ir[j] gk[it,jt]=gk[it,jt]+sm[i,j] # treatment of boundary conditions for i in range(0,npoin): for j in range(0,nfree): if mpfix[j,i]==1: iz=i*nfree+j fp[iz]=0.0 for i in range(0,npoin): for j in range(0,nfree): if mpfix[j,i]==1: iz=i*nfree+j fp=fp-rdis[j,i]*gk[:,iz] gk[:,iz]=0.0 gk[iz,iz]=1.0 # solution of simultaneous linear equations #disg = np.linalg.solve(gk, fp) gk = csr_matrix(gk) disg = spsolve(gk, fp, use_umfpack=True) # recovery of restricted displacements for i in range(0,npoin): for j in range(0,nfree): if mpfix[j,i]==1: iz=i*nfree+j disg[iz]=rdis[j,i] # calculation of element stress wd =np.zeros(8,dtype=np.float64) strs =np.zeros([4,nele],dtype=np.float64) # Section force vector for ne in range(0,nele): i=node[0,ne]-1 j=node[1,ne]-1 k=node[2,ne]-1 l=node[3,ne]-1 m=node[4,ne]-1 z1=x[0,i]; r1=x[1,i] z2=x[0,j]; r2=x[1,j] z3=x[0,k]; r3=x[1,k] z4=x[0,l]; r4=x[1,l] wd[0]=disg[2*i]; wd[1]=disg[2*i+1] wd[2]=disg[2*j]; wd[3]=disg[2*j+1] wd[4]=disg[2*k]; wd[5]=disg[2*k+1] wd[6]=disg[2*l]; wd[7]=disg[2*l+1] E =ae[0,m] po =ae[1,m] alpha=ae[2,m] dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] strs[:,ne]=calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4) # print out of result prout_ax4(fnameW,npoin,nele,disg,strs) # information dtime=time.time()-start print('n={0} time={1:.3f}'.format(nfree*npoin,dtime)+' sec') fout=open(fnameW,'a') print('n={0} time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout) fout.close() #============== # Execution #============== if __name__ == '__main__': main_ax4()
以 上
トランジション(矩形断面から円形断面)での損失水頭
圧力水路の設計をしていると、取水口周辺で矩形断面から円形断面に変化する部分が必要となるが、その時の形状変化による損失(摩擦損失を除く)の計算について書いておく。 これまでは、矩形断面の断面積に等しい等価円を考え、漸縮管の損失水頭として扱ってきた。
今回同じ損失水頭計算が必要となり、しらべたら「矩形管の等価円」という考え方があったので、それについても記載しておく。
等価円の計算法は以下のページによる、
http://www.mekatoro.net/digianaecatalog/ebara-fan50/book/ebara-fan50-P0380.pdf
上記ページによれば等価円の算定式は以下の通り。
上記式および慣用的に用いてきた矩形断面と同一断面積を有する円の直径を計算してみる。
Upstream edge | Length | Downstream edge |
---|---|---|
Rectangular shape (a x b = 8.4m x 8.4m) |
L=5.5m | Circular shape (diameter D=8.4m) |
換算式を用いた等価半径
矩形断面と同一断面積を有する円管の直径
上記より、矩形断面と同一断面積を有する円管の直径は以下の通りとなる。
損失水頭計算
この事例では、矩形断面から換算された円形断面の直径は、2つの手法でほぼ同じとなった。 よってここでは、これまで慣用的に用いてきた、矩形断面と同一断面積を有する円の直径を用いて、漸縮管としての損失水頭計算を行う。
漸縮管の損失水頭は、以下に示すGardelの図(出典:発電水力演習、千秋信一著)による。
計算条件を整理すると以下のとおり。ここで設計流量はQ=290m3/sとする。
Q (m3/s) | D1 (m) | A1 (m2) | v1 (m/s) | D2 (m) | A2 (m2) | v2 (m/s) | A2/A1 | L (m) | 290 | 9.481 | 70.560 | 4.110 | 8.400 | 55.390 | 5.236 | 0.785 | 5.5 |
---|
角度 は、次式で計算する。
断面積比 、角度 より、損失係数は、上記図より となる。 よって、損失水頭は以下の通り計算できる。
以 上