!统一截面LINK33
/prep7
! 定义环境参数
SHPP,ON !显示启用
N = 172 ! 日序数(年中的第85天)
F = 0.4 ! 辐射吸收系数(白色氯化橡胶面漆)
current_time= 13.8 ! 地方时间(小时)
T1 = 31.9 ! 结构表面01:00的结构表面温度(℃)
T2 = 32 ! 大气环境温度
T3 = T1 + 273.15 ! 结构表面开尔文温度
T4 = T2 + 273.15 ! 大气环境开尔文温度
F1 = 0.92 ! 表面辐射率(白色氯化橡胶面漆)
V = 3 ! 风速参量 单位(m/s)若转为毫米制
! 常数
z = 5.667e-8 ! 斯蒂芬-玻尔兹曼常数 单位 米制-8,毫米制-11
D = 0.2 ! 地(水)面反射系数
PI = 3.1415926535
W = 26.04 ! 纬度
P = 0.624 ! 大气透明系数
MP,HF,1,1e-10 !定义热生成率,由于太阳辐射,内部没有生成热,定义为极小
MP,KXX,1,45 !单位 米是45 毫米是0.045
! 创建局部坐标系(绕Z轴顺时针旋转122度)
LOCAL, 11, 0, 0, 0, 0, -122, 0, 0
CSYS, 11
! 验证新正北方向
CSYS, 0
K, 100, 0, 1, 0
CSYS, 11
*GET, KP_Y_DIR, KP, 100, LOC, Y
! 定义单元类型
*DO, I, 1, 1379, 1
ET, I, LINK33
*ENDDO
! 定义统一截面(外径450mm,内径420mm,壁厚15mm)
*SET, OD, 0.45 ! 外径 !单位
*SET, ID, 0.42 ! 内径 !单位
*SET, Area, PI*(OD**2 - ID**2)/4 ! 横截面积
!*SET,RID33,1000
*DIM,R33_ARR,ARRAY,4980
*DO,I,1,4980
*SET,R33_ARR(I),0
*ENDDO
*SET,RID33,0
*DO, ELEM_NUM, 1, 4980
R, RID33, Area ! 仅需传递截面积
REAL,RID33
*SET,R33_ARR(ELEM_NUM), RID33
*SET,RID33, RID33+1
*ENDDO
ESEL,S,ELEM,,1,4980
ESEL,ALL
! 计算换热系数
*SET, H, F1*z*(T3**2 + T4**2)*(T3 + T4)
delta_T = ABS(T1 - T2)
*SET, H1, 2.6*(delta_T**0.25 + 1.54*V)
*SET, H2, H + H1
! 太阳辐射计算
M = 1370*(1 + 0.034*COS( (360*N/365)*PI/180 ))
B1=0.2051 - 4.05369e-4*N + 3.5186e-5*N**2
B2= 2.8939e-10*N**4 - 1.9832e-7*N**3
B=B1+B2
*SET, C1, 0.078763 - 4.2177e-4*N
*SET, C2, 1.9908e-5*N**2 - 1.0607e-7*N**3 + 1.5024e-10*N**4
C=C1+C2
T = (current_time - 12)*15
delta = 23.45*SIN( (360*(284 + N)/365)*PI/180 )
sinh = SIN(W*PI/180)*SIN(delta*PI/180) + COS(W*PI/180)*COS(delta*PI/180)*COS(T*PI/180)
I1 = P*M / EXP(B/sinh)
!*if, I1, gt, 1200, then
! I1 = 0
!*endif
*SET, cosh, SQRT(1 - sinh**2)
*SET, sin_A, (COS(delta*PI/180)*SIN(T*PI/180)) / cosh
*SET, A_rad, ASIN(sin_A)
*SET, A_deg, A_rad*180/PI
N,1000000,0,30000,20000
D,1000000,TEMP,T4
! 定义LINK34单元
ET,5000,LINK34
KEYOPT,5000,1,0
! 定义实常数,尝试初始化为0
!*SET, RID34,2000
*DIM, R34_ARR, ARRAY,4980
*DO,I,1,4980
*SET,R34_ARR(I),0
*ENDDO
*DIM, NODE_AREA, ARRAY, 200000
*VOPER, NODE_AREA(1), NODE_AREA(1), ADD, 0
*DO, ELEM_NUM, 1, 4980
*GET, NODE_I, ELEM, ELEM_NUM, NODE, 1
*GET, NODE_J, ELEM, ELEM_NUM, NODE, 2
*GET, XI, NODE, NODE_I, LOC, X
*GET, YI, NODE, NODE_I, LOC, Y
*GET, ZI, NODE, NODE_I, LOC, Z
*GET, XJ, NODE, NODE_J, LOC, X
*GET, YJ, NODE, NODE_J, LOC, Y
*GET, ZJ, NODE, NODE_J, LOC, Z
*SET, LENGTH, SQRT((XJ-XI)**2 + (YJ-YI)**2 + (ZJ-ZI)**2)/1e6 !单位
! 统一截面参数
*SET, OD ,0.45 ! 外径(m) !单位
*SET, ID , 0.42 ! 内径(m) !单位
*SET, PERIMETER , PI*OD ! 周长(m)
*SET, A_SURF , PERIMETER*LENGTH/2 !平均表面积分给节点
NODE_AREA(NODE_I) = NODE_AREA(NODE_I) + A_SURF
NODE_AREA(NODE_J) = NODE_AREA(NODE_J) + A_SURF
*ENDDO
*SET, RID34, 5000
*DO, NODE_NUM, 1, 200000
*IF, NODE_AREA(NODE_NUM), GT, 0, THEN
R, RID34, NODE_AREA(NODE_NUM), H2, 1000000 !为节点定义实常数,编号、表面积、综合换热系数、环境节点编号
REAL, RID34
TYPE, 5000
E, NODE_NUM, 1000000 ! 只连接一个结构节点和温度接待你
*SET, RID34, RID34+1
*ENDIF
*ENDDO
! 温度计算
*DIM, T5_ARRAY, ARRAY,4980
*DIM, NODE_TEMP_SUM, ARRAY,200000
*DIM, NODE_TEMP_COUNT, ARRAY,200000
*VOPER, NODE_TEMP_SUM(1), NODE_TEMP_SUM(1), ADD, 0
*VOPER, NODE_TEMP_COUNT(1), NODE_TEMP_COUNT(1), ADD, 0
*DIM,Q_ARRAY,ARRAY,4980
*SET,kkk,0.001 !单位,转为米
*DO, ELEM_NUM,1,4980
*GET, NODE_I, ELEM, ELEM_NUM, NODE, 1
*GET, NODE_J, ELEM, ELEM_NUM, NODE, 2
! 计算杆件方向
*GET, XI, NODE, NODE_I, LOC, X
*GET, YI, NODE, NODE_I, LOC, Y
*GET, ZI, NODE, NODE_I, LOC, Z
*GET, XJ, NODE, NODE_J, LOC, X
*GET, YJ, NODE, NODE_J, LOC, Y
*GET, ZJ, NODE, NODE_J, LOC, Z
VECTOR_X = (XJ - XI)*kkk !单位
VECTOR_Y = (YJ - YI)*kkk !单位
VECTOR_Z = (ZJ - ZI)*kkk !单位
VECTOR_MAG = SQRT(VECTOR_X**2 + VECTOR_Y**2)
*IF, VECTOR_MAG, LT, 1e-6, THEN
THETA_DEG = 0
*ELSE
DOT_PRODUCT = VECTOR_Y*(-1)
THETA_RAD = ACOS(DOT_PRODUCT/VECTOR_MAG)
THETA_DEG = THETA_RAD*180/PI
*IF, VECTOR_X, LT, 0, THEN
THETA_DEG = 180 - THETA_DEG
*ENDIF
*ENDIF
LENGTH = SQRT(VECTOR_X**2 + VECTOR_Y**2 + VECTOR_Z**2) !单位,前方向量转为米制了
*IF, LENGTH, LT, 1e-6, THEN
BETA_DEG = 0
*ELSE
BETA_RAD = ASIN(ABS(VECTOR_Z)/LENGTH)
BETA_DEG = BETA_RAD*180/PI
*ENDIF
! 计算太阳辐射
*SET, temp_beta_rad , BETA_DEG*(PI/180)
*SET, temp_A_theta_diff , ABS(A_deg - THETA_DEG)
*SET, temp_diff_rad , temp_A_theta_diff*(PI/180)
*SET, R , COS(temp_beta_rad)*sinh + SIN(temp_beta_rad)*cosh*COS(temp_diff_rad)
*IF, R, GT, 0, THEN
I2 = R*I1
*ELSE
I2 = 0
*ENDIF
*SET, S , (C*I2*(1 + COS(BETA_DEG*PI/180)))/2
*SET, L , 0.55 + 0.437*R + 0.313*R**2
*IF, R, GT, -0.2, THEN
*SET, G , L
*ELSE
*SET, G , 0.45
*ENDIF
*SET, S1 , C*I1*G
*SET, D1 , (I2 + S1)*D*(1 - COS(BETA_DEG*PI/180))/2
*SET, Q , (I2 + S1 + D1)*F !单位
*SET,Q_ARRAY(ELEM_NUM),Q
*SET, T5 , Q/H2 + T2
*SET, PERIMETER , PI*0.45 ! 周长(m)单位
*SET, HGEN , Q*PERIMETER ! 单位长度热生成率 (W/m)
BFE, ELEM_NUM, HGEN, , HGEN !通过BFE施加热流
*IF, NODE_I, GE, 1, AND, NODE_I, LE, 1000000, THEN
*SET, NODE_TEMP_SUM(NODE_I), NODE_TEMP_SUM(NODE_I) + T5
*SET, NODE_TEMP_COUNT(NODE_I), NODE_TEMP_COUNT(NODE_I) + 1
*ENDIF
*IF, NODE_J, GE, 1, AND, NODE_J, LE, 1000000, THEN
*SET, NODE_TEMP_SUM(NODE_J), NODE_TEMP_SUM(NODE_J) + T5
*SET, NODE_TEMP_COUNT(NODE_J), NODE_TEMP_COUNT(NODE_J) + 1
*ENDIF
*SET, T5_ARRAY(ELEM_NUM), T5
*ENDDO
*DO, NODE_NUM, 1, 200000
*IF, NODE_TEMP_COUNT(NODE_NUM), GT, 0, THEN
*SET, AVG_TEMP, NODE_TEMP_SUM(NODE_NUM)/NODE_TEMP_COUNT(NODE_NUM)
IC, NODE_NUM, TEMP, AVG_TEMP
*ENDIF
*ENDDO !提取节点温度作为初始温度
! 清理孤立节点
NSEL,S,LOC,Z,0,30000
NSEL,U,NODE,,1000000
ESLN,U,0
NDELE,ALL
ALLSEL
! 尝试选择目标节点
NSEL, S, NODE, , 1000000
*GET, node_count, NODE, 0, COUNT ! 获取当前选择集中的节点数量
*IF, node_count, EQ, 0, THEN
N, 1000000, 1e5, 1e5, 1e5
*ENDIF
D, 1000000, TEMP, T4
/INPUT,F:\\AREA_ID.txt
*SET, SAVED_delta, delta
*SET, SAVED_M, M
*SET, SAVED_B, B
*SET, SAVED_C, C
*SET, SAVED_kkk, kkk
FINISH
/SOLU
SOLCONTROL,ON
NROPT,FULL,,ON !完全牛顿法
!EQSLV,PCG !预条件共轭梯度
EQSLV,SPARSE !稀疏矩阵求解器
PIVCHECK,ON
CONV,HEAT,,0.05,2 !收敛准则放宽至0.1
PRED,ON
LNSRCH,ON
NEQIT,300
ANTYPE,4
TIME,3600
AUTOTS,ON
DELTIM,120,10,600
KBC,0
OUTRES,ALL,ALL
ESEL,S,LIVE
NSLE,S
NSEL,R,LOC,Z,0,30000
ALLSEL
SOLVE
FINISH
/POST1
/INPUT,F:\\kilometerss.txt
!/POST1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
/prep7
!RSYS,11
*DO, hour, 2, 24
/POST1
SET, LAST
*DIM, hour_temp, ARRAY, 2046
*VGET, hour_temp, NODE, 1, TEMP
/PREP7
! 恢复固定参数
delta = SAVED_delta
M = SAVED_M
B = SAVED_B
C = SAVED_C
kkk = SAVED_kkk
/INPUT,F:\\READ_DATA.txt
t_hour = 0.8 + hour
T = (t_hour - 12)*15
sinh_hour = SIN(26.04*PI/180)*SIN(delta*PI/180) + COS(26.04*PI/180)*COS(delta*PI/180)*COS(T*PI/180)
I1 = P*M / EXP(B/sinh_hour)
cosh_hour = SQRT(1 - sinh_hour**2)
sin_A_hour = (COS(delta*PI/180)*SIN(T*PI/180)) / cosh_hour
A_rad = ASIN(sin_A_hour)
A_deg = A_rad*180/PI
N, 1000000, 0, 30000, 20000
D, 1000000, TEMP, T2 + 273.15 ! 使用开尔文温度
*DIM,H2_ARR,ARRAY,2046 !综合换热系数因为边界条件的不一样,设数组
*DO, NODE_NUM, 1, 2046
*IF, R34_ARR(NODE_NUM), GT, 0, THEN ! 检查有效实常数ID
T1_hour = hour_temp(NODE_NUM)
T3_hour = T1_hour + 273.15
T4_hour = T2 + 273.15
H_hour = F1*z*(T3_hour**2 + T4_hour**2)*(T3_hour + T4_hour)
H1_hour = 2.6*(ABS(T1_hour-T2)**0.25 + 1.54*V)
H2_hour(NODE_NUM) = H_hour + H1_hour
! 使用存储的实常数ID修改
RMODIF, R34_ARR(NODE_NUM), 2, H2_hour
*ENDIF
*ENDDO
*DO, ELEM_NUM,1,4980
*GET, NODE_I, ELEM, ELEM_NUM, NODE, 1
*GET, NODE_J, ELEM, ELEM_NUM, NODE, 2
! 计算杆件方向
*GET, XI, NODE, NODE_I, LOC, X
*GET, YI, NODE, NODE_I, LOC, Y
*GET, ZI, NODE, NODE_I, LOC, Z
*GET, XJ, NODE, NODE_J, LOC, X
*GET, YJ, NODE, NODE_J, LOC, Y
*GET, ZJ, NODE, NODE_J, LOC, Z
VECTOR_X = (XJ - XI)*kkk !单位
VECTOR_Y = (YJ - YI)*kkk !单位
VECTOR_Z = (ZJ - ZI)*kkk !单位
VECTOR_MAG = SQRT(VECTOR_X**2 + VECTOR_Y**2)
*SET, PERIMETER , PI*OD
*SET, A_SURF , PERIMETER*LENGTH/2
NODE_AREA(NODE_I) = NODE_AREA(NODE_I) + A_SURF
NODE_AREA(NODE_J) = NODE_AREA(NODE_J) + A_SURF
*IF, VECTOR_MAG, LT, 1e-6, THEN
THETA_DEG = 0
*ELSE
DOT_PRODUCT = VECTOR_Y*(-1)
THETA_RAD = ACOS(DOT_PRODUCT/VECTOR_MAG)
THETA_DEG = THETA_RAD*180/PI
*IF, VECTOR_X, LT, 0, THEN
THETA_DEG = 180 - THETA_DEG
*ENDIF
*ENDIF
LENGTH = SQRT(VECTOR_X**2 + VECTOR_Y**2 + VECTOR_Z**2) !单位,前方向量转为米制了
*IF, LENGTH, LT, 1e-6, THEN
BETA_DEG = 0
*ELSE
BETA_RAD = ASIN(ABS(VECTOR_Z)/LENGTH)
BETA_DEG = BETA_RAD*180/PI
*ENDIF
! 计算太阳辐射
*SET, temp_beta_rad , BETA_DEG*(PI/180)
*SET, temp_A_theta_diff , ABS(A_deg - THETA_DEG)
*SET, temp_diff_rad , temp_A_theta_diff*(PI/180)
*SET, R_hour, COS(temp_beta_rad)*sinh_hour + SIN(temp_beta_rad)*cosh_hour*COS(temp_diff_rad)
*IF, R_hour, GT, 0, THEN
I2 = R_hour*I1
*ELSE
I2 = 0
*ENDIF
*SET, S_hour, (C*I2*(1 + COS(BETA_DEG*PI/180)))/2
*SET, L , 0.55 + 0.437*R_hour + 0.313*R_hour**2
*IF, R_hour, GT, -0.2, THEN
*SET, G_hour, L
*ELSE
*SET, G_hour, 0.45
*ENDIF
*SET, S1_hour , C*I1*G_hour
*SET, D1_hour , (I2 + S1_hour)*D*(1 - COS(BETA_DEG*PI/180))/2
*SET, Q_hour , (I2 + S1_hour+ D1_hour)*F !单位
*SET,Q_ARRAY(ELEM_NUM),Q_hour
*SET, H2_hour,H2_ARR(NODE_I) !单元左节点的H2
*SET, T5_hour , Q_hour/H2_hour + T2
*SET, PERIMETER , PI*0.45 ! 周长(m)单位
*SET, HGEN , Q_hour*PERIMETER ! 单位长度热生成率 (W/m)
BFE, ELEM_NUM, HGEN, , HGEN !通过BFE施加热流
*ENDDO上述命令流报错信息为 *** ERROR *** CP = 16.609 TIME= 09:53:59
Attempt to divide by zero in parameter expression.
*** ERROR *** CP = 16.609 TIME= 09:54:01
The above error occurred processing field= Q_HOUR/H2_HOUR+T2
Line= *SET, T5_hour, (Q_hour/H2_hour + T.
PARAMETER T5_HOUR = 0.000000000 查看H2_Hours为0的原因
最新发布