clear;clc;
%………………………………基本参数设定…………………工况以Li et al.UHPC板孔隙压力测试中PPST组分为例
h1=10; %沿板厚方向的单元尺寸，单位：mm
h2=1; %时间步长，单位：s
h3=120; %板厚,单位：mm
m=h3/h1; %沿板厚的单元数，个
h4=6; %受火时间，h
n=h4*3600/h2; %荷载步数
fcm=154.8; %HRPC立方体抗压强度，MPa
Pc=833; %每单位体积混凝土含水泥的质量,kg
hf=10; %迎火面热传导率
hc=9; %背火面热传导率
r4=0.12; %辐射率
To=25; %初始温度
z=0.72; %初始相对湿度
Po=101325; %大气压,Pa
M=0.018; %水的摩尔质量,kg
Pso=3169; %25℃时升温饱和水蒸气压力,Pa
Q0=0.0578; %HRPC的初始孔隙率
E1=0.00151; %HRPC的孔隙率系数A1
E2=2.65606*10^(-6); %HRPC的孔隙率系数A2
Vso=1-Q0; %固体体积分数
ktop=2.0*10^(-18); %初始渗透系数
Ad=4.6; %渗透性常数
a1=46.9364*10^6; %饱和度计算公式中：高性能混凝土常数a，pa
a2=2.0601; %饱和度计算公式中：高性能混凝土常数b
R=8.31; %气态方程参数
RV=461.5; %水蒸气气体参数，单位：J/（kg*K）
r1=h2/(h1*h1); % 网比，要求r<=1/2差分格式计算才能稳定
%………………………………液态水表面张力相关数据…………………
r0=-0.1965; %20℃时液态水表面张力
A1=3.55*10^(-17);
A2=9.07*10^(-14);
A3=9.61*10^(-11);
A4=5.35*10^(-8);
A5=1.62*10^(-5);
A6=2.66*10^(-3);
A7=2.71*10^(-3);
%…………………………………………………
%………………………………饱和水蒸汽压力相关数据…………………
D1=15.972;
D2=-13.09;
D3=0.029089;
D4=-5.5926*10^(-5);
D5=4.7229*10^(-8);
D6=1.7916;
%…………………………………………………
%………………………………赋存储空间,用于存储各节点不同时刻的值…………………
Tf=zeros(n+1); %火灾环境温度
T=zeros(n+1,m+1); %温度
Pv=zeros(n+1,m+1); %孔隙压力
Ps=zeros(n,m+1); %饱和水蒸气压力
Ps1=zeros(n,m+1); %饱和水蒸汽压力因温度增量
mv=zeros(n,m+1); %水蒸汽质量
Vv=zeros(n,m+1); %水蒸汽体积
Uv=zeros(n,m+1); %水蒸汽动力粘度
PL=zeros(n,m+1); %液态水密度
dPLdT=zeros(n,m+1); %dPL/dT
mL=zeros(n,m+1); %液态水质量
mD=zeros(n,m+1); %结合水转化为液态水的质量
dmLdPv=zeros(n,m+1); %dmL/dPv ，mL对孔隙压力求导
dmLdT=zeros(n,m+1); %dmL/dT，mL对温度求导
dmDdT=zeros(n,m+1); %dmD/dT，mD对温度求导
k1=zeros(n,m+1); %HRPC高温下的渗透系数
k3=zeros(n,m+1); %HRPC高温下的渗透系数因温度增量引起的变化
k4=zeros(n,m+1); %HRPC高温下的渗透系数因孔隙压增量引起的变化
Q=zeros(n,m+1); %HRPC高温下孔隙率变化
Q1=zeros(n,m+1); %HRPC高温下孔隙率因温度增量引起的变化
S=zeros(n,m+1); %饱和度
S1=zeros(n,m+1); %饱和度因孔隙压力增量引起的变化
S2=zeros(n,m+1); %饱和度因温度增量引起的变化
r=zeros(n,m+1); %高温下液态水表面张力
r5=zeros(n,m+1); %高温下液态水表面张力因温度增量引起的变化
A=zeros(n,m+1); %差分参数
B=zeros(n,m+1); %差分参数
C=zeros(n,m+1); %差分参数
e=zeros(n,m+1); %相对湿度
he=zeros(n,m+1); %辐射传热系数
kc=zeros(n,m+1); %导热系数
Cv=zeros(n,m+1); %比热容
w=zeros(n,m+1); %热扩散率
a=zeros(n,m+1); %差分项系数
b=zeros(n,m+1); %差分项系数
v=zeros(n,m+1); %差分项系数
Pv1=zeros(n); %某节点不同时刻的孔隙压力
Pv2=zeros(n); %某节点不同时刻的孔隙压力
%……………………………………温度场计算…………………………………………
for i=1:m+1
T(1,i)=To; %初始温度
end
for i=1:n+1 %火灾条件
Tf(i)=h2*i/30+To; %升温速度，2度/min
if Tf(i)>600; %温度达到600℃后温度恒定
Tf(i)=600;
end
end
for i=1:n
for j=2:m
kc(i,j)=1.44+1.85*exp(-1*T(i,j)/242.95); %导热系数
if 20<=T(i,j)&&T(i,j)<=100
Cv(i,j)=950;
else if 100<T(i,j)&&T(i,j)<=300
Cv(i,j)=950+(T(i,j)-100); %比热容
else if 300<T(i,j)&&T(i,j)<=600
Cv(i,j)=1150+(T(i,j)-300)/2;
else Cv(i,j)=1300;
end
end
end
v(i,m)=(hc*To)/kc(i,m); %差分项系数
b(i,m)=-hc/kc(i,m); %差分项系数
he(i,1)=r4*(((Tf(i)+273)^2+(T(i,1)+273)^2)*((Tf(i)+273)+(T(i,1)+273))*5.67*10^(-8));
a(i,1)=-(hf+he(i,1))/kc(i,2); %差分项系数
w(i,j)=1000000*kc(i,j)/(2500*Cv(i,j)); %热扩散率，其中2500为RPC的体积密度
T(i+1,1)=(1-2*w(i,2)*r1*(1-(a(i,1)*h1/1000)))*T(i,1)+2*w(i,2)*r1*T(i,2)+2*w(i,2)*r1*h1*((hf+he(i,1))*Tf(i))/(1000*kc(i,2)); %迎火面边界温度差分格式
T(i+1,m+1)=(1-2*w(i,m)*r1*(1-(b(i,m)*h1/1000)))*T(i,m+1)+2*w(i,m)*r1*T(i,m)+2*w(i,m)*r1*h1*v(i,m)/1000; %背火面边界温度差分格式
T(i+1,j)=T(i,j)+w(i,j)*r1*(T(i,j+1)-2*T(i,j)+T(i,j-1)); %截面内部温度差分格式
end
end
%……………………………………孔隙压力计算…………………………………………
for j=1:m+1;
Pv(1,j)=z*Pso; %初始孔隙压力
end
for i=1:n;
Pv(i+1,1)=z*Pso; %边界孔隙压力
Pv(i+1,m+1)=z*Pso; %边界孔隙压力
end
for i=1:n;
for j=2:m;
%……………………与液态水或者水蒸气有关的参数……………………
%PL(i,j)=(4.8863*10^(-7)-1.6528*10^(-9)*T(i,j)+1.8621*10^(-12)*T(i,j)^2+2.4266*10^(-13)*T(i,j)^3-1.5996*10^(-15)*T(i,j)^4+3.3703*10^(-18)*T(i,j)^5)*(-10*10^6)+(1.0213*10^3-7.7377*10^(-1)*T(i,j)+8.7696*10^(-3)*T(i,j)^2-9.2118*10^(-5)*T(i,j)^3+3.3534*10^(-7)*T(i,j)^4-4.4034*10^(-10)*T(i,j)^5); %液态水密度,考虑密度随温度的变化
%dPLdT(i,j)=(-1.6528*10^(-9)+2*1.8621*10^(-12)*T(i,j)+3*2.4266*10^(-13)*T(i,j)^2-4*1.5996*10^(-15)*T(i,j)^3+5*3.3703*10^(-18)*T(i,j)^4)*(-10*10^6)+(-7.7377*10^(-1)+2*8.7696*10^(-3)*T(i,j)-3*9.2118*10^(-5)*T(i,j)^2+4*3.3534*10^(-7)*T(i,j)^3-5*4.4034*10^(-10)*T(i,j)^4); %dPL/dT，对温度求导
PL(i,j)=1000; %液态水密度，忽略密度随温度的变化
dPLdT(i,j)=0; %PL对温度求导
Ps(i,j)=100*6.1121*exp(17.67*T(i,j)/(243.5+T(i,j))); %饱和水蒸汽压力,Pa
Ps1(i,j)=100*6.1121*exp(17.67*(1.0001*T(i,j))/(243.5+1.0001*T(i,j))); %饱和水蒸汽压力因温度增量引起的变化,Pa
Uv(i,j)=8.85*10^(-6)+3.53*10^(-8)*(T(i,j)-To); %水蒸汽动力粘度
r(i,j)=A1*(T(i,j)+273)^6-A2*(T(i,j)+273)^5+A3*(T(i,j)+273)^4-A4*(T(i,j)+273)^3+A5*(T(i,j)+273)^2-A6*(T(i,j)+273)+A7 ; %高温下液态水的表面张力
r5(i,j)=A1*(1.0001*T(i,j)+273)^6-A2*(1.0001*T(i,j)+273)^5+A3*(1.0001*T(i,j)+273)^4-A4*(1.0001*T(i,j)+273)^3+A5*(1.0001*T(i,j)+273)^2-A6*(1.0001*T(i,j)+273)+A7 ; %高温下液态水表面张力因温度增量引起的变化
%……………………与混凝土有关的参数……………………
if T(i,j)<=200 %随温度变化，结合水转化成液态水质量
mD(i,j)=0;
else if 200<T(i,j)&&T(i,j)<=300
mD(i,j)=0.07*Pc*((T(i,j)-200)/100); %Pc：每立方米混凝土中含水泥质量
else if 300<T(i,j)&&T(i,j)<=800
mD(i,j)=0.004*Pc*((T(i,j)-300)/100)+0.07*Pc;
else mD(i,j)=0.09*Pc;
end
end
end
if T(i,j)<=200 %……………………高温下结合水转化成液态水的质量对温度求导……………………
dmDdT(i,j)=0;
else if 200<T(i,j)&&T(i,j)<=300
dmDdT(i,j)=0.0007*Pc;
else if 300<T(i,j)&&T(i,j)<=800
dmDdT(i,j)=0.00004*Pc;
else dmDdT(i,j)=0;
end
end
end
Q(i,j)=Q0*(1+E1*(T(i,j)-To)+E2*(T(i,j)-To)^2); %高温下HRPC孔隙率
Q1(i,j)=Q0*(1+E1*(1.0001*T(i,j)-To)+E2*(1.0001*T(i,j)-To)^2); %高温下HRPC孔隙率因温度增量引起的变化
k1(i,j)=ktop*10^(Ad*(1+0.051-1.118*exp(-0.00311*T(i,j)))); %高温下HRPC渗透系数
k3(i,j)=ktop*10^(Ad*(1+0.051-1.118*exp(-0.00311*1.0001*T(i,j)))); %高温下HRPC渗透系数因温度增量引起的变化
k4(i,j)= k1(i,j); %高温下HRPC渗透系数因孔隙压增量引起的变化
%……………………高温下混凝土饱和度S的计算公式……………………
if 0<Pv(i,j)/Ps(i,j)<=1
S(i,j)=((-1/a1*RV*(T(i,j)+273)*PL(i,j)*log(Pv(i,j)/Ps(i,j))*r0/r(i,j)*(Q0/ktop)^(1/2)*(k1(i,j)./Q(i,j))^(1/2))^(1/(1-1/a2))+1)^(-1/a2);
else S(i,j)=((1/a1*RV*(T(i,j)+273)*PL(i,j)*log(1)*r0/r(i,j)*(Q0/ktop)^(1/2)*(k1(i,j)/Q(i,j))^(1/2))^(1/(1-1/a2))+1)^(-1/a2);
end
if 0<1.0001*Pv(i,j)/Ps(i,j)<=1
S1(i,j)=((-1/a1*RV*(T(i,j)+273)*PL(i,j)*log((1.0001*Pv(i,j))/Ps(i,j))*r0/r(i,j)*(Q0/ktop)^(1/2)*(k4(i,j)/Q(i,j))^(1/2))^(1/(1-1/a2))+1)^(-1/a2);
else S1(i,j)=((1/a1*RV*(T(i,j)+273)*PL(i,j)*log(1)*r0/r(i,j)*(Q0/ktop)^(1/2)*(k4(i,j)/Q(i,j))^(1/2))^(1/(1-1/a2))+1)^(-1/a2);
end
if 0<Pv(i,j)/Ps1(i,j)<=1
S2(i,j)=((-1/a1*RV*(1.0001*T(i,j)+273)*PL(i,j)*log(Pv(i,j)/Ps1(i,j))*r0/r5(i,j)*(Q0/ktop)^(1/2)*(k3(i,j)/Q1(i,j))^(1/2))^(1/(1-1/a2))+1)^(-1/a2);
else S2(i,j)=((1/a1*RV*(1.0001*T(i,j)+273)*PL(i,j)*log(1)*r0/r5(i,j)*(Q0/ktop)^(1/2)*(k3(i,j)/Q1(i,j))^(1/2))^(1/(1-1/a2))+1)^(-1/a2);
end
%…………………………………水蒸汽体积和水蒸汽质量……………………………………
Vv(i,j)=1-Q(i,j)*S(i,j)-(Vso-mD(i,j)/PL(i,j)); %水蒸汽体积
mv(i,j)=Pv(i,j)*M*Vv(i,j)/(R*(T(i,j)+273)); %水蒸汽质量
%……………………………………液态水质量表达式………………………………………
mL(i,j)=PL(i,j)*Q(i,j)*S(i,j);
%……………………………………液态水质量对孔隙压力求导……………………………
dmLdPv(i,j)=PL(i,j)*Q(i,j)*(S1(i,j)-S(i,j))/(0.0001*Pv(i,j));
%……………………………………液态水质量对温度求导…………………………………
dmLdT(i,j)=(PL(i,j)*Q1(i,j)*S2(i,j)-PL(i,j)*Q(i,j)*S(i,j))/(0.0001*T(i,j));
%……………………………………相关参数表达式…………………………………………
A(i,j)=(1-mv(i,j)/(Vv(i,j)*PL(i,j)))*dmLdPv(i,j)+Vv(i,j)*M/(R*(T(i,j)+273));
B(i,j)=mv(i,j)*k1(i,j)/Uv(i,j);
C(i,j)=((1-mv(i,j)/(Vv(i,j)*PL(i,j)))*(-dmLdT(i,j)+dmDdT(i,j))+mv(i,j)/(T(i,j)+273)+mv(i,j)/(Vv(i,j)*PL(i,j)^2)*dPLdT(i,j)*(mD(i,j)-mL(i,j)))*(T(i+1,j)-T(i,j))/h2;
%………………………………内部孔隙压力差分格式………………………………………
Pv(i+1,j)=Pv(i,j)+(1000000*r1*B(i,j)*(Pv(i,j+1)-2*Pv(i,j)+Pv(i,j-1)))/A(i,j)+C(i,j)*h2/A(i,j); %内部孔隙压差分格式
if Pv(i+1,j)<0 %消除复数影响
Pv(i+1,j)=z*Pso;
end
end
end
%……………………………………数据提取…………………………………………
for i=1:n
% Pv1(i)=Pv(i,2)/1000000; %P10测点孔隙压力，MPa
Pv2(i)=Pv(i,3)/1000000; %P20测点孔隙压力，MPa
end
这段Matlab代码主要实现了对超高性能混凝土（UHPC）板在火灾条件下的温度场和孔隙压力场的计算与分析。具体功能如下：
- 参数设定：定义了UHPC板的几何尺寸、材料属性、热物理参数、时间步长、荷载步数等基本参数。
- 存储空间分配：为存储各节点在不同时刻的温度、孔隙压力、饱和水蒸气压力、水蒸汽质量、液态水密度等参数分配了相应的矩阵空间。
- 温度场计算：
- 根据给定的升温速度和温度上限，计算火灾环境温度随时间的变化。
- 初始化UHPC板的初始温度分布。
- 基于热传导理论和差分格式，计算UHPC板在不同时刻沿板厚方向的温度分布，包括迎火面、背火面和截面内部的温度。
- 孔隙压力计算：
- 初始化UHPC板的初始孔隙压力分布。
- 根据温度场计算结果，计算与液态水或水蒸气有关的参数，如液态水密度、饱和水蒸气压力、水蒸汽动力粘度、液态水表面张力等。
- 计算与混凝土有关的参数，如结合水转化为液态水的质量、高温下HRPC孔隙率、渗透系数、饱和度等。
- 计算水蒸汽体积和质量、液态水质量及其对孔隙压力和温度的导数。
- 基于上述参数，采用差分格式计算UHPC板在不同时刻沿板厚方向的孔隙压力分布，并处理孔隙压力为负的情况。
- 数据提取：提取特定节点（如P20测点）在不同时刻的孔隙压力数据，单位转换为MPa。
通过上述计算，该代码能够模拟UHPC板在火灾条件下的温度和孔隙压力变化，为进一步研究UHPC板的热-水-力耦合行为提供了基础数据。
扫一扫
401
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
