> restart;

> with(DEtools):

> with(linalg):

Warning, new definition for adjoint

Warning, new definition for norm

Warning, new definition for trace

> StartTemp:=310:

> RorDiameter:=1.02:

> Massestrom:=1000:

> StartTrykk:=20000000:

> RorRuhet:=0.00001:

> DeltaTrykk:=50000:

> DeltaTemp:=1:

> Cp:=2000:

> VannTemp:=275:

> KonduktivitetBetong:=2.8:

> KonduktivitetGass:=0.025:

> KonduktivitetStaal:=50:

> YtreVarmeovergangstall:=25:

> YtreDiameterBetong:=1.350:

> YtreDiameterStaal:=1.050:

>

>

> KomponentVektor:=array(1..5):
KomponentVektor[1]:=Metan: KomponentVektor[2]:=Etan: KomponentVektor[3]:=Propan: KomponentVektor[4]:=Butan: KomponentVektor[5]:=Pentan:

> MolfraksjonVektor:=array(1..5):
MolfraksjonVektor[1]:=0.85: MolfraksjonVektor[2]:=0.05: MolfraksjonVektor[3]:=0.02:
MolfraksjonVektor[4]:=0.06: MolfraksjonVektor[5]:=0.02:

> MolmasseVektor:=array(1..5):
MolmasseVektor[1]:=16.04: MolmasseVektor[2]:=30.08: MolmasseVektor[3]:=46.76:
MolmasseVektor[4]:=60.02: MolmasseVektor[5]:=76.3:

> TCVektor:=array(1..5):
TCVektor[1]:=190.555: TCVektor[2]:=305.45: TCVektor[3]:=369.82:
TCVektor[4]:=425.15: TCVektor[5]:=469.65:

> PCVektor:=array(1..5):
PCVektor[1]:=4600100: PCVektor[2]:=4883900: PCVektor[3]:=4249600:
PCVektor[4]:=3850046: PCVektor[5]:=3413741:

> Molmasse:=multiply(MolfraksjonVektor,MolmasseVektor):

> TC:=multiply(MolfraksjonVektor,TCVektor):

> PC:=multiply(MolfraksjonVektor,PCVektor):

> TR:=Temp->Temp/TC:

> PR:=Pres->Pres/PC:

> GasGravity:=Molmasse/28.97:

> g:=9.81:

> R:=8314.41:

> T0=298:

> P0=101325:

> fZ:=1:

> Z:=proc(Ppr, Tpr)
option remember;
local T2, T3, T4, T5, rho, Zl, F, FP,
A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,Z1;
global fZ;
if not(type([args],list(numeric)))
then 'procname(args)';
elif nargs>0 then
A1 := 0.3265: A2 := -1.07: A3 := -0.5339: A4 := 0.01569:
A5 := -0.05165: A6 := 0.5475: A7 := -0.7361: A8 := 0.1844:
A9 := 0.1056: A10 := 0.6134: A11 := 0.721;
Z1 := 2;
if (Ppr > 0) then
T2 := A1 * Tpr + A2 + A3 / (Tpr ^ 2) + A4 / (Tpr ^ 3)
+ A5 / (Tpr ^ 4);
T3 := A6 * Tpr + A7 + A8 / Tpr;
T4 := -A9 * (A7 + A8 / Tpr);
T5 := A10 / (Tpr ^ 2);
rho := 0.27 * Ppr / Tpr;
while (abs(fZ - Z1) > 0.005) do
Z1 := fZ;
F := rho * (Tpr + T2 * rho + T3 * rho ^ 2 + T4 * rho ^ 5
+ T5 * rho ^ 2 * (1 + A11 * rho ^ 2)
* exp(-A11 * rho ^ 2)) - 0.27 * Ppr;
FP := Tpr + 2 * T2 * rho + 3 * T3 * rho ^ 2
+ 6 * T4 * rho ^ 5 + T5 * rho ^ 2
* exp(-A11 * rho ^ 2) * (3 + A11 * rho ^ 2
* (3 - 2 * A11 * rho ^ 2));
rho := rho - F / FP;
fZ := 0.27 * Ppr / Tpr / rho;
od;
fi;
RETURN(fZ);
fi;
end:

> Viskositet:=proc(Temp,Pres,Gg,TC,PC)
option remember;
local PZ, Tf, TR, PR, Visk;
if not(type([args], list(numeric)))
then 'procname(args)';
elif nargs>0 then
Tf:=Temp*1.8-241;
TR:=Temp/TC;
PR:=Pres/PC;
PZ:=Z(PR,TR);
Visk:=Viskos(PZ,Gg,Tf);
RETURN(Visk);
fi;
end:

> Viskositet(Temp,Pres,GasGravity,TC,PC):

> Viskos:=proc(PZ, Gg, Tf)
local Mg,Tr,fMug;
Mg := Gg * 28.9625;
Tr := Tf + 459.67;
fMug := (7.77 + 0.0063 * Mg) * Tr ^ 1.5
/ (122.4 + 12.9 * Mg + Tr) / 10000
* ViskRatio(PZ, Gg, Tf, Tr, Mg);
RETURN (fMug*0.001);
end:

> ViskRatio:=proc(PZ, Gg, Tf, Tr, Mg)
local SGg, B, C, fMugr;
SGg := fRhog(PZ, Gg, Tf) / 62.368;
B := 2.57 + 1914.5 / Tr + 0.0095 * Mg;
C := 1.11 + 0.04 * B;
fMugr := exp(B * SGg ^ C);
RETURN (fMugr);
end:

> fRhog:=proc(PZ, Gg, Tf)
local Mg, Tr,fRhog;
Mg := Gg * 28.9625;
Tr := Tf + 459.67;
fRhog := PZ * Mg / 10.7315 / Tr;
RETURN (fRhog);
end:

> Tetthet:=(Trykk,Temp,Molmasse)->(Trykk*Molmasse)/(Z(Trykk/PC,Temp/TC)*R*Temp):

> Areal:=(Diameter)->3.14*Diameter^2/4:

> Hastighet:=(Trykk,Temp,Molmasse,Diameter,Massestrom)->Massestrom/(Tetthet(Trykk,Temp,Molmasse)/Areal(Diameter)):

> Reynoldstall:=(Trykk,Temp,Molmasse, Massestrom, Diameter,TC,PC)->Tetthet(Trykk,Temp,Molmasse)*Hastighet(Trykk,Temp,Molmasse,Diameter,Massestrom)* Diameter/Viskositet(Temp,Trykk,GasGravity,TC,PC):

> friksjonsfaktor:=(Trykk,Temp,Molmasse,Diameter, epsilon,Massestrom)->simplify(1/(2*log10(3.7/(epsilon/Diameter))))^2:

> friksjonsfaktor2:=(Trykk,Temp,Diameter, epsilon)->simplify(1/(-1.8*log10(6.9/Reynoldstall(Trykk,Temp,Molmasse, Massestrom, Diameter,TC,PC)+(epsilon/Diameter/3.7)^1.11)))^2:

> DiffHastighet:=(Trykk,Temp,Molmasse,Diameter,Massestrom)->(Hastighet(Trykk-2*DeltaTrykk,Temp,Molmasse,Diameter,Massestrom)-8*Hastighet(Trykk-DeltaTrykk,Temp,Molmasse,Diameter,Massestrom)+8*Hastighet(Trykk+DeltaTrykk,Temp,Molmasse,Diameter,Massestrom)-Hastighet(Trykk+2*DeltaTrykk,Temp,Molmasse,Diameter,Massestrom))/12/DeltaTrykk:

> Prandtltall:=(Temp, Trykk)->Cp*Viskositet(Temp,Trykk,GasGravity,TC,PC)/KonduktivitetGass:

IndreVarmeovergangstall:=(Temp, Trykk)->0.023*Reynoldstall(Trykk,Temp,Molmasse, Massestrom, RorDiameter,TC,PC)^0.8*Prandtltall(Temp, Trykk)^0.33*KonduktivitetGass/RorDiameter:

TotaltVarmeovergangstall:=(Temp, Trykk)->simplify(1/(1/IndreVarmeovergangstall(Temp, Trykk)+(RorDiameter/(2*KonduktivitetBetong)*log(YtreDiameterBetong/YtreDiameterStaal))+(RorDiameter/(2*KonduktivitetStaal)*log(YtreDiameterStaal/RorDiameter))+RorDiameter/(YtreDiameterBetong*YtreVarmeovergangstall))):

> JouleThomson:=(Temp,Trykk)->(((Z(Trykk/PC,(Temp-2*DeltaTemp)/TC))-8*(Z(Trykk/PC,(Temp-DeltaTemp)/TC))+8*(Z(Trykk/PC,(Temp+DeltaTemp)/TC))-(Z(Trykk/PC,(Temp+2*DeltaTemp)/TC)))/12/DeltaTemp)*(Temp^2)*(R/Molmasse)/(Cp*Trykk):

Setter opp bevaringsloven for impuls og energi:

> Impulsligning:=(Trykk,Temp,Molmasse,Diameter,Massestrom,epsilon)->Hastighet(Trykk,Temp,Molmasse,Diameter,Massestrom)*DiffHastighet(Trykk,Temp,Molmasse,Diameter,Massestrom)*diff(Trykk(x),x)=-1/Tetthet(Trykk,Temp,Molmasse)*diff(Trykk(x),x)-friksjonsfaktor(Trykk,Temp,Molmasse,Diameter, epsilon,Massestrom)*(Hastighet(Trykk,Temp,Molmasse,Diameter,Massestrom))^2/Diameter/2:

Energiligning:=(Trykk,Temp,Molmasse,Diameter,Massestrom,epsilon,JouleThomson,Cp,Varmeovergangstall,VannTemp)->Cp*(diff(Temp(x),x)-JouleThomson(Temp,Trykk)*diff(Trykk(x),x))=(4*TotaltVarmeovergangstall(Temp, Trykk)*(VannTemp-Temp))/(Tetthet(Trykk,Temp,Molmasse)*Diameter*Hastighet(Trykk,Temp,Molmasse,Diameter,Massestrom)):

> Digits:=5:

### WARNING: `dsolve` has been extensively rewritten, many new result forms can occur and options are slightly different, see help page for details
t:=dsolve({Impulsligning(Trykk,Temp,Molmasse,RorDiameter,Massestrom,RorRuhet),Energiligning(Trykk,Temp,Molmasse,RorDiameter,Massestrom,Rorruhet,JouleThomson,Cp,Varmeovergangstall,VannTemp),Trykk(0)=StartTrykk, Temp(0)=StartTemp},{Trykk(x),Temp(x)}, type=numeric);

> t(0);
t(1);
t(50);
t(100000);
t(600000);
t(728000);

> with(plots):
plotsetup(default):
odeplot(t,[x,Trykk(x)],0..728000,title=Trykkfall);
odeplot(t, [x,Temp(x)],0..728000, title=Temperatur);