Моделирование теплового поля при рассеивании ультразвуковых акустических колебаний в вязкой среде
Данный пример демонстрирует моделирование распространения ультразвукового излучения в вязкой среде. Моделируется затухание акустических колебаний и нагрев среды засчёт распределённых источников тепла.
Энергия переносимая ультразвуковым пучком, затухает при прохождении её через вязкую среду. Если в начале координат при интенсивность плоской бегущей волны равна , то на расстоянии x она уменьшается до величины :
где (k) – коэффициент затухания по интенсивности. Таким образом, потери энергии исходного ультразвукового пучка при его распространении на единичное расстояние, $\frac{dI}{dx}$, равны . Коэффициент затухания, , состоит из двух компонент, одной, обусловленной поглощением, и другой, обусловленной рассеянием. Энергия, рассеянная из основного пучка, может поглощаться в других областях вязкой среды.
Используя простые уравнения распространения тепла, можно рассчитать ожидаемый нагрев. Распределение температуры вязкой среды описывается уравнением:
где – температура, – время воздействия, – координата по оси распространения излучения, - начальное значение мощности объёмных источников тепла, – теплоёмкость среды, - плотность среды. Мощность объёмных источников в зависимости от координаты:
Расчёт температурного поля
Подключение необходимых библиотек:
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 для решения систем дифференциальных уравнений.