Modelling of the thermal field at dissipation of ultrasonic acoustic oscillations in viscous medium

This example demonstrates the modelling of ultrasonic radiation propagation in a viscous medium. Damping of acoustic oscillations and heating of the medium due to distributed heat sources are modelled.

The energy carried by an ultrasonic beam decays as it passes through a viscous medium. If at the origin at $x=0$ the intensity of a plane travelling wave is equal to $I_0$, then at a distance x it decreases to the value $I(x)$: $$I(x) = I_0 \cdot e^{-kx},$$ where (k) is the intensity damping factor. Thus, the energy loss of the original ultrasonic beam during its propagation over a unit distance, $\frac{dI}{dx}\$, is equal to $k \cdot I$. The attenuation coefficient, $k$, consists of two components, one due to absorption and the other due to scattering. The energy scattered from the main beam can be absorbed in other regions of the viscous medium.

Using simple heat propagation equations, the expected heating can be calculated. The temperature distribution of a viscous medium is described by Eq: $$\frac{\partial T(x,t)}{\partial \tau} = \frac{\partial}{\partial x} \frac{\partial T(x,t)}{\partial x} + \frac{q_0}{c\gamma} e^{-kx}$$

where $T$ is the temperature, $t$ is the exposure time, $x$ is the coordinate along the axis of radiation propagation, $q_0$ is the initial value of the power of volumetric heat sources, $c$ is the heat capacity of the medium, $\gamma$ is the density of the medium. Power of volumetric sources as a function of coordinate: $q(x) = I_0 k e^{-kx}.$

Calculation of the temperature field

Connection of necessary libraries:

using .DifferentialEquations, Plots
plotly()

Plots.PlotlyBackend()

Determination of radiation and viscous medium parameters:

α = 0.5 # коэффициент теплопроводности Вт / (м*К)
c = 3360.0 # теплоёмкость, Дж/кг*К
γ = 1050.0 # плотность, кг/м3
L = 0.1 # глубина рассматриваемого участка, м
k = 0.35 # Нп/м
I0 = 10000.0 # интенсивность излучения Вт / м2
F = 0.0001 #м2
m = 0.001; #кг

Realisation of the temperature calculation function by time and coordinate:

function heat_equation!(du, u, p, t)
    α, I0, c, γ, k, x = p
    dx = x[2] - x[1]
    N = length(x)
    # Первый участок
    du[1] = α * F * (1 / (c * m)) * (u[2] - 2u[1] + 320) / dx^2 + I0 / (c * γ) * exp(-k*x[1])
    # Промежуточные участки
    for i in 2:(N-1)
        du[i] = α * F * (1 / (c * m)) * (u[i+1] - 2u[i] + u[i-1]) / dx^2 + I0 / (c * γ) * exp(-k*x[i])
    end
    # Конечный участок
    du[N] = α * F * (1 / (c * m)) * (u[N-1] - u[N]) / dx^2 + I0 / (c * γ) * exp(-k*x[N])
end

heat_equation! (generic function with 1 method)

Defining the calculation parameters:

x = range(0, stop=1, length=100)  # задаем координаты
tspan = (0.0, 2000.0)  # временной интервал
u0 = 310.0 .* ones(length(x))  # начальные условия
p = (0.5, I0, 3360.0, 1050.0, k, x)  # параметры (α, I0, c, γ, k, L, x)  

(0.5, 10000.0, 3360.0, 1050.0, 0.35, 0.0:0.010101010101010102:1.0)

Realisation of the coordinate-dependent heat flux calculation function:

function teplo(x)
    return I0 / (c * γ) * exp.(-k*x)
end

teplo (generic function with 1 method)

Visualisation of the heat flux value as a function of coordinate, taking into account the attenuation coefficient:

plot(teplo(x))

Solving the problem using DifferentialEquations library methods:

prob = ODEProblem(heat_equation!, u0, tspan, p)
sol = solve(prob);

Determination of the temperature field at the final step of the calculation:

t = 2000.0
idx = argmin(abs.(sol.t .- t))
u_at_t = sol.u[idx];

Visualisation of the temperature field at the final calculation step:

plot(u_at_t)

Conclusions

This case study demonstrated the modelling of the propagation of acoustic oscillations in a viscous medium, as well as their damping and energy conversion to heat. The temperature field of a semi-infinite plate has been constructed and the methods of the DifferentialEquations library for solving systems of differential equations have been demonstrated.