均匀温度场网壳计算

FINISH
/CLEAR
/FILENAME,K6analysis
*AFUN,DEG ! 在角度函数的输入与输出中使用度为单位
/PREP7 ! 进入前处理模块PREP7
!输入基本几何参数
f=8 ! 输入矢高f
Span=40 ! 输入跨度Span
Kn=8 !输入沿环向分割的份数Kn
Nx=6 !输入环杆圈数Nx

!定义单元类型及实常数
ET,1,beam188
SECTYPE,1,BEAM,CTUBE
SECDATA,0.071,0.076
!为材料特性指定温度列表中11个位置处的温度值
MPTEMP,1,20,100,200,300,400,500
MPTEMP,7,600,700,800,900,1000
!指定各个温度下钢管材料的弹性模量
MPDATA,EX,1,1,2.065E+11,2.065E+11,1.853E+11,1.649E+11,1.450E+11,1.240E+11
MPDATA,EX,1,7,6.389E+10,2.679E+10,1.858E+10,1.389E+10,9.27E+09
!指定各个温度下钢管材料的热膨胀系数
MPDATA,ALPX,1,1,1.4E-05,1.4E-05,1.4E-05,1.4E-05,1.4E-05,1.4E-05
MPDATA,ALPX,1,7,1.4E-05,1.4E-05,1.4E-05,1.4E-05,1.4E-05
!指定各个温度下钢管材料的密度
MPDATA,DENS,1,1,7850,7850,7850,7850,7850,7850
MPDATA,DENS,1,7,7850,7850,7850,7850,7850
!指定各个温度下钢管材料的泊松比
MPDATA,PRXY,1,1,0.3,0.3,0.3,0.3,0.3,0.3
MPDATA,PRXY,1,7,0.3,0.3,0.3,0.3,0.3
!指定各个温度下钢管材料的参考温度
MPDATA,REFT,1,1,20,20,20,20,20,20
MPDATA,REFT,1,7,20,20,20,20,20
!指定各个温度下钢管材料的应力应变曲线，采用多线段非线性材料模型
!定义多线性钢材的应力应变模型参数，每个温度值下的应力应变曲线采用8个点近似模拟
TB,MISO,1,11,7
TBTEM,20,1 !定义温度为20时钢管的应力应变曲线
TBPT,DEFI,1.67E-03,3.45E+08
TBPT,DEFI,1.79E-03,3.45E+08
TBPT,DEFI,1.99E-03,3.45E+08
TBPT,DEFI,2.30E-03,3.45E+08
TBPT,DEFI,2.62E-03,3.45E+08
TBPT,DEFI,2.94E-02,3.45E+08
TBPT,DEFI,0.15,3.45E+08
TBTEM,100,2 !定义温度为100时钢管的应力应变曲线
TBPT,DEFI,1.67E-03,3.45E+08
TBPT,DEFI,1.79E-03,3.45E+08
TBPT,DEFI,1.99E-03,3.45E+08
TBPT,DEFI,2.30E-03,3.45E+08
TBPT,DEFI,2.62E-03,3.45E+08
TBPT,DEFI,2.94E-02,3.45E+08
TBPT,DEFI,0.15,3.45E+08
TBTEM,200,3 !定义温度为200时钢管的应力应变曲线
TBPT,DEFI,1.50E-03,2.78E+08
TBPT,DEFI,2.37E-03,2.94E+08
TBPT,DEFI,4.98E-03,3.10E+08
TBPT,DEFI,8.46E-03,3.22E+08
TBPT,DEFI,1.54E-02,3.36E+08
TBPT,DEFI,2.94E-02,3.38E+08
TBPT,DEFI,0.15,3.45E+08
TBTEM,300,4 !定义温度为300时钢管的应力应变曲线
TBPT,DEFI,1.28E-03,2.11E+08
TBPT,DEFI,2.37E-03,2.45E+08
TBPT,DEFI,4.83E-03,2.74E+08
TBPT,DEFI,8.37E-03,2.99E+08
TBPT,DEFI,2.23E-02,3.41E+08
TBPT,DEFI,2.94E-02,3.42E+08
TBPT,DEFI,0.15,3.45E+08
TBTEM,400,5 !定义温度为400时钢管的应力应变曲线
TBPT,DEFI,1.0E-03,1.45E+08
TBPT,DEFI,2.53E-03,2.03E+08
TBPT,DEFI,8.09E-03,2.74E+08
TBPT,DEFI,1.52E-02,3.17E+08
TBPT,DEFI,2.23E-02,3.38E+08
TBPT,DEFI,2.94E-02,3.41E+08
TBPT,DEFI,0.15,3.45E+08
TBTEM,500,6 !定义温度为500时钢管的应力应变曲线
TBPT,DEFI,1.0E-03,1.24E+08
TBPT,DEFI,

2.53E-03,1.67E+08
TBPT,DEFI,5.15E-03,1.97E+08
TBPT,DEFI,1.52E-02,2.49E+08
TBPT,DEFI,2.23E-02,2.64E+08
TBPT,DEFI,2.94E-02,2.66E+08
TBPT,DEFI,0.15,2.69E+08
TBTEM,600,7 !定义温度为600时钢管的应力应变曲线
TBPT,DEFI,9.72E-04,6.21E+07
TBPT,DEFI,2.75E-03,9.33E+07
TBPT,DEFI,8.07E-03,1.26E+08
TBPT,DEFI,1.52E-02,1.48E+08
TBPT,DEFI,2.23E-02,1.59E+08
TBPT,DEFI,2.94E-02,1.60E+08
TBPT,DEFI,0.15,1.62E+08
TBTEM,700,8 !定义温度为700时钢管的应力应变曲线
TBPT,DEFI,9.66E-04,2.59E+07
TBPT,DEFI,2.50E-03,4.07E+07
TBPT,DEFI,8.06E-03,5.99E+07
TBPT,DEFI,1.52E-02,7.16E+07
TBPT,DEFI,2.23E-02,7.75E+07
TBPT,DEFI,2.94E-02,7.78E+07
TBPT,DEFI,0.15,7.94E+07
TBTEM,800,9 !定义温度为800时钢管的应力应变曲线
TBPT,DEFI,9.31E-04,1.73E+07
TBPT,DEFI,2.48E-03,2.34E+07
TBPT,DEFI,8.04E-03,3.07E+07
TBPT,DEFI,1.52E-02,3.51E+07
TBPT,DEFI,2.23E-02,3.73E+07
TBPT,DEFI,2.94E-02,3.74E+07
TBPT,DEFI,0.15,3.80E+07
TBTEM,900,10 !定义温度为900时钢管的应力应变曲线
TBPT,DEFI,9.29E-04,1.29E+07
TBPT,DEFI,2.48E-03,1.54E+07
TBPT,DEFI,8.04E-03,1.80E+07
TBPT,DEFI,1.52E-02,1.95E+07
TBPT,DEFI,2.23E-02,2.05E+07
TBPT,DEFI,2.94E-02,2.06E+07
TBPT,DEFI,0.15,2.07E+07
TBTEM,1000,11 !定义温度为1000时钢管的应力应变曲线
TBPT,DEFI,9.31E-04,8.63E+06
TBPT,DEFI,2.48E-03,1.02E+07
TBPT,DEFI,8.04E-03,1.20E+07
TBPT,DEFI,1.52E-02,1.31E+07
TBPT,DEFI,2.23E-02,1.36E+07
TBPT,DEFI,2.94E-02,1.37E+07
TBPT,DEFI,0.15,1.38E+07


!计算节点坐标位置,并定义节点
CSYS,2 !转换为球面坐标系
R=(Span*Span/4+f*f)/(2*f) !根据矢高和跨度计算曲面半径R
DPha=Atn(Span/2/Sqrt(R*R-Span*Span/4))/Nx !计算相邻两圈环杆对应的球心夹角Dpha
N,1,R,0,90 ! 定义顶点1号节点的位置坐标为（R,0,90）
*DO,i,1,Nx ! 第1～Nx圈的节点循环
*DO,j,1,Kn*i ! 对称区的节点循环,第i圈可分为Kn*i份
x=R ! 输入x坐标
y=(j-1)*360/(Kn*i) ! 输入y坐标，各圈上的节点数目为Kn倍圈号
z=90-i*DPha ! 输入z坐标
N,1+Kn*(i-1)*i/2+j,x,y,z ! 依次计算并定义节点的编号和位置
*ENDDO
*ENDDO

!定义单元连接
type,1
mat,1
real,1
!环向杆连接
*DO,i,1,Nx ! 第1～Nx圈的节点循环
*DO,j,1,Kn*i-1 ! 第1～Kn*i-1对称区的节点循环
E,1+Kn*(i-1)*i/2+j,1+Kn*(i-1)*i/2+j+1 ! 连接相邻两节点生成单元
*ENDDO
E,1+Kn*(i-1)*i/2+1,1+Kn*(i-1)*i/2+Kn*i ! 连接最后一对称区的环杆
*ENDDO
! 径向的单元连接
*DO,i,1,Kn ! 对各对称区循环
E,1,1+i ! 定义中心处径向杆单元
*ENDDO

*DO,i,1,Nx-1 ! 从里圈第2圈开始向最外圈循环
*DO,j,1,Kn ! 对各个对称区循环
*DO,k,1,i+1 ! 对每个对称区内的杆循环
*IF,k,EQ,i+1,THEN ! 判断是否是当前对称区内的最后一根杆
*IF,j,EQ,K

n,THEN ! 判断是否是最后一个对称区
E,1+Kn*(i-1)*i/2+1,1+Kn*(i+1)*i/2+(j-1)*(i+1)+k ! 第1个节点和最后一个节点的连杆
*ELSE
E,1+Kn*(i-1)*i/2+(j-1)*i+k,1+Kn*(i+1)*i/2+(j-1)*(i+1)+k !一般区内的最后一根连杆
*ENDIF
*ELSE
E,1+Kn*(i-1)*i/2+(j-1)*i+k,1+Kn*(i+1)*i/2+(j-1)*(i+1)+k ! 一般正向连杆
*ENDIF
*ENDDO
*DO,k,1,i
E,1+Kn*(i-1)*i/2+(j-1)*i+k,1+Kn*(i+1)*i/2+(j-1)*(i+1)+k+1 ! 一般负向连杆
*ENDDO
*ENDDO
*ENDDO


















!获取温度值
*dim,t_s,,1440 !定义两个小时内，时间间隔为5s的1440个点的钢管温度数组
*dim,t_deta,,1440!定义两个小时内，1440个时间点的钢管温度增量数组
*dim,t_pipe,,360 !定义两个小时内，时间间隔为20s的360个点的钢管温度数组
desi=7850 !定义钢管材料的密度
f_v=206.8 !定义钢管的截面形状系数
deta_t=5 !定义计算钢管1440个时间点的温度值时的时间间隔为5s
c_f=45 !定义钢管材料的热传到系数
*set,t_s(1),20 !设置钢管温度的初值（即常温下钢管的温度）为20摄氏度
*do,i,2,1440,1 !利用循环计算其余1439个时间点的钢管温度值
!计算当前时刻空气温度值（基于ISO834标准升温曲线）
Tg=20+345*LOG10(8*(i-1)*5/60+1)
!计算以辐射方式由空气向构件表面传热的传热系数
a_r=0.5*5.67e-8*((tg+273)**4-(t_s(i-1)+273)**4)/(tg-t_s(i-1))
!计算综合传热系数
k_f=a_r+25
!计算当前钢管温度的温度增量
*set,t_deta(i-1),k_f*f_v*(tg-t_s(i-1))*deta_t/(desi*c_f)
!计算当前的钢管温度
*set,t_s(i),t_s(i-1)+t_deta(i-1)
*enddo
!提取间隔20s的钢管温度值并赋值给t_pipe数组
*do,i,1,360,1
*set,t_pipe(i),t_s((i-1)*4+1)
*enddo





CSYS,0
/SOLU !进入求解器
!定义边界约束
*DO,i,1, 1+Kn*(Nx-1)*Nx/2+Kn*Nx ! 所有节点循环
*If,i,GT, 1+Kn*(Nx-1)*Nx/2,THEN ! 选择非边界节点
D,i,UX,0 ! 定义约束
D,i,UY,0
D,i,UZ,0 ! 给边界节点定义边界约束
*ENDIF
*ENDDO



allsel
Numnode=1+Kn*(Nx-1)*Nx/2+Nx*Kn
*DIM,Force1,,Numnode,3 !定义节点等效荷载数组
/input,hz1,txt
allsel
*DO,i,1, 1+Kn*(Nx-1)*Nx/2+Kn*Nx
F,i,FX,Force1(i,1)
F,i,FY,Force1(i,2)
F,i,FZ,Force1(i,3)
*ENDDO









ANTYPE,0 ! 设置求解类型为静力分析
PSTRES,ON ! 打开预应力选项
SOLVE ! 静力求解
FINISH ! 退出求解器
!获得特征值屈曲解
/SOLU ! 进入求解器
ANTYPE,1 ! 设置分析类型为特征值屈曲分析
BUCOPT,LANB,10 ! 设置屈曲模态提取方

法及模态提取数
MXPAND,10,0,0,1,0.001 ! 设置屈曲模态扩展数及扩展算法选项
SOLVE ! 特征值屈曲分析
FINISH
/PREP7 !进入前处理模块，考虑初始缺陷
UPGEOM,Span/300,1,1,K6analysis,RST !以一阶模态（即第一个荷载步的第一个子步）位移值的Span/300修改原模型节点坐标
FINISH !退出前处理器
/SOLU
tunif,20
tref,20
ANTYPE,STATIC !指定为静力分析
ERESX,NO
AUTOTS,ON
NLGEOM,ON
OUTRES,BASIC,1
NSUBST,1,500,1
TIME,1
ACEL,0,9.8 !施加重力加速度
*DIM,Force2,,Numnode,3
/input,hz2,txt
allsel
*DO,i,1, 1+Kn*(Nx-1)*Nx/2+Kn*Nx
F,i,FX,Force2(i,1)
F,i,FY,Force2(i,2)
F,i,FZ,Force2(i,3)
*ENDDO
SOLVE !求解
*DO,tm,20,7200,20
time,tm
*do,i,1,1+Kn*(Nx-1)*Nx/2+Kn*Nx
bfe,all,temp,,t_pipe(tm/20)
*enddo
allsel
lswrite,tm/20+1
solve
*ENDDO
finish