matlab_2_matlab_粘度_温度_瞬变流_

% 加入了能量方程，对Cd也进行了插值 % 用上一时刻的数据来计算密度、粘度 % 采用新的温度变化曲线，水深95米，管道长度计算长度为1300m。 % 温度、压力对密度、粘度的影响 % 初始数据： clear; tic N=1000; % 分段数，运行6057.110521 秒。 Hr=30; % 水库高度，MPa NS=N+1; L1=95; L2=1206; L=L1+L2; % 管线长度L dx=L/N; % 分段长度 e=0.001651; % 壁厚,0.065'' D=0.00635-2*e; % 管道内径 0.00635mm Area=pi*D^2/4; % 横截面积 G=9.806; % 重力加速度 Tem=zeros(1,NS); % 温度矩阵 eta=zeros(1,NS); % 动力/绝对粘度矩阵 Rho=zeros(1,NS); % 密度矩阵 a=zeros(1,NS); % 波速矩阵 P=zeros(1,NS); % 管道内部压力矩阵 P(1,1)=0.1; % 管道内部初始压力 % ---------------能量方程中的参数-------------------- a1=13.705314387398905019151955113057;% 以a1表示，与波速a区分 b=0.0000000000027166961307500002704551085603983; c=-0.000012203798517995766309649779997993; d=-27326.165179604799509539992131829; e=0.012166193517894900889946923805305; % -------------------------------------------------- for j=1:N %---------根据温度变化计算一个初始值------------- %---------------包括 密度、粘度------------------ if j<=(L1/L)*N Tem(1,1)=30; Tem(1,j+1)=Tem(1,1)-j*(30-16)/(N*(L1/L)); elseif j>=(L1/L)*N Tem(1,j+1)=16+(j-L1/dx)*(66/(N*(L2/L))); end % 计算每个节点的温度 Rho(1,1)=md(Tem(1,1),P(1,1)); % 计算第一点的密度 P(1,j+1)=0.1; Rho(1,j+1)=md(Tem(1,j+1),P(1,j+1)); % 计算密度 eta(1,1)=eta_d(Tem(1,1),P(1,1)); % 计算第一点的动力/绝对粘度 eta(1,j+1)=eta_d(Tem(1,j+1),P(1,1+1)); % 计算动力/绝对粘度 a(1,1)=sqrt((P(1,1)/Rho(1,1)^2-2*a1*Rho(1,1)-c*P(1,1)-d)/(2*b*P(1,1)+c*Rho(1,1)+e));% 计算第一点的波速 a(1,j+1)=sqrt((P(1,j+1)/Rho(1,j)^2-2*a1*Rho(1,j+1)-c*P(1,j+1)-d)/(2*b*P(1,j+1)+c*Rho(1,j+1)+e));% 计算波速 end dt=dx/1502.8901734104; % 时间步长 t_max=2400; % 总计算时间 k=1; % 节点计数 t=0; % 时间 w=0; % 用于计数，w每增加1表示t增加一个dt int=1000; % 表示每过1000个dt,记录 一行数据 time=linspace(1,t_max,round(t_max/int/dt)); % 时间矩阵 V=zeros(2,NS); % 流量矩阵 H=zeros(2,NS); % 压力矩阵 H_end=zeros(round(t_max/int/dt),N/5); % 末端压力矩阵 V(1,:)=0; % 第一行流量 H(1,:)=0; % 第一行压力 H(1,1)=Hr; % H(2,1)=Hr; while t<t_max if t<=t_max/2 % 对时间进行判断 H(2,1)=Hr; else H(2,1)=0; end for k=2:N % ------------------------------------------------------------------------------------------------- Rho(1,1)=md(Tem(1,1),P(1,1)); P(1,1)=Hr; P(1,k)=H(1,k); Rho(1,k)=md(Tem(1,k),P(1,k)); % 计算密度 eta(1,1)=eta_d(Tem(1,1),P(1,1)); eta(1,k)=eta_d(Tem(1,k),P(1,k)); % 计算动力/绝对粘度 a(1,1)=sqrt((H(1,1)/Rho(1,1)^2-2*a1*Rho(1,1)-c*H(1,1)-d)/(2*b*H(1,1)+c*Rho(1,1)+e)); a(1,k)=sqrt((H(1,k)/Rho(1,k)^2-2*a1*Rho(1,k)-c*H(1,k)-d)/(2*b*H(1,k)+c*Rho(1,k)+e)); %--------------------------------------------------------------------------------------------------- %--------------------------------------------------------------------------------------------------- aR=a(1,k)/(1+(a(1,k)-a(1,k-1))*(dt/dx)); % 插值法计算波速aR aS=a(1,k)/(1+(a(1,k)-a(1,k+1))*(dt/dx)); % 插值法计算波速aS VR=(V(1,k)-(V(1,k)-V(1,k-1))*(dt/dx)*aR)/(1+(dt/dx)*(V(1,k)-V(1,k-1))); % 插值法计算流速VR VS=(V(1,k)-(V(1,k)-V(1,k+1))*(dt/dx)*aS)/(1-(dt/dx)*(V(1,k)-V(1,k+1))); % 插值法计算流速VS HR=H(1,k)-(VR*(dt/dx)+(dt/dx)*aR)*(H(1,k)-H(1,k-1)); % 插值法计算压力HR HS=H(1,k)+(VS*(dt/dx)-(dt/dx)*aS)*(H(1,k)-H(1,k+1)); % 插值法计算压力HS RhoR=Rho(1,k)-(VR+aR)*(dt/dx)*(Rho(1,k)-Rho(1,k-1)); % 插值法计算密度RhoR RhoS=Rho(1,k)+(VS-aS)*(dt/dx)*(Rho(1,k)-Rho(1,k+1)); % 插值法计算密度RhoS etaR=eta(1,k)-(VR+aR)*(dt/dx)*(eta(1,k)-eta(1,k-1)); % 插值法计算粘度etaR etaS=eta(1,k)+(VS-aS)*(dt/dx)*(eta(1,k)-eta(1,k+1)); % 插值法计算粘度etaS % -------------------------------------------------------------------------------------------------- % -------根据含有能量方程的特征线法计算V↓------- Xr=HR/RhoR^2-2*a1*RhoR-c*HR-d; Yr=2*b*HR+c*RhoR+e; Xs=HS/RhoS^2-2*a1*RhoS-c*HS-d; Ys=2*b*HS+c*RhoS+e; A1=HR-HS+(RhoR/Yr)*sqrt(Xr*Yr)*VR+(RhoS/Ys)*sqrt(Xs*Ys)*VS; f=f_nu(etaR,RhoR,VR); B1=(1/Yr)*((-RhoR*sqrt(Xr*Yr)*f*VR*abs(VR)/(2*D))+(RhoR/etaR)*(f*VR*abs(VR)/8)^2+G*VR)*dt; f=f_nu(etaS,RhoS,VS); C1=(1/Ys)*(( RhoS*sqrt(Xs*Ys)*f*VS*abs(VS)/(2*D))+(RhoS/etaS)*(f*VS*abs(VS)/8)^2+G*VS)*dt; V(2,k)=(1/(((RhoR/Yr)*sqrt(Xr*Yr))+((RhoS/Ys)*sqrt(Xs*Ys))))*(A1+B1-C1); % -------------------------------------------------------------------------------------------------- % -------根据含有能量方程的特征线法计算H↓------- A2=HR+HS-(RhoR/Yr)*sqrt(Xr*Yr)*(V(2,k)-VR)+(RhoS/Ys)*sqrt(Xs*Ys)*(V(2,k)-VS); f=f_nu(etaR,RhoR,VR); B2=(1/Yr)*((-RhoR*sqrt(Xr*Yr)*f*VR*abs(VR)/(2*D))+(RhoR/etaR)*(f*VR*abs(VR)/8)^2+G*VR)*dt; f=f_nu(etaS,RhoS,VS); C2=(1/Ys)*(( RhoS*sqrt(Xs*Ys)*f*VS*abs(VS)/(2*D))+(RhoS/etaS)*(f*VS*abs(VS)/8)^2+G*VS)*dt; H(2,k)=0.5*(A2+B2+C2); % ------------------------------------------------------------------------------------------------- % -------根据含有能量方程的特征线法计算Rho↓------- VD=V(1,k)/(1+(V(1,k)-V(1,k-1))*(dt/dx)); HD=H(1,k)-(VD*(dt/dx))*(H(1,k)-H(1,k-1)); RhoD=Rho(1,k)-(VD*(dt/dx))*(Rho(1,k)-Rho(1,k-1)); etaD=eta(1,k)-(VD*(dt/dx))*(eta(1,k)-eta(1,k-1)); Xd=-2*a1*RhoD-c*HD-d+HD/RhoD^2 ; Yd=2*b*HD+c*RhoD+e; f=f_nu(etaD,RhoD,VD); Rho(2,k)=RhoD+(Yd/Xd)*(H(2,k)-HD)-(1/Xd)*((RhoD/etaD)*(f*VD*abs(VD)/8)^2+G*VD)*dt; end % -------边界条件↓--------------------------------------------------------------------------------- % -------边界条件↓--------------------------------------------------------------------------------- % -------边界条件↓--------------------------------------------------------------------------------- P(1,NS)=H(1,NS); Rho(1,NS)=md(Tem(1,NS),P(1,NS)); % 计算边界密度 eta(1,NS)=eta_d(Tem(1,NS),P(1,NS)); % 计算边界动力/绝对粘度 a(1,NS)=sqrt((H(1,NS)/Rho(1,NS)^2-2*a1*Rho(1,NS)-c*H(1,NS)-d)/(2*b*H(1,NS)+c*Rho(1,NS)+e));% 计算边界波速 % 对于H（2，1），等于Hr % H(2,1)=Hr; % ------------------------------------------------------------------------------------------------- % 对于V（2，1），利用C-以及（1，2）计算 aS=a(1,1)/(1+(a(1,1)-a(1,2))*(dt/dx)); % 插值法计算波速aS VS=(V(1,1)-(V(1,1)-V(1,2))*(dt/dx)*aS)/(1-(dt/dx)*(V(1,1)-V(1,2))); % 插值法计算流速VS HS=H(1,1)+(VS*(dt/dx)-(dt/dx)*aS)*(H(1,1)-H(1,2)); % 插值法计算压力HS RhoS=Rho(1,1)+(VS-aS)*(dt/dx)*(Rho(1,1)-Rho(1,2)); % 插值法计算密度RhoS etaS=eta(1,1)+(VS-aS)*(dt/dx)*(eta(1,2)-eta(1,2)); % 插值法计算粘度etaS Xs=HS/RhoS^2-2*a1*RhoS-c*HS-d; Ys=2*b*HS+c*RhoS+e; f=f_nu(etaS,RhoS,VS); C4=((RhoS*sqrt(Xs*Ys)*f*VS*abs(VS)/(2*D))+(RhoS/etaS)*(f*VS*abs(VS)/8)^2+G*VS)*dt; V(2,1)=VS+(1/(RhoS*sqrt(Xs*Ys)))*(Ys*(H(2,1)-HS)-C4); % -----------------