Embedded Files
matlab contains cross terms
excel (no cross terms), uses linest function
/***************************************************************************************************
/* the following UDF source code contains three functions for tep2190oil density and viscosity
/* as a function of temperature and pressure
/* March 11th, 2020
1) a custom speed of sound function named sound_speed that is defined in terms of temperature and pressure
2) a custom density function named cell_density that is defined in terms of temperature and pressure
3) a custom viscosity function named cell_viscosity that is defined in terms of temperature and pressure
Michael Thompson
**************************************************************************************************/
#include "udf.h"
DEFINE_PROPERTY(sound_speed, cell,thread)
{
real density;
real temp = C_T(cell, thread);
real press= C_P(cell, thread);
real a;
real d_den_dp;
real BMODULUS;
density= -0.003549783549783796*temp*temp+9.086811117268821e-10*temp*press+ 1.725469755469909*temp+-2.190068559062531e-15;
d_den_dp=9.086811117268821e-10*temp+ -2.190068559062531e-15*2*press+ 2.86986682925744e-07;
BMODULUS= density/d_den_dp;
a = sqrt(BMODULUS/density);
return a;
}
DEFINE_PROPERTY(cell_density,cell,thread)
{real density;
real temp=C_T(cell,thread);
real press=C_P(cell,thread);
density= -0.00380952*temp*temp+1.86209524*temp+643.4579142+(press/6895000)*(0.00009524*temp*temp-0.05155238*temp+9.70430214);
return density;}
DEFINE_PROPERTY(cell_viscosity,cell, thread)
{real mu_lam;
real temp=C_T(cell,thread);
real press=C_P(cell,thread);
mu_lam= 0.00012520571*temp*temp-0.083207*temp+13.846+(press/6895000)*(0.00001895*temp*temp-0.01262935*temp+2.10840086);
return mu_lam;}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main Matlab File for:
% 1. Read in Reference NAVSEA 2190 measurement data
% 2. Create Sargent 2190 measurment data
% 3. Create analytical equations density and viscosity for sargent 2190
% March 06 2020
% Paul Kalbfleisch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% addtional Required Files:
% NAVSEA.xlsx, 2190_dec_data.xlsx
% Reset Matlab Workspace
close all; clear; clc;
local = pwd;
font_size = 18;
font_name = 'times';
linewidth = 0.8;
scrsz = get(0,'screensize');
fig_width = 800;
fig_height = 600;
% INPUTS
navsea.loc = [local '\NAVSEA.xlsx'];
sargent.loc = [local '\2190_dec_data.xlsx'];
sargent.eqs = [local '\Eqs_dec_data.xlsx'];
sargent.txt = [local '\Eqs_dec_data.txt'];
sargent.tables = [local '\sargent_output_tables.xlsx'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% clean data / Make structs
navsea.kin = xlsread(navsea.loc,'kin_visc');
navsea.den = xlsread(navsea.loc,'density');
navsea.dyn = xlsread(navsea.loc,'dyn_visc');
sargent.data = xlsread(sargent.loc);
sargent.F = sargent.data(:,1);
sargent.cSt = sargent.data(:,2);
%% Convert to SI
navsea.K = navsea.kin(1,2:end)+273; % c to K
navsea.F = (navsea.K - 273.15)* 1.8 + 32; %K to F
navsea.psig = navsea.kin(2:end,1);
navsea.pa = navsea.psig*6895; %psig to pa
navsea.kinSI = navsea.kin(2:end,2:end) * 1E-6; %Cst to m^2/s
navsea.denSI = navsea.den(2:end,2:end);
navsea.dynSI = navsea.dyn(2:end,2:end) / 1000; %cp to kg/m-s
sargent.K = (sargent.F + 459.67) * 5/9;
sargent.kinSI = sargent.cSt * 1E-6; %cSt to m^2/s
%% Create Meshgrid
[navsea.xg,navsea.yg] = meshgrid(273:5:100+273,0:1E5:10000*6895);
sargent.xg = navsea.xg; sargent.yg = navsea.yg;
%% Fit Navsea Density Data
count = 1;
for i = 1:length(navsea.K)
for ii = 1:length(navsea.pa)
navsea.den_fit(count,1) = navsea.K(i);
navsea.den_fit(count,2) = navsea.pa(ii);
navsea.den_fit(count,3) = navsea.denSI(ii,i);
count = count + 1;
end
end
% Polyfit
navsea.P_den = polyfitn([navsea.den_fit(:,1) navsea.den_fit(:,2)],navsea.den_fit(:,3),2);
navsea.z_den = polyvaln(navsea.P_den,[navsea.xg(:),navsea.yg(:)]);
% Plot Polyfit
fig_name = ['Navsea Density'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(navsea.den_fit(:,1),navsea.den_fit(:,2),navsea.den_fit(:,3));
hold on;
title(fig_name)
surf(navsea.xg,navsea.yg,reshape(navsea.z_den,size(navsea.xg)),'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Density [$\frac{kg}{m^3}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Fit Navsea Dynamic Viscosity Data
count = 1;
for i = 1:length(navsea.K)
for ii = 1:length(navsea.pa)
navsea.dyn_fit(count,1) = log(navsea.K(i));
navsea.dyn_fit(count,2) = navsea.pa(ii);
navsea.dyn_fit(count,3) = log(navsea.dynSI(ii,i));
count = count + 1;
end
end
navsea.P_dyn = polyfitn([navsea.dyn_fit(:,1) navsea.dyn_fit(:,2)],navsea.dyn_fit(:,3),2);
navsea.z_dyn = polyvaln(navsea.P_dyn,[log(navsea.xg(:)),navsea.yg(:)]);
fig_name = ['NavSea Dynamic Viscosity'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(navsea.dyn_fit(:,1),navsea.dyn_fit(:,2),navsea.dyn_fit(:,3));
title(fig_name)
hold on;
surf(log(navsea.xg),navsea.yg,reshape(navsea.z_dyn,size(navsea.xg)),'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('ln(Temp) [ln(K)]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('ln(Viscosity) [$ln(\frac{kg}{ms})$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Fit Navsea Kinematic Viscosity Data
count = 1;
for i = 1:length(navsea.K)
for ii = 1:length(navsea.pa)
navsea.kin_fit(count,1) = log(navsea.K(i));
navsea.kin_fit(count,2) = navsea.pa(ii);
navsea.kin_fit(count,3) = log(navsea.kinSI(ii,i));
count = count + 1;
end
end
navsea.P_kin = polyfitn([navsea.kin_fit(:,1) navsea.kin_fit(:,2)],navsea.kin_fit(:,3),2);
navsea.z_kin = polyvaln(navsea.P_kin,[log(navsea.xg(:)),navsea.yg(:)]);
fig_name = ['NavSea Kinematic Viscosity'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(navsea.kin_fit(:,1),navsea.kin_fit(:,2),navsea.kin_fit(:,3));
title(fig_name)
hold on;
surf(log(navsea.xg),navsea.yg,reshape(navsea.z_kin,size(navsea.xg)),'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('ln(Temp) [ln(K)]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('ln(Viscosity) [$ln(\frac{m^2}{s}$)]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Create Sargent Kinemtic Viscosity
sargent.pa = linspace(0,10000*6895,11).';
sargent.psig = sargent.pa / 6895; %pa to psig
% for each kin temp, fit a line through pressures
for i = 1:size(navsea.kinSI,2)
[sargent.navP(i,:), sargent.navS(i,:)] = polyfit(navsea.pa,navsea.kinSI(:,i),1);
end
%Calc Slope as function of y-intercept
sargent.kin_slope_eq = polyfit(sargent.navP(:,2),sargent.navP(:,1),1);
count = 1;
for i = 1:size(sargent.K,1)
for ii = 1:size(sargent.pa,1)
sargent.kin_fit(count,1) = log(sargent.K(i,1));
sargent.kin_fit(count,2) = sargent.pa(ii,1);
sargent.kin_slope = sargent.kin_slope_eq(1) * sargent.kinSI(i,1) + sargent.kin_slope_eq(2);
sargent.kin_fit(count,3) = log(sargent.kin_slope .* sargent.pa(ii,1) + sargent.kinSI(i,1));
count = count + 1;
end
end
sargent.P_kin = polyfitn([sargent.kin_fit(:,1) sargent.kin_fit(:,2)],sargent.kin_fit(:,3),2);
sargent.z_kin = polyvaln(sargent.P_kin,[log(sargent.xg(:)),sargent.yg(:)]);
sargent.P_kin.R2
fig_name = ['Sargent Kinematic Viscosity'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(sargent.kin_fit(:,1),sargent.kin_fit(:,2),sargent.kin_fit(:,3));
title(fig_name)
hold on;
surf(log(sargent.xg),sargent.yg,reshape(sargent.z_kin,size(sargent.xg)),'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('ln(Temp) [ln(K)]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('ln(Viscosity) [$ln(\frac{m^2}{s}$)]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Copy Navsea Density into Sargent
sargent.den_fit(:,1) = exp(sargent.kin_fit(:,1));
sargent.den_fit(:,2) = sargent.kin_fit(:,2);
sargent.den_fit(:,3) = polyvaln(navsea.P_den,[sargent.den_fit(:,1),sargent.den_fit(:,2)]);
sargent.P_den = polyfitn([sargent.den_fit(:,1) sargent.den_fit(:,2)],sargent.den_fit(:,3),2);
sargent.z_den = polyvaln(sargent.P_den,[sargent.xg(:),sargent.yg(:)]);
sargent.P_den.R2
fig_name = ['Sargent Density'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
title(fig_name)
hold on;
surf(sargent.xg,sargent.yg,reshape(sargent.z_den,size(sargent.xg)),'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Density [$\frac{kg}{m^3}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Create Sargent Dynamic Viscosity
sargent.dyn_fit(:,1:2) = sargent.kin_fit(:,1:2);
sargent.dyn_fit(:,3) = log(exp(sargent.kin_fit(:,3)) .* sargent.den_fit(:,3));
sargent.P_dyn = polyfitn([sargent.dyn_fit(:,1) sargent.dyn_fit(:,2)],sargent.dyn_fit(:,3),2);
sargent.z_dyn = polyvaln(sargent.P_dyn,[log(sargent.xg(:)),sargent.yg(:)]);
sargent.P_dyn.R2
fig_name = ['Sargent Dynamic Viscosity'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(sargent.dyn_fit(:,1),sargent.dyn_fit(:,2),sargent.dyn_fit(:,3));
title(fig_name)
hold on;
surf(log(sargent.xg),sargent.yg,reshape(sargent.z_dyn,size(sargent.xg)),'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('ln(Temp) [ln(K)]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('ln(Viscosity) [$ln(\frac{kg}{ms})$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Test Density Surf
sargent.denSURF = sargent.P_den.Coefficients(1).*sargent.xg.^2 ...
+ sargent.P_den.Coefficients(2).*sargent.xg.*sargent.yg ...
+ sargent.P_den.Coefficients(3).*sargent.xg ...
+ sargent.P_den.Coefficients(4).*sargent.yg.^2 ...
+ sargent.P_den.Coefficients(5).*sargent.yg ...
+ sargent.P_den.Coefficients(6);
fig_name = ['Sargent Density Check'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
title(fig_name)
surf(sargent.xg,sargent.yg,reshape(sargent.z_den,size(sargent.xg)),'EdgeColor','none')
hold on;
surf(sargent.xg,sargent.yg,sargent.denSURF,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Density [$\frac{kg}{m^3}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
% Calculate speed of sound
sargent.d_den_dp = ...
sargent.P_den.Coefficients(2) .* sargent.xg + ...
sargent.P_den.Coefficients(4) .*2.*sargent.yg + ...
sargent.P_den.Coefficients(5);
sargent.bulk = sargent.denSURF ./ sargent.d_den_dp;
sargent.c = sqrt(sargent.bulk ./ sargent.denSURF);
fig_name = ['Sargent Bulk Modulus'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
title(fig_name)
surf(sargent.xg,sargent.yg,reshape(sargent.bulk.*1E-9,size(sargent.xg)),'EdgeColor','none')
hold on;
surf(sargent.xg,sargent.yg,sargent.bulk.*1E-9,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Bulk Modulus [$Gpa$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Bulk Modulus Contour (SI)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
% title(fig_name)
% surf(sargent.xg,sargent.yg,reshape(sargent.bulk.*1E-9,size(sargent.xg)),'EdgeColor','none')
% hold on;
contour(sargent.xg,sargent.yg,sargent.bulk.*1E-9,20,'ShowText','on')
set(gca,'fontsize',font_size,'fontname',font_name);
title([fig_name ' [Gpa]'],'Interpreter','latex')
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Bulk Modulus [$Gpa$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Bulk Modulus Contour(English)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
% title(fig_name)
% surf(sargent.xg,sargent.yg,reshape(sargent.bulk.*1E-9,size(sargent.xg)),'EdgeColor','none')
% hold on;
sargent.xg_F = ((sargent.xg - 273.15) .* 9 ./ 5) + 32;
sargent.yg_psi = sargent.yg .* 0.00014503773;
sargent.bulk_psi = sargent.bulk .* 0.00014503773 .*1E-3;
contour(sargent.xg_F,sargent.yg_psi,sargent.bulk_psi,20,'ShowText','on')
set(gca,'fontsize',font_size,'fontname',font_name);
title([fig_name ' [KSI]'],'Interpreter','latex')
xlabel('Temp [F]','Interpreter','latex')
ylabel('pressure [psi]','Interpreter','latex')
zlabel('Bulk Modulus [$ksi$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Isothermal Wave speed (SI)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
title(fig_name)
surf(sargent.xg,sargent.yg,reshape(sargent.c,size(sargent.xg)),'EdgeColor','none')
hold on;
surf(sargent.xg,sargent.yg,sargent.c,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Wave speed [$\frac{m}{s}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Isothermal Wave speed contour (SI)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
% title(fig_name)
% surf(sargent.xg,sargent.yg,reshape(sargent.c,size(sargent.xg)),'EdgeColor','none')
% hold on;
contour(sargent.xg,sargent.yg,sargent.c,20,'ShowText','on')
set(gca,'fontsize',font_size,'fontname',font_name);
title([fig_name ' [$\frac{m}{s}$]'],'Interpreter','latex')
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Wave speed [$\frac{m}{s}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Isothermal Wave speed (English)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
% title(fig_name)
sargent.xg_F = ((sargent.xg - 273.15) .* 9 ./ 5) + 32;
sargent.yg_psi = sargent.yg .* 0.00014503773;
sargent.c_ft = sargent.c .* 3.2808398950131;
surf(sargent.xg_F,sargent.yg_psi,reshape(sargent.c_ft,size(sargent.xg_F)),'EdgeColor','none')
hold on;
surf(sargent.xg_F,sargent.yg_psi,sargent.c_ft,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [F]','Interpreter','latex')
ylabel('pressure [psi]','Interpreter','latex')
zlabel('Wave speed [$\frac{ft}{s}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Isothermal Wave speed contour (English)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
% title(fig_name)
sargent.xg_F = ((sargent.xg - 273.15) .* 9 ./ 5) + 32;
sargent.yg_psi = sargent.yg .* 0.00014503773;
sargent.c_ft = sargent.c .* 3.2808398950131;
% contour(sargent.xg_F,sargent.yg_psi,reshape(sargent.c_ft,size(sargent.xg_F)),'EdgeColor','none')
% hold on;
contour(sargent.xg_F,sargent.yg_psi,sargent.c_ft,20,'ShowText','on')
set(gca,'fontsize',font_size,'fontname',font_name);
title([fig_name ' [$\frac{ft}{s}$]'],'Interpreter','latex')
xlabel('Temp [F]','Interpreter','latex')
ylabel('pressure [psi]','Interpreter','latex')
zlabel('Wave speed [$\frac{ft}{s}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
fig_name = ['Sargent Lump Parameter length (20 kHz)(English)'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
% title(fig_name)
sargent.xg_F = ((sargent.xg - 273.15) .* 9 ./ 5) + 32;
sargent.yg_psi = sargent.yg .* 0.00014503773;
sargent.c_ft = sargent.c .* 3.2808398950131;
% contour(sargent.xg_F,sargent.yg_psi,reshape(sargent.c_ft,size(sargent.xg_F)),'EdgeColor','none')
% hold on;
contour(sargent.xg_F,sargent.yg_psi,sargent.c_ft./20000.*0.1.*12.*1000,20,'ShowText','on')
set(gca,'fontsize',font_size,'fontname',font_name);
title([fig_name ' [Thou]'],'Interpreter','latex')
xlabel('Temp [F]','Interpreter','latex')
ylabel('pressure [psi]','Interpreter','latex')
zlabel('Wave speed [$Thou$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
sum(sum(sargent.denSURF - reshape(sargent.z_den,size(sargent.xg))))
fig_name = ['Sargent Density Error'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
% scatter3(sargent.den_fit(:,1),sargent.den_fit(:,2),sargent.den_fit(:,3));
title(fig_name)
surf(sargent.xg,sargent.yg,sargent.denSURF - reshape(sargent.z_den,size(sargent.xg)),'EdgeColor','none')
hold on;
% surf(sargent.xg,sargent.yg,sargent.denSURF,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Density [$\frac{kg}{m^3}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Test Dynamic Viscosity SURF
sargent.dynSURF = exp(sargent.P_dyn.Coefficients(1).*log(sargent.xg).^2 ...
+ sargent.P_dyn.Coefficients(2).*log(sargent.xg).*sargent.yg ...
+ sargent.P_dyn.Coefficients(3).*log(sargent.xg) ...
+ sargent.P_dyn.Coefficients(4).*sargent.yg.^2 ...
+ sargent.P_dyn.Coefficients(5).*sargent.yg ...
+ sargent.P_dyn.Coefficients(6));
fig_name = ['Sargent Dynamic Viscosity Check'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
scatter3(exp(sargent.dyn_fit(:,1)),sargent.dyn_fit(:,2),exp(sargent.dyn_fit(:,3)));
title(fig_name)
hold on;
scatter3(exp(navsea.dyn_fit(:,1)),navsea.dyn_fit(:,2),exp(navsea.dyn_fit(:,3)),'r+');
surf(sargent.xg,sargent.yg,sargent.dynSURF,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('Temp [K]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('Viscosity [$\frac{kg}{ms}$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
% sum(sum(sargent.denSURF - reshape(sargent.z_den,size(sargent.xg))))
fig_name = ['Sargent Dynamic Viscosity Error'];
figure('name',fig_name,'position',[scrsz(3)/16 scrsz(4)/16 fig_width fig_height])
surf(sargent.xg,sargent.yg,sargent.dynSURF - reshape(sargent.z_dyn,size(sargent.xg)),'EdgeColor','none')
title(fig_name)
hold on;
% surf(log(sargent.xg),sargent.yg,sargent.dynSURF,'EdgeColor','none')
set(gca,'fontsize',font_size,'fontname',font_name);
xlabel('ln(Temp) [ln(K)]','Interpreter','latex')
ylabel('pressure [pa]','Interpreter','latex')
zlabel('ln(Viscosity) [$ln(\frac{kg}{ms})$]','Interpreter','latex')
saveas(gcf,[local '\figure\' fig_name],'png')
%% Write Sargent Data to Excel
sargent.kinSI = zeros(length(sargent.pa),length(sargent.K));
sargent.kinSI = exp(reshape(sargent.kin_fit(:,3),size(sargent.kinSI)));
sargent.kin = sargent.kinSI * 1E6; % m^2/s to Cst
xlswrite(sargent.tables,{'T[K]','p[pa]','m^2/s'},'kin_visc','A1');
xlswrite(sargent.tables,sargent.K.','kin_visc','B2');
xlswrite(sargent.tables,sargent.pa,'kin_visc','A3');
xlswrite(sargent.tables,sargent.kinSI,'kin_visc','B3');
xlswrite(sargent.tables,{'T[F]','p[psig]','Cst'},'kin_visc','H1');
xlswrite(sargent.tables,sargent.F.','kin_visc','I2');
xlswrite(sargent.tables,sargent.psig,'kin_visc','H3');
xlswrite(sargent.tables,sargent.kin,'kin_visc','I3');
sargent.navsea_K = navsea.K;
sargent.navsea_C = sargent.navsea_K - 273;
[sargent.navsea_xg,sargent.navsea_yg] = meshgrid(sargent.navsea_K,sargent.pa);
sargent.navsea_kin = polyvaln(sargent.P_kin,[log(sargent.navsea_xg(:)),sargent.navsea_yg(:)]);
sargent.navsea_kinSI = exp(reshape(sargent.navsea_kin,size(sargent.navsea_xg)));
sargent.navsea_kin = sargent.navsea_kinSI * 1E6; % m^2/s to Cst
xlswrite(sargent.tables,{'T[C]','p[psig]','Cst'},'kin_visc','O1');
xlswrite(sargent.tables,sargent.navsea_C,'kin_visc','P2');
xlswrite(sargent.tables,sargent.psig,'kin_visc','O3');
xlswrite(sargent.tables,sargent.navsea_kin,'kin_visc','P3');
sargent.dynSI = zeros(length(sargent.pa),length(sargent.K));
sargent.dynSI = exp(reshape(sargent.dyn_fit(:,3),size(sargent.dynSI)));
sargent.dyn = sargent.dynSI * 1000; % kg/m-s to cp
xlswrite(sargent.tables,{'T[K]','p[pa]','kg/m-s'},'dyn_visc','A1');
xlswrite(sargent.tables,sargent.K.','dyn_visc','B2');
xlswrite(sargent.tables,sargent.pa,'dyn_visc','A3');
xlswrite(sargent.tables,sargent.dynSI,'dyn_visc','B3');
xlswrite(sargent.tables,{'T[F]','p[psig]','cp'},'dyn_visc','H1');
xlswrite(sargent.tables,sargent.F.','dyn_visc','I2');
xlswrite(sargent.tables,sargent.psig,'dyn_visc','H3');
xlswrite(sargent.tables,sargent.dyn,'dyn_visc','I3');
% xlswrite(sargent.tables,sargent.K.','Density','B1');
% xlswrite(sargent.tables,sargent.pa,'Density','A2');
% xlswrite(sargent.tables,.denSI,'Density','B2');
sargent.denSI = zeros(length(sargent.pa),length(sargent.K));
sargent.denSI = reshape(sargent.den_fit(:,3),size(sargent.denSI));
sargent.den = sargent.denSI; %kg/m^3
xlswrite(sargent.tables,{'T[K]','p[pa]','kg/m^3'},'Density','A1');
xlswrite(sargent.tables,sargent.K.','Density','B2');
xlswrite(sargent.tables,sargent.pa,'Density','A3');
xlswrite(sargent.tables,sargent.denSI,'Density','B3');
xlswrite(sargent.tables,{'T[F]','p[psig]','kg/m^3'},'Density','H1');
xlswrite(sargent.tables,sargent.F.','Density','I2');
xlswrite(sargent.tables,sargent.psig,'Density','H3');
xlswrite(sargent.tables,sargent.den,'Density','I3');
%% Write to Density and Viscosity to Excel
% xlswrite(sargent.eqs,{'xg T[K]','yg p[pa]'},'Density','A1')
% xlswrite(sargent.eqs,sargent.P_den.ModelTerms,'Density','A2')
% xlswrite(sargent.eqs,sargent.P_den.Coefficients,'Density','A9')
% xlswrite(sargent.eqs,{'sargent.denSURF = sargent.P_den.Coefficients(1).*sargent.xg.^2 + sargent.P_den.Coefficients(2).*sargent.xg.*sargent.yg + sargent.P_den.Coefficients(3).*sargent.xg + sargent.P_den.Coefficients(4).*sargent.yg.^2 + sargent.P_den.Coefficients(5).*sargent.yg + sargent.P_den.Coefficients(6)'},'Density','A11')
%
% xlswrite(sargent.eqs,{'xg T[K]','yg p[pa]'},'Dynamic','A1')
% xlswrite(sargent.eqs,sargent.P_dyn.ModelTerms,'Dynamic','A2')
% xlswrite(sargent.eqs,sargent.P_dyn.Coefficients,'Dynamic','A9')
% xlswrite(sargent.eqs,{'sargent.dynSURF = exp(sargent.P_dyn.Coefficients(1).*log(sargent.xg).^2 + sargent.P_dyn.Coefficients(2).*log(sargent.xg).*sargent.yg + sargent.P_dyn.Coefficients(3).*log(sargent.xg) + sargent.P_dyn.Coefficients(4).*sargent.yg.^2 + sargent.P_dyn.Coefficients(5).*sargent.yg + sargent.P_dyn.Coefficients(6))'},'Dynamic','A11')
%% Write to Txt file
%temp,press,mu,den
format longE
fileid = fopen(sargent.txt,'w');
fprintf(fileid,'WARNING: IGNORING MORLINA AND AERATION EFFECTS \r\n \r\n');
fprintf(fileid,'Equation for Density[m^3] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['density = ' ...
num2str(sargent.P_den.Coefficients(1),16) '*temp*temp + ' ...
num2str(sargent.P_den.Coefficients(2),16) '*temp*press + ' ...
num2str(sargent.P_den.Coefficients(3),16) '*temp + ' ...
num2str(sargent.P_den.Coefficients(4),16) '*press*press + ' ...
num2str(sargent.P_den.Coefficients(5),16) '*press + ' ...
num2str(sargent.P_den.Coefficients(6),16) '; \r\n \r\n']);
fprintf(fileid,'Equation for Dyanmic Viscosity [kg/ms] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['mu_lam = exp(' ...
num2str(sargent.P_dyn.Coefficients(1),16) '*log(temp)*log(temp) + ' ...
num2str(sargent.P_dyn.Coefficients(2),16) '*log(temp)*press + ' ...
num2str(sargent.P_dyn.Coefficients(3),16) '*log(temp) + ' ...
num2str(sargent.P_dyn.Coefficients(4),16) '*press^2 + ' ...
num2str(sargent.P_dyn.Coefficients(5),16) '*press + ' ...
num2str(sargent.P_dyn.Coefficients(6),16) '); \r\n \r\n']);
fprintf(fileid,'Equation for Density pressure partial derivative [m^3/pa] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['d_den_dp = ' ...
num2str(sargent.P_den.Coefficients(2),16) '*temp + ' ...
num2str(sargent.P_den.Coefficients(4),16) '*2*press + ' ...
num2str(sargent.P_den.Coefficients(5),16) '; \r\n \r\n']);
fprintf(fileid,'Equation for Isothermal tangent Bulk Modulus [pa] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['K = density / d_den_dp; \r\n \r\n']);
fprintf(fileid,'Equation for Isothermal phase speed [m/s] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['c = sqrt(K / density); \r\n \r\n']);
fprintf(fileid,'Equation for Density temperature partial derivative [m^3/pa] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['d_den_dT = ' ...
num2str(sargent.P_den.Coefficients(1),16) '*2*temp + ' ...
num2str(sargent.P_den.Coefficients(2),16) '*press + ' ...
num2str(sargent.P_den.Coefficients(3),16) '; \r\n \r\n']);
fprintf(fileid,'Equation for cubical expansion coefficient [m^3/K] from Temp[K], pressure[pa] \r\n');
fprintf(fileid,['alpha = density / d_den_dT; \r\n \r\n']);
fclose(fileid);
Page updated
Google Sites
Report abuse