Trying to solve a second order differential equation. Your advices

Specific Domains Optimization (Mathematical)
5

Oh, this is a debugging technique I always hope will work in Julia, and I think it does here! The Unitful.jl package provides the tools to define quantities with units in the type. I added types to all your parameters, evaluated the expression to be integrated, and got this error:

using OrdinaryDiffEq, Plots, Unitful                                             
                                                                                 
#SI units                                                                        
#external parameters                                                             
Pₕ = 1.01325e5u"Pa"                  #hydrostatic pressure (Pa)                  
                                                                                 
#Parameters (liquid = water)                                                     
Pᵥ = 2.5e3u"Pa"                      #vapour pressure (Pa)                       
σ = 73.0e-3u"N/m"                     #surface tension (N/m)                     
ρ = 1000.0u"kg*m^-3"                      #volumic mass (Kg.m-3)                 
η = 1.0e-3u"Pa*s"                      #dynamic viscosity (Pa.s)                 
                                                                                 
#Parameter (gas = Ar)                                                            
γ = 1.666                       #polytropic ratio (-)                               
                                                                                    
#Acoustic Parameters                                                             
Pₐ = 1.3e5u"Pa"                      #acoustic pressure (Pa)                     
fₐ = 2.0e4u"Hz"                      #acoustic frequency (s-1)                   
                                                                                 
#Initial conditions                                                              
R₀ = [4.0e-6u"m"]                   #initial radius (m)                             
dR₀ = [1.0e-2u"m/s"]                  #initial radial velocity (ms-1)            
tspan = (0.0u"s",150.0e-6u"s")          #(s) -- 1 acoustic cycle = 50 μs         
                                                                                 
                                                                                 
function test_rpnnp(dR, R, p, t)
    @. (1/ρ*R)*((Pₕ-Pᵥ-(2*σ/R₀))*(R₀/R)^(3*γ)-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end

P = t->Pₐ*sin(2*pi*fₐ*t + pi/2)

test_rpnnp(dR₀, R₀, P, 0u"s")

ERROR: DimensionError: -0.0008120400000000001 m^3 s^-2 and 37.5 m s^-2 are not dimensionally compatible.

I think the first term needs a divide-by-area.

Unfortunately, OrdinaryDiffEq does not support units internally since it assumes velocity and acceleration have the same type, so you will need to go back to your version with no units to actually run the solver.

1 Like