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

Скрипт Open Live Script

Этот пример показывает квазистатический осесимметричный тепловой рабочий процесс расчета напряжений путем репродуцирования результатов упрощенной модели дискового тормоза, обсужденной в [1]. Дисковые тормоза поглощают механическую энергию посредством трения и преобразовывают его в тепловую энергию, которая затем рассеивается. Пример использует упрощенную модель дискового тормоза в одном процессе торможения от постоянной начальной угловой скорости до бездействия. Рабочий процесс имеет два шага:

Переходный тепловой анализ, чтобы вычислить температурное распределение в диске с помощью тепла течет от тормозных колодок
Квазистатический структурный анализ, чтобы вычислить тепловые усилия в несколько раз решения с помощью ранее полученного температурного распределения, чтобы задать тепловые нагрузки

Получившиеся графики показывают температурное распределение, радиальное напряжение, напряжение обруча и напряжение фон Мизеса в течение соответствующих времен решения.

Свойства дискового тормоза и геометрия

На основе предположений, используемых в [1], пример уменьшает аналитическую область до прямоугольной области, соответствующей осесимметричному разделу кольцевого диска. Из-за геометрической симметрии и симметрии загрузки диска, модели в качестве примера только половина толщины диска и эффекта одной клавиатуры. В следующем рисунке левый край соответствует внутреннему радиусу диска $r_{d}$ . Правый край соответствует внешнему радиусу диска $R_{d}$ и также совпадает с внешним радиусом клавиатуры $R_{p}$ . Диск испытывает давление клавиатуры, которая генерирует поток тепла. Вместо того, чтобы моделировать клавиатуру явным образом, включайте ее эффект в тепловой анализ путем определения этого потока тепла как граничного условия от внутреннего радиуса клавиатуры $r_{p}$ к внешнему радиусу клавиатуры $R_{p}$ .

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

Создайте переходную осесимметричную тепловую модель.

modelT = createpde('thermal','transient-axisymmetric');

Создайте геометрию с двумя смежными прямоугольниками. Верхний край более длинного прямоугольника (справа) представляет область контакта клавиатуры диска.

R1 = [3,4, [  66,  76.5,  76.5,   66, -5.5, -5.5, 0, 0]/1000]';
R2 = [3,4, [76.5, 113.5, 113.5, 76.5, -5.5, -5.5, 0, 0]/1000]';

gdm = [R1 R2];
ns = char('R1','R2');
g = decsg(gdm,'R1 + R2',ns');

Присвойте геометрию тепловой модели.

geometryFromEdges(modelT,g);

Постройте геометрию с ребром и столкнитесь с метками.

figure
pdegplot(modelT,'EdgeLabels','on','FaceLabels','on')

Figure contains an axes object. The axes object contains 10 objects of type line, text.

Сгенерируйте mesh. Чтобы совпадать с mesh, используемой в [1], используйте линейный геометрический порядок вместо квадратичного порядка по умолчанию.

generateMesh(modelT,'Hmax',0.5E-04,'GeometricOrder','linear');

Задайте тепловые свойства материала диска.

alphad = 1.44E-5; % Diffusivity of disc
Kd = 51;
rhod = 7100;
cpd = Kd/rhod/alphad;
thermalProperties(modelT,'ThermalConductivity',Kd, ...
                         'MassDensity',rhod, ...
                         'SpecificHeat',cpd);

Задайте граничное условие потока тепла с учетом области клавиатуры. Для определения qFcn функционируйте, смотрите Функцию Потока Тепла.

thermalBC(modelT,'Edge',6,'HeatFlux',@qFcn);

Установите начальную температуру.

thermalIC(modelT,20);

Решите модель в течение времен, используемых в [1].

tlist = [0 0.1 0.2 1.0 2.0 3.0 3.96];
Rt = solve(modelT,tlist);

Постройте температурное изменение со временем в трех ключевых радиальных местах. Получившийся график сопоставим с графиком, полученным в [1].

iTRd = interpolateTemperature(Rt,[0.1135;0],1:numel(Rt.SolutionTimes));
iTrp = interpolateTemperature(Rt,[0.0765;0],1:numel(Rt.SolutionTimes));
iTrd = interpolateTemperature(Rt,[0.066;0],1:numel(Rt.SolutionTimes));

figure
plot(tlist,iTRd)
hold on
plot(tlist,iTrp)
plot(tlist,iTrd)
title('Temperature Variation with Time at Key Radial Locations')
legend('R_d','r_p','r_d')
xlabel 't, s'
ylabel 'T,^{\circ}C'

Figure contains an axes object. The axes object with title Temperature Variation with Time at Key Radial Locations contains 3 objects of type line. These objects represent R_d, r_p, r_d.

Структурный анализ: вычислите тепловое напряжение

Создайте осесимметричную статическую модель структурного анализа.

model = createpde('structural','static-axisymmetric');

Присвойте геометрию и mesh, используемую для тепловой модели.

model.Geometry = modelT.Geometry;
model.Mesh = modelT.Mesh;

Задайте структурные свойства диска.

structuralProperties(model,'YoungsModulus',99.97E9, ...
                           'PoissonsRatio',0.29, ...
                           'CTE',1.08E-5);

Ограничьте модель предотвращать твердое движение.

structuralBC(model,'Edge',[3,4],'ZDisplacement',0);

Задайте ссылочную температуру, которая соответствует состоянию нулевого теплового напряжения модели.

model.ReferenceTemperature = 20;

Задайте тепловую нагрузку при помощи переходных тепловых результатов Rt. Времена решения эквивалентны в тепловом анализе модели. В течение каждого раза решения решите соответствующую статическую задачу структурного анализа и постройте температурное распределение, радиальное напряжение, напряжение обруча и напряжение фон Мизеса. Для определения plotResults функционируйте, смотрите Функцию Результатов Графика. Результаты сопоставимы с фигурой 5 от [1].

for n = 2:numel(Rt.SolutionTimes)
structuralBodyLoad(model,'Temperature',Rt,'TimeStep',n);
R = solve(model);
plotResults(model,R,modelT,Rt,n);
end

Нагрейте функцию потока

Эта функция помощника вычисляет переходное значение потока тепла от клавиатуры до диска. Это использует эмпирическую формулу от [1].

function q = qFcn(r,s)
alphad = 1.44E-5; % Diffusivity of disc
Kd = 51; % Conductivity of disc
rhod = 7100; % Density of disc
cpd = Kd/rhod/alphad; % Specific heat capacity of disc

alphap = 1.46E-5; % Diffusivity of pad
Kp = 34.3; % Conductivity of pad
rhop = 4700; % Density of pad
cpp = Kp/rhop/alphap; % Specific heat capacity of pad

f = 0.5; % Coefficient of friction
omega0 = 88.464; % Initial angular velocity
ts = 3.96; % Stopping time
p0 = 1.47E6*(64.5/360); % Pressure only spans 64.5 deg occupied by pad

omegat = omega0*(1 - s.time/ts); % Angular speed over time

eta = sqrt(Kd*rhod*cpd)/(sqrt(Kd*rhod*cpd) + sqrt(Kp*rhop*cpp));
q = (eta)*f*omegat*r.r*p0;
end

Постройте функцию результатов

Этот помощник графики функций температурное распределение, радиальное напряжение, напряжение обруча и напряжение фон Мизеса.

function plotResults(model,R,modelT,Rt,tID)
figure
subplot(2,2,1)
pdeplot(modelT,'XYData',Rt.Temperature(:,tID), ...
               'ColorMap','jet','Contour','on')
title({'Temperature'; ...
      ['max = ' num2str(max(Rt.Temperature(:,tID))) '^{\circ}C']})
xlabel 'r, m'
ylabel 'z, m'

subplot(2,2,2)
pdeplot(model,'XYData',R.Stress.srr, ...
              'ColorMap','jet','Contour','on')
title({'Radial Stress'; ...
       ['min = ' num2str(min(R.Stress.srr)/1E6,'%3.2f') ' MPa']; ...
       ['max = ' num2str(max(R.Stress.srr)/1E6,'%3.2f') ' MPa']})
xlabel 'r, m'
ylabel 'z, m'

subplot(2,2,3)
pdeplot(model,'XYData',R.Stress.sh, ...
              'ColorMap','jet','Contour','on')
title({'Hoop Stress'; ...
      ['min = ' num2str(min(R.Stress.sh)/1E6,'%3.2f') ' MPa']; ...
      ['max = ' num2str(max(R.Stress.sh)/1E6,'%3.2f') ' MPa']})
xlabel 'r, m'
ylabel 'z, m'

subplot(2,2,4)
pdeplot(model,'XYData',R.VonMisesStress, ...
              'ColorMap','jet','Contour','on')
title({'Von Mises Stress'; ...
      ['max = ' num2str(max(R.VonMisesStress)/1E6,'%3.2f') ' MPa']})
xlabel 'r, m'
ylabel 'z, m'

sgtitle(['Time = ' num2str(Rt.SolutionTimes(tID)) ' s'])
end

Ссылки

[1] Adamowicz, Адам. "Осесимметричная Модель FE к Анализу Тепловых Усилий в Тормозном диске". Журнал Теоретической и Прикладной Механики 53, выпуск 2 (апрель 2015): 357–370. https://doi.org/10.15632/jtam-pl.53.2.357.

Документация