Моделирование теплового поля при рассеивании ультразвуковых акустических колебаний в вязкой среде

Открыть пример в Engee

Данный пример демонстрирует моделирование распространения ультразвукового излучения в вязкой среде. Моделируется затухание акустических колебаний и нагрев среды засчёт распределённых источников тепла.

Энергия переносимая ультразвуковым пучком, затухает при прохождении её через вязкую среду. Если в начале координат при $x=0$ интенсивность плоской бегущей волны равна $I_0$, то на расстоянии x она уменьшается до величины $I(x)$: $$I(x) = I_0 \cdot e^{-kx},$$ где (k) – коэффициент затухания по интенсивности. Таким образом, потери энергии исходного ультразвукового пучка при его распространении на единичное расстояние, $\frac{dI}{dx}\$, равны $k \cdot I$. Коэффициент затухания, $k$, состоит из двух компонент, одной, обусловленной поглощением, и другой, обусловленной рассеянием. Энергия, рассеянная из основного пучка, может поглощаться в других областях вязкой среды.

Используя простые уравнения распространения тепла, можно рассчитать ожидаемый нагрев. Распределение температуры вязкой среды описывается уравнением: $$\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}$$

где $T$ – температура, $t$ – время воздействия, $x$ – координата по оси распространения излучения, $q_0$ - начальное значение мощности объёмных источников тепла, $c$ – теплоёмкость среды, $\gamma$ - плотность среды. Мощность объёмных источников в зависимости от координаты: $q(x) = I_0 k e^{-kx}.$

Расчёт температурного поля

Подключение необходимых библиотек:

using .DifferentialEquations, Plots
plotly()

Plots.PlotlyBackend()

Определение параметров излучения и вязкой среды:

α = 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; #кг

Реализация функции расчёта температуры по времени и координате:

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)

Определение параметров расчёта:

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)

Реализация функции расчёта теплового потока в зависимости от координаты:

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

teplo (generic function with 1 method)

Визуализация величины теплового потока в зависимости от координаты с учётом коэффициента затухания:

plot(teplo(x))

Решение задачи методами библиотеки DifferentialEquations:

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

Определение температурного поля на конечном шаге расчёта:

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

Визуализация температурного поля на конечном шаге расчёта:

plot(u_at_t)

Выводы

В данном примере было продемонстрировано моделирование распространения акустических колебаний в вязкой среде, а также их затухание и преобразование энергии в тепло. Было построено температурное поле полубесконечной пластины, а также продемонстрированы методы библиотеки DifferentialEquations для решения систем дифференциальных уравнений.