Hello all. I have a sinusoidal excitational model of the hydraulic pumps.
I have an expression which calculates the density, viscosity and bulk modulus as a fluid(P,T) ,P meaning pressure is a variable.
function fluid(Pre::Float64,Temp::Float64)
χ₀=0.001; # Air content percentage in the oil at STP
Pₛₐₜ=4e6; # Saturation pressure of the oil at stp
Pᵥₐₚ=1e4; # Vapor pressure of the oil at stp (Shell's data)
αₜ=7.2e-4; # Volumetric thermal expansion coeff. (/K)
ρₗ₀=872.0; # oil density at stp
ρₐ₀=1.225; # air density at stp
νₚₐ=60.3095; # Kinematic viscosity of liquid oil at operating temperature.(cSt) (Solved by params.m file)
bulk_liq=(10^((0.3766*(log10(νₚₐ)^0.3307))-0.2766)
+(5.851-0.01382*Temp)*(Pre/1e9))*1e9;
μₐ=(1.458e-6*(Temp+273.15)^1.5)/(Temp+110.4+273.15); #Sutherland Expression for air
μₗ₀=(νₚₐ*1e-6)*ρₗ₀*exp(-αₜ*(Temp-T_atm::Float64)); #Dynamic viscosity of liquid oil
if Pre>=Pₛₐₜ
rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)*exp((Pre-p_atm::Float64 )/bulk_liq-αₜ*(Temp-T_atm::Float64));
Beta_f=bulk_liq;
mu_f=(νₚₐ*1e-6)*(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)*exp((Pre-p_atm::Float64 )/bulk_liq-αₜ*(Temp-T_atm::Float64));
elseif Pre<Pₛₐₜ && Pre>Pᵥₐₚ
phi=(Pₛₐₜ-Pre)/(Pₛₐₜ-Pᵥₐₚ);
rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)/
((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*(p_atm::Float64 /Pre)^(1/1.42))+exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
#Beta_f=(1+(χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*((p_atm::Float64 /Pre)^(1/1.42))*phi)/
#(bulk_liq^-1 + (χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*(p_atm::Float64 /Pre)^(1/1.42)*(1/(Pₛₐₜ-p_atm::Float64 ))*((Pₛₐₜ-Pre)/1.42/Pre+1) );
Beta_f=((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*((p_atm::Float64 /Pre)^(1/1.42)))+exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)))/
(((χ₀)/(1-χ₀)*((Temp+273.15)/(T_atm::Float64 +273.15))*(p_atm::Float64 /Pre)^(1/1.42))*((1/(Pₛₐₜ-Pᵥₐₚ))+(phi/(1.42*Pre)))+
((exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)))/bulk_liq));
mu_f=((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pre)^(1/1.42))*μₐ)+
(exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64))*μₗ₀))/((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pre)^(1/1.42)))+
exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
elseif Pre<=Pᵥₐₚ
rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)/
((χ₀/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pᵥₐₚ)^(1/1.42)))+exp((p_atm::Float64 -Pᵥₐₚ)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
Beta_f=1.42*Pre;
mu_f=μₐ;
end
return rho_f,mu_f,Beta_f
end
Creating my ode function, I did this:
function damper!(dPdt,P,pars,t)
dₚ,dᵣ,Aₚ,Aᵣ,Aₜₚ,hₚ,hᵣ,Δₚ,Δᵣ,lₚₛ,δₚₛ,lᵣₛ,δᵣₛ,
Hₜ,
dₒco,dₒre,dₒba,lₒco,lₒre,lₒba,Cfmax_co,Cfmax_re,Cfmax_ba,
rᵥₐco,lᵥₐco,rₑₓₜco,Aₑₓₜco,kᵥₐco,x₀co,xₛₜᵣco,lₛₚₒco,δco,
rᵥₐre,lᵥₐre,rₑₓₜre,Aₑₓₜre,kᵥₐre,x₀re,xₛₜᵣre,lₛₚₒre,δre,
rᵥₐba,lᵥₐba,rₑₓₜba,Aₑₓₜba,kᵥₐba,x₀ba,xₛₜᵣba,lₛₚₒba,δba,
dₚₐ,lₚₐ,Cfmaxₕ,d_f₁,d_f₂,A_bot,Aₜₒₚ,Cfₗₐₙ,kᵢₙ,x₀,
V₀co,V₀re,Vg₀,Pg₀,Vres₀,
amp,freq,=pars;
#--Pressure in the 3 Chambers--#
Pᵣₑ,Pcₒ,Pᵣₑₛ=P; # Pressure in the rebound chamber,compression chamber and reservoir
#--Imposed Sinusodial Excitation-#
x=amp*sin(2*pi*freq*t); # Displacement from the middle of the piston
u=amp*2*pi*freq*cos(2*pi*freq*t); # Velocity of the piston rod
ρ₁,μ₁,β₁=fluid(Pᵣₑ,T);
ρ₂,μ₂,β₂=fluid(Pcₒ,T);
ρ₃,μ₃,β₃=fluid(Pᵣₑₛ,T);
Q_12=ori(Pᵣₑ,Pcₒ,dₒco,lₒco,Cfmax_co,ρ₁,μ₁)+
ori(Pᵣₑ,Pcₒ,dₒre,lₒre,Cfmax_re,ρ₁,μ₁)+
sum(spring.(Pᵣₑ,Pcₒ,ρ₁,rᵥₐre,lᵥₐre,rₑₓₜre,Aₑₓₜre,kᵥₐre,x₀re,xₛₜᵣre,lₛₚₒre,δre,3))+
pi*dₚ*δₚₛ^3/12/lₚₛ/μ₁*abs(Pᵣₑ-Pcₒ)*(Pᵣₑ>=Pcₒ); #Leakage of the piston seals
Q_21=ori(Pcₒ,Pᵣₑ,dₒco,lₒco,Cfmax_co,ρ₂,μ₂)+
ori(Pcₒ,Pᵣₑ,dₒre,lₒre,Cfmax_re,ρ₂,μ₂)+
sum(spring.(Pcₒ,Pᵣₑ,ρ₂,rᵥₐco,lᵥₐco,rₑₓₜco,Aₑₓₜco,kᵥₐco,x₀co,xₛₜᵣco,lₛₚₒco,δco,3))+
pi*dₚ*δₚₛ^3/12/lₚₛ/μ₂*abs(Pcₒ-Pᵣₑ)*(Pcₒ>=Pᵣₑ); #Leakage of the piston seals
Q_23=ori(Pcₒ,Pᵣₑₛ,dₒba,lₒba,Cfmax_ba,ρ₂,μ₂)+
sum(spring.(Pcₒ,Pᵣₑₛ,ρ₂,rᵥₐba,lᵥₐba,rₑₓₜba,Aₑₓₜba,kᵥₐba,x₀ba,xₛₜᵣba,lₛₚₒba,δba,3));
Q_32=ori(Pᵣₑₛ,Pcₒ,dₒba,lₒba,Cfmax_ba,ρ₃,μ₃)+
check(Pᵣₑₛ,Pcₒ,Pg₀,dₚₐ,lₚₐ,Cfmaxₕ,d_f₁,d_f₂,A_bot,Aₜₒₚ,Cfₗₐₙ,kᵢₙ,x₀);
dPdt[1]=β₁/(V₀re-Aₜₚ*x)*(Q_21*(ρ₂/ρ₁)-Q_12
+Aₜₚ*u);
dPdt[2]=β₂/(V₀co+Aₚ*x)*(Q_12*(ρ₁/ρ₂)-Q_21+Q_32*(ρ₃/ρ₂)-Q_23
-Aₚ*u);
dPdt[3]=(Q_23*(ρ₂/ρ₃)-Q_32)/
(((Vres₀+Vg₀-((Pg₀/Pᵣₑₛ)^(1/1.42))*Vg₀)/β₃)+(Pg₀^(1/1.42)*Vg₀/1.42/Pᵣₑₛ^(1+1/1.42)));
nothing
#Force=P_ex*A_tp-P_co*A_p+pi*vis_dyna((P_co+P_ex)/2)*(d_p*h_p/Delta_p+0.5*d_r*h_gu/Delta_r)*abs(u);
end
While solving the ode with Rodas5(), I get an error with a huge page, few lines of it being:
MethodError: no method matching fluid(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}, ::Float64)
Closest candidates are:
fluid(::Float64, ::Float64) at In[89]:1
Stacktrace:
[1] damper!(dPdt::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, P::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, pars::Vector{Any}, t::Float64)
@ Main .\In[94]:20
[2] (::ODEFunction{true, SciMLBase.AutoSpecialize, typeof(damper!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, ::Vararg{Any})
I guess I called the fluid function in a matlab way which would give real values of density, viscosity and bulk modulus. I have checked coupled of tutorials with this kind of a problem where a parameter is function of the variable but I didn’t get a grasp.
Your help will be highly appreciated.
