Python 岩盤内埋設式水圧鉄管設計プログラム（旧版）
概要
岩盤内埋設式水圧鉄管の設計プログラムをフルで書いてみたので残しておく。
プログラム前半の非常に長い関数「def xlswrite」は、エクセルに書式指定して結果を書き出すもの。 最終的にワードの報告書に貼り付けなくてはならないので、エクセルにしておくと便利である。このため、面倒ではあったがこの部分も作っておいた。
入力
プログラムでは 'load2.xlsx' を読み込んでいるが、これは出力例「荷重計算」と同じ内容（欄外の変数凡例行は除く）のファイルである。
出力
プログラムでは出力は'penstock2.xlsx' というエクセルファイルに行われ、以下のシートに結果が出力される。
| Sheet | Outout |
|---|---|
| Load | 荷重計算（入力ファイルと同じ） |
| Pin | 内圧設計結果 |
| Pex | 外圧設計結果 |
なお設計では、関数「calc1」の中で、水門鉄管技術基準で定める最小板厚とは別に、施工性や耐外厚設計などを考慮し、これより大きな推奨最小板厚（recommended minimum thickness: t0t）というものを導入し、板厚がこれ以下にならないようにしている。
プログラム上での使用鋼材はJIS：SM400, SM490, SM570である。
荷重計算
内圧設計
外圧設計
プログラム
import numpy as np import pandas as pd from scipy import optimize import openpyxl from openpyxl.styles.borders import Border, Side from openpyxl.styles import PatternFill from openpyxl.styles.fonts import Font def xlline(ws, cmax): def wline(ws, row_num, cmax, bl, bc, br): col_num=1; ws.cell(row=row_num ,column=col_num).border = bl for col_num in range(2,cmax): ws.cell(row=row_num ,column=col_num).border = bc col_num=cmax; ws.cell(row=row_num ,column=col_num).border = br # font font = Font(name='Courier New', size=12, bold=False, italic=False, underline='none', strike=False, color='000000') for row_num in range(1,53): for col_num in range(1,cmax+1): ws.cell(row=row_num ,column=col_num).font = font # cell height and width lcol=['A','B','C','D','E','F','G','H','I','J','K','L','M','N'] for row_num in range(1,53): ws.row_dimensions[row_num].height = 20 for i in range(0,cmax): ws.column_dimensions[lcol[i]].width = 12 ws.column_dimensions[lcol[cmax-1]].width = 30 # set color for cell filling fill = PatternFill(patternType='solid', fgColor='ffffaa') # write in sheet for row_num in range(6,12): for col_num in range(1,cmax+1): ws.cell(row=row_num ,column=col_num).fill = fill for row_num in range(20,26): for col_num in range(1,cmax+1): ws.cell(row=row_num ,column=col_num).fill = fill bul = Border(top=Side(style='medium', color='000000'), bottom=Side(style='dashed', color='000000'), left=Side(style='medium', color='000000'), right=Side(style='dashed', color='000000') ) buc = Border(top=Side(style='medium', color='000000'), bottom=Side(style='dashed', color='000000'), left=Side(style='dashed', color='000000'), right=Side(style='dashed', color='000000') ) bur = Border(top=Side(style='medium', color='000000'), bottom=Side(style='dashed', color='000000'), left=Side(style='dashed', color='000000'), right=Side(style='medium', color='000000') ) bcl = Border(top=Side(style='dashed', color='000000'), bottom=Side(style='dashed', color='000000'), left=Side(style='medium', color='000000'), right=Side(style='dashed', color='000000') ) 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') ) bcr = Border(top=Side(style='dashed', color='000000'), bottom=Side(style='dashed', color='000000'), left=Side(style='dashed', color='000000'), right=Side(style='medium', color='000000') ) bll = Border(top=Side(style='dashed', color='000000'), bottom=Side(style='medium', color='000000'), left=Side(style='medium', color='000000'), right=Side(style='dashed', color='000000') ) blc = Border(top=Side(style='dashed', color='000000'), bottom=Side(style='medium', color='000000'), left=Side(style='dashed', color='000000'), right=Side(style='dashed', color='000000') ) blr = Border(top=Side(style='dashed', color='000000'), bottom=Side(style='medium', color='000000'), left=Side(style='dashed', color='000000'), right=Side(style='medium', color='000000') ) row_num=1; wline(ws, row_num, cmax, bul, buc, bur) row_num=2; wline(ws, row_num, cmax, bul, buc, bur) for row_num in range(3,29): wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=29; wline(ws, row_num, cmax, bul, buc, bur) row_num=30; wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=31; wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=32; wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=33; wline(ws, row_num, cmax, bul, buc, bur) row_num=34; wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=35; wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=36; wline(ws, row_num, cmax, bul, buc, bur) row_num=37; wline(ws, row_num, cmax, bcl, bcc, bcr) row_num=38; wline(ws, row_num, cmax, bll, blc, blr) if cmax==11: row_num=39; wline(ws, row_num, cmax, bll, blc, blr) def xlswrite(df0,df1,df2): # Write results to Excel fnameW='penstock2.xlsx' wb = openpyxl.Workbook() # ws=wb.create_sheet(index=1,title='Load') ws.cell(row=1,column=1).value='No' ws.cell(row=1,column=2).value='L(m)' ws.cell(row=1,column=3).value='Sum(L)' ws.cell(row=1,column=4).value='EL(m)' ws.cell(row=1,column=5).value='D0(m)' ws.cell(row=1,column=6).value='GWL(m)' ws.cell(row=1,column=7).value='Hst(m)' ws.cell(row=1,column=8).value='Hsg(m)' ws.cell(row=1,column=9).value='Hwh(m)' ws.cell(row=1,column=10).value='Hin(m)' ws.cell(row=1,column=11).value='Hex(m)' ws.cell(row=1,column=12).value='Remarks' # j=1; temp=np.array(df0['No'],dtype=np.int) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=2; temp=np.array(df0['L(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=3; temp=np.array(df0['Sum(L)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=4; temp=np.array(df0['EL(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=5; temp=np.array(df0['D0(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=6; temp=np.array(df0['GWL(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=7; temp=np.array(df0['Hst(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=8; temp=np.array(df0['Hsg(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=9; temp=np.array(df0['Hwh(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=10; temp=np.array(df0['Hin(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=11; temp=np.array(df0['Hex(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=12; temp=df0['Remarks'] for i in range(len(temp)): ws.cell(row=i+2,column=j).value = temp[i] # border line cmax=12; xlline(ws, cmax) # legend i=38 i=i+1; ws.cell(row=i,column=1).value = 'No: section number' i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'Sum(L): cumulative length' i=i+1; ws.cell(row=i,column=1).value = 'EL: elevation of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'GWL: ground water level' i=i+1; ws.cell(row=i,column=1).value = 'Hst: hydrostatic head' i=i+1; ws.cell(row=i,column=1).value = 'Hsg: pressure head rise by surging' i=i+1; ws.cell(row=i,column=1).value = 'Hwh: pressure head rise by water hammering' i=i+1; ws.cell(row=i,column=1).value = 'Hin: design internal pressure head' i=i+1; ws.cell(row=i,column=1).value = 'Hex: design external pressure head' # ws=wb.create_sheet(index=1,title='Pin') ws.cell(row=1,column=1).value='No' ws.cell(row=1,column=2).value='Pi(MPa)' ws.cell(row=1,column=3).value='L(m)' ws.cell(row=1,column=4).value='D0(mm)' ws.cell(row=1,column=5).value='t0(mm)' ws.cell(row=1,column=6).value='steel' ws.cell(row=1,column=7).value='lam' ws.cell(row=1,column=8).value='sig(MPa)' ws.cell(row=1,column=9).value='siga(MPa)' ws.cell(row=1,column=10).value='Weight(t)' ws.cell(row=1,column=11).value='Remarks' # j=1; temp=np.array(df1['No'],dtype=np.int) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=2; temp=np.array(df1['Pi(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=3; temp=np.array(df1['L(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=4; temp=np.array(df1['D0(mm)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=5; temp=np.array(df1['t0(mm)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=6; temp=df1['steel'] for i in range(len(temp)): ws.cell(row=i+2,column=j).value = temp[i] j=7; temp=np.array(df1['lam'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=8; temp=np.array(df1['sig(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=9; temp=np.array(df1['siga(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=10; temp=np.array(df1['weight(t)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] ws.cell(row=len(temp)+2,column=1).value = 'Sum' ws.cell(row=len(temp)+2,column=10).number_format = "0.000" ws.cell(row=len(temp)+2,column=10).value = np.sum(temp) j=11; temp=df0['Remarks'] for i in range(len(temp)): ws.cell(row=i+2,column=j).value = temp[i] # border line cmax=11; xlline(ws, cmax) # legend i=39 i=i+1; ws.cell(row=i,column=1).value = 'No: section number' i=i+1; ws.cell(row=i,column=1).value = 'Pi: design internal pressure' i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section' i=i+1; ws.cell(row=i,column=1).value = 't0: plate thickness of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'steel: kind of steel' i=i+1; ws.cell(row=i,column=1).value = 'lam: internal pressure shared ratio be bedrock' i=i+1; ws.cell(row=i,column=1).value = 'sig: stress of pipe shell' i=i+1; ws.cell(row=i,column=1).value = 'siga: allowable stress' i=i+1; ws.cell(row=i,column=1).value = 'Weight: weight of pipe section' # ws=wb.create_sheet(index=1,title='Pex') ws.cell(row=1,column=1).value='No' ws.cell(row=1,column=2).value='Pe(MPa)' ws.cell(row=1,column=3).value='L(m)' ws.cell(row=1,column=4).value='D0(mm)' ws.cell(row=1,column=5).value='t0(mm)' ws.cell(row=1,column=6).value='steel' ws.cell(row=1,column=7).value='pk_0(MPa)' ws.cell(row=1,column=8).value='SF_0' ws.cell(row=1,column=9).value='pk_s(MPa)' ws.cell(row=1,column=10).value='SF_s' ws.cell(row=1,column=11).value='sc_r(MPa)' ws.cell(row=1,column=12).value='sc_c(MPa)' ws.cell(row=1,column=13).value='SF_c' ws.cell(row=1,column=14).value='stiffener' # j=1; temp=np.array(df2['No'],dtype=np.int) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=2; temp=np.array(df2['Pe(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=3; temp=np.array(df2['L(m)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=4; temp=np.array(df2['D0(mm)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=5; temp=np.array(df2['t0(mm)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0" ws.cell(row=i+2,column=j).value = temp[i] j=6; temp=df2['steel'] for i in range(len(temp)): ws.cell(row=i+2,column=j).value = temp[i] j=7; temp=np.array(df2['pk0(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=8; temp=np.array(df2['sf0'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=9; temp=np.array(df2['pks(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=10; temp=np.array(df2['sfs'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=11; temp=np.array(df2['scr(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=12; temp=np.array(df2['scc(MPa)'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=13; temp=np.array(df2['sfc'],dtype=np.float64) for i in range(len(temp)): ws.cell(row=i+2,column=j).number_format = "0.000" ws.cell(row=i+2,column=j).value = temp[i] j=14; temp=df2['stiffener'] for i in range(len(temp)): ws.cell(row=i+2,column=j).value = temp[i] # border line cmax=14; xlline(ws, cmax) # legend i=38 i=i+1; ws.cell(row=i,column=1).value = 'No: section number' i=i+1; ws.cell(row=i,column=1).value = 'Pe: design external pressure' i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section' i=i+1; ws.cell(row=i,column=1).value = 't0: plate thickness of pipe section' i=i+1; ws.cell(row=i,column=1).value = 'steel: kind of steel' i=i+1; ws.cell(row=i,column=1).value = 'pk_0: critical buckling pressure of pipe section without stiffener' i=i+1; ws.cell(row=i,column=1).value = 'SF_0: safety factor against external pressure without pressure' i=i+1; ws.cell(row=i,column=1).value = 'pk_s: critical buckling pressure of pipe section with stiffener' i=i+1; ws.cell(row=i,column=1).value = 'SF_s: safety factor against external pressure with stiffener' i=i+1; ws.cell(row=i,column=1).value = 'sc_r: critical buckling stress of stiffener' i=i+1; ws.cell(row=i,column=1).value = 'sc_c: compressive stress of stiffener' i=i+1; ws.cell(row=i,column=1).value = 'SF_c: safety factor against compressive stress in stiffener' i=i+1; ws.cell(row=i,column=1).value = 'size of stiffener plate: 75mm height x 20mm thickness' # wb.save(fnameW) def stcr(dd0,t0,eps,hr,tr,el,pp,sigF): # Calculation of critical stress and compressive stress of stiffener def calrg(dd0,t0,eps,hr,tr): 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): # 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 eps=2.0 # margin of plate thickness Es=206000.0 # elastic modulus of steel plate ns=0.3 # Poisson's ratio of steel plate #D0=2100.0 # design internal diameter of penstock #t0=30.0 # design plate thickness of penstock #sigF=315.0 # yield stress of steel plate (SM490) #siga=175.0 # allowable stress of steel plate (SM490) 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-ns**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-ns**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sinh(beta*el))/(np.cosh(beta*el)-np.cosh(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,eps): # Calculate 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 # t: plate thickness pi=np.pi Es=206000 # elastic modulus of steel ns=0.3 # Poisson's ratio of steel 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-ns**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-ns**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sinh(beta*el))/(np.cosh(beta*el)-np.cosh(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-ns**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-ns)/(1+n**2*ell**2/pi**2/r0d**2)) _pk=(c1+c2)*Es*t/(1-ns**2)/r0d apk=apk+[_pk] pk=min(apk) return pk def ams(dd0,t0,sigF,siga): # Calculate criticsl buckling pressure # of steel pipe without stiffener by Amstutz's formula def func(x,params): # 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): aa=rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess) pk = sigN/aa return pk def cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag): aa = alphas*deltaT+betag*siga*eta/Es k0 = aa*r0d/(1+betag) return k0 # main routine eps=2.0 # margin of plate thickness Es=206000.0 # elastic modulus of steel plate pos=0.3 # Poisson's ratio of steel plate eta=0.85 # welding efficiency alphas=1.2e-5 # thermal expansion coefficient deltaT=20.0 # temperature decrease betag=1.0 # plastic deformation coefficient of bedrock #D0=2100.0 # design internal diameter of penstock #t0=30.0 # design plate thickness of penstock #sigF=315.0 # yield stress of steel plate (SM490) #siga=175.0 # allowable stress of steel plate (SM490) 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 r0d=(dd0+2*t0)/2 k0=cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag) 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 datainp(): df0 = pd.read_excel('load2.xlsx') return df0 def calc1(num,pp,dd0,sa,eps): # Design for internal pressure # Calculate required plate thickness Es=206000.0 alpha=1.2e-5 deltaT=20.0 Ec=20600.0 Bc=0.0 Eg=2000.0 Bg=1.0 mg=4.0 dr=5.2 t0t=25 # 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 5<=num<=10 or 20<=num<=25: 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 main(): df0=datainp() no=np.array(df0['No'],dtype=np.int) al=np.array(df0['L(m)'],dtype=np.float64) # 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) eps=2.0 ef=0.85 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,pp,dd0,ppe in zip(no,pi,d0,pe): sname='SM400' sa=130*ef t0,sig,lam=calc1(num,pp,dd0,sa,eps) if sname=='SM400' and 40<t0: sname='SM490' sa=175*ef t0,sig,lam=calc1(num,pp,dd0,sa,eps) if sname=='SM490' and 40<t0: sname='SM570' sa=240*ef t0,sig,lam=calc1(num,pp,dd0,sa,eps) if sname=='SM570' and 40<t0: sname='SM570' sa=235*ef t0,sig,lam=calc1(num,pp,dd0,sa,eps) 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,sa/ef) sf0=pk0/ppe if sf0<1.5: # stiffener hr=75; tr=20 # stiffener interval=2000mm el=3000.0 pk1=timo(el,hr,tr,dd0,t0,eps) sf1=pk1/ppe scr1,sc1=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF) sfc1=scr1/sc1 # stiffener interval=1500mm el=1500.0 pk2=timo(el,hr,tr,dd0,t0,eps) sf2=pk2/ppe scr2,sc2=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF) sfc2=scr2/sc2 # stiffener interval=1000mm el=1000.0 pk3=timo(el,hr,tr,dd0,t0,eps) sf3=pk3/ppe scr3,sc3=stcr(dd0,t0,eps,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, '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, 'pk0(MPa)': a_pk0, 'sf0': a_sf0, 'pks(MPa)': a_pks, 'sfs': a_sfs, 'scr(MPa)': a_scr, 'scc(MPa)': a_scc, 'sfc': a_sfc, 'stiffener': a_stiff }) with pd.ExcelWriter('penstock.xlsx') as writer: df0.to_excel(writer,sheet_name='Load',index=False) df1.to_excel(writer,sheet_name='Internal Pressure',index=False) df2.to_excel(writer,sheet_name='External Pressure',index=False) xlswrite(df0,df1,df2) if __name__ == '__main__': main()
以 上
