Анализ на собственные частоты с учетом температурных эффектов

22/05/2017

В некоторых случаях, особенно в контекте МЭМС-приложений, важно знать чувствительность собственных частот устройства к изменению температуры. В этой статье мы покажем, как данное исследование можно провести в программном пакете COMSOL Multiphysics® начиная с версии 5.3. Мы также обсудим эффекты, вызванные снижением жёсткости под механическим напряжением, геометрическими деформациями и учетом температурных зависимостей свойств материала.

Исследование собственных частот с учетом температурных эффектов

К некоторым устройствам предъявляются повышенные требования по стабильности частот при изменении условий окружающей среды. Наиболее распространённым параметром является температура. Однако, изменение условий также может быть связано с гигроскопическим расширением из-за изменения влажности. В устройствах сверхвысокой точности, требования к стабильности частот могут достигать чувствительности в одну миллиардную часть (parts-per-billion или ppb), т.е. 10-9. Создание и настройка подобных моделей, в которых будут учитываться даже такие небольшие эффекты, может стать довольно сложной задачей, поскольку необходимо учитывать взаимосвязь нескольких явлений.

Пример: расчёт балки прямоугольного сечения

Давайте рассмотрим балку квадратного сечения со следующими параметрами:

Свойство Символьное обозначение Значение
Длина L 10 мм
Ширина a 1 мм
Высота b 0.5 мм
Модуль Юнга E 100 ГПа
Коэффициент Пуассона ν 0
Плотность ρ 1000 кг/м3
Коэффициент теплового расширения в направлении оси x αx 1·10-5 1/K
Коэффициент теплового расширения в направлении оси y αy 2·10-5 1/K
Коэффициент теплового расширения в направлении оси z αz 3·10-5 1/K
Температурный сдвиг ΔT 10 K

Геометрия и сетка исследуемой балки.
Геометрия и сетка рассматриваемого примера балки.

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

Для анализа эффектов теплового расширения выберем исследование Eigenfrequency, Prestressed (Анализ на собственные частоты, с преднапряжением).

Скриншот исследования Eigenfrequency, Prestressed.
Добавление исследования Eigenfrequency, Prestressed.

Данное исследование состоит из двух шагов:

  1. Шаг Stationary (Стационарный), в рамках которого рассчитываются деформации и напряжения, вызванные тепловым расширением
  2. Шаг Eigenfrequency (Исследование на собственные частоты), в котором используется результат предыдущего шага.

Скриншот окна Построителя моделей COMSOL Multiphysics® с исследованием Eigenfrequency, Prestressed, которое содержит два шага.
Скриншот двух шагов исследования Eigenfrequency, Prestressed в дереве модели.

Для расчёта референсного решения можно добавить ещё одно исследование Eigenfrequency, либо выполнить ту же последовательность исследований, но без учета теплового расширения.

Собственные частоты балки, закреплённой с двух сторон, и консольной балки, закрепленной с одной стороны

Собственные частоты балки были рассчитаны для двух разных постановок (граничных условий):

  1. Балка, закрепленная с двух сторон
  2. Консольная балка, закрепленая с одной стороны и свободная с другой стороны (Кантилевер)

Результаты расчёта балки, закреплённой с двух сторон, приведены ниже.

Тип моды Собственная частота, Гц Собственная частота, Гц при ΔT = 10 K Отношение частот
Первая изгибная, в направлении оси z 50713.9 50425.1 0.9943
Первая изгибная, в направлении оси y 97659.6 97526.2 0.9986
Первая крутильная 266902 266917 1.00006
Первая продольная 500000 500025 1.00005

Изображение форм колебаний для балки, закреплённой с двух сторон.

Формы мод колебаний балки, закреплённой с двух сторон. В следующей таблице показаны результаты расчёта консольной балки.

Тип моды Собственная частота, Гц Собственная частота, Гц при ΔT = 10 K Отношение частот
Первая изгибная, в направлении оси z 8063.79 8066.92 1.00039
Первая изгибная, в направлении оси y 16049.1 16053.7 1.00028
Первая крутильная 132233 132265 1.00025
Первая продольная 250000 250050 1.0002

Изображение форм колебаний для консольной балки.
Формы мод колебаний консольной балки.

Первое, на что следует обратить внимание: для изгибных форм колебаний расчетные значения собственных частот достаточно сильно изменяются при учете термических деформаций. На первой моде изменение составляет 0.6%. Для всех других мод относительный сдвиг частоты значительно меньше. Но если балку сделать тоньше, разница будет ещё более заметной. Причины такого явления мы подробно разберём в следующих разделах.

Изучение эффекта снижения жёсткости под механическим напряжением

В случае балки, закреплённой с двух сторон, тепловое расширение приводит к сжимающему осевому (продольному) напряжению. С учётом исходных данных, напряжение составляет -10 MPa (в соответствии с формулой xΔT). Это напряжение приводит к значительному снижению жёсткости балки. Данный эффект часто называют упрочнением в напряженно-деформированном состоянии (stress stiffening), обычно он возникает в конструкциях при растягивающих напряжениях (tensile stresses). При сжимающих напряжениях жёсткость конструкции уменьшается, иногда такой эффект называют разупрочнением (stress softening).

Давайте взглянем на это с другой стороны и проведем линейный анализ устойчивости для балки. Для этого добавим в модель исследование Linear Buckling (Анализ потери устойчивости) и используем в качестве нагрузки тепловое расширение, вызванное разницей температур ΔT = 10 K. В результате не сложно получить, что коэффициент критической нагрузки равен 80.

Изображение первой моды потери устойчивости балки.
Первая мода потери устойчивости.

При линейном допущении получается, что увеличение температуры до 800 K приводит к потери устройчивости балки. При критической нагрузке, вызывающей потерю устойчивости, жёсткость достигает нуля. Полагая, что жёсткость линейно уменьшается со сжимающим напряжением, при разнице температур ΔT = 10 K, её относительное изменение можно оценить по формуле

1-\frac{1}{80} = 0.9875

Так как собственная частота пропорциональная квадратному корню из величины жёсткости, то можно оценить и её уменьшение \sqrt{0.9875}=0.9937, что хорошо соответствует расчётному значению 0.9943. Снижение жёсткости под напряжением также влияет на крутильные и продольные моды, но не так явно и заметно, как в случае с изгибными колебаниями.

Влияние изменения геометрических размеров на собственные частоты

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

Собственные частоты изгибных, крутильных и поперечных мод можно оценить по следующим формулам:

f_b \propto \frac{1}{L^2} \sqrt{\frac{EI}{\rho A}}
f_t \propto \frac{1}{L} \sqrt{\frac{GK}{\rho J}}
f_a \propto \frac{1}{L} \sqrt{\frac{E}{\rho}}

В уравнениях использовались следующие переменные:

  • I = Осевой момент инерции вокруг оси изгиба
  • G = Модуль сдвига
  • K = Модуль упругости при кручении
  • J = Полярный момент инерции вокруг оси балки

Предположим, что начальные размеры балки — L0 x a0 x b0, где a0 > b0. После теплового расширения размеры – L x a x b.

Расширения (деформации) в трёх ортогональных направлениях — εx, εy и εz, соответственно. В данном случае они линейно зависят от теплового расширения — εx = αxΔT, εy = αyΔT, and εz = αzΔT. В общем случае, это может быть любой тип неупругой деформации.

Тогда геометрические размеры масштабируются по следующим формулам:

\begin{align} L& =L_0(1+\epsilon_x)\\ A& = ab = a_0b_0(1+\epsilon_y)(1+\epsilon_z)\\ I_y& = \frac{ab^3}{12} = \frac{a_0b_0^3}{12}(1+\epsilon_y)(1+\epsilon_z)^3\\ I_z& = \frac{ba^3}{12} = \frac{b_0a_0^3}{12}(1+\epsilon_z)(1+\epsilon_y)^3\\ K& = \frac{ab^3}{3}F_1(a/b) \approx \frac{a_0b_0^3}{3}F_1(a_0/b_0)(1+\epsilon_y)(1+\epsilon_z)^3\\ J& =\frac{ab^3+ba^3}{12}=\frac{a_0b_0^3(1+\epsilon_y)(1+\epsilon_z)^3+b_0a_0^3(1+\epsilon_z)(1+\epsilon_y)^3}{12} \end{align}

Плотность материала также изменяется. Поскольку та же масса теперь ограничена большим объёмом, получаем:

\rho = \frac{\rho_0}{(1+\epsilon_x)(1+\epsilon_y)(1+\epsilon_z)}

Подставляя эти выражения в формулы для оценки собственных частот, можно рассчитать следующие ожидаемые сдвиги собственных частот:

\frac{f_{b,z}}{f_{b0,z}} = \sqrt{\frac{(1+\epsilon_y)(1+\epsilon_z)^3}{(1+\epsilon_x)^3}} \approx 1-\frac{3\epsilon_x}{2}+\frac{\epsilon_y}{2}+\frac{3\epsilon_z}{2}
\frac{f_{b,y}}{f_{b0,y}} = \sqrt{\frac{(1+\epsilon_z)(1+\epsilon_y)^3}{(1+\epsilon_x)^3}} \approx 1-\frac{3\epsilon_x}{2}+\frac{3\epsilon_y}{2}+\frac{\epsilon_z}{2}
\frac{f_a}{f_{a0}} = \sqrt{\frac{(1+\epsilon_y)(1+\epsilon_z)}{(1+\epsilon_x)}} \approx 1-\frac{\epsilon_x}{2}+\frac{\epsilon_y}{2}+\frac{\epsilon_z}{2}

Так как тепловые расширения очень малы, можно полагать, что приближённые разложения в ряды первого порядка будут вполне актуальными и точными.

При крутильных колебаниях ситуация несколько сложнее, так как в выражении для расчёта полярного момента J разные степени у переменных a и b. Но, если учесть, что a = 2b для данной геометрии, можно получить похожее выражение.

\frac{f_{t}}{f_{t0}} = \sqrt{\frac{5(1+\epsilon_z)(1+\epsilon_y)^3}{(1+\epsilon_x)((1+\epsilon_z)^2+ 4(1+\epsilon_y)^2))}} \approx 1-\frac{\epsilon_x}{2}-\frac{3\epsilon_y}{10}+\frac{6\epsilon_z}{5}

Теперь, давайте сравним расчётные значения сдвига частот с выведенными аналитическими оценками для консольной балки. Результаты приведены в таблице ниже, и соответствие между ними достаточно высокое.

Тип моды Расчетное значение Аналитическая оценка
Первая изгибная, в направлении оси z 1.00039 1.00040
Первая изгибная, в направлении оси y 1.00028 1.00030
Первая крутильная 1.00025 1.00025
Первая продольная 1.00020 1.00020

Эффект о задания ограничений

Граничные условия полной фиксации на концах балки приводят к локальным концентрациям напряжений при увеличении температуры, поскольку поперечные смещения ограничены.

График осевого напряжения, вызванного увеличением температуры, в балке, закреплённой с двух сторон.
Осевое напряжение, вызванное увеличением температуры на 10 K, в балке, закреплённой с двух сторон.

Это может привести к двум эффектам:

  1. Упрочнение в напряженно-деформированном состоянии может возникнуть в элементе, в котором, как ожидается, будут только объёмные изменения
  2. Размер поперечного сечения больше не будет постоянным из-за ограниченного поперечного смещения (как в вышеприведённом примере)

Чтобы определить, к каким эффектам могут привести те или иные ограничения, вы должны полагаться на свой инженерный опыт. Как правило, сам рассматриваемый компонент и его окружающие объекты зависят от изменений температуры. В данном случае в COMSOL Multiphysics у нас есть хорошая возможность добавить тепловое расширение к границам с ограничениями в модели. Давайте посмотрим, как это повлияет на решение.

Скриншот окна настроек подузла Thermal Expansion для граничного условия Fixed Constraint.
Тепловое расширение учитывается в граничных условиях Fixed Constraint (Фиксированное ограничение) для балки, закреплённой с двух сторон.

Для консольной балки результаты изменились таким образом, что теперь полностью совпадают с аналитическими оценками.

Тип моды Стандартные фиксации Ограничения со снятыми внутренними напряжениями (Stress-Free) Аналитическая оценка
Первая изгибная, в направлении оси z 1.00039 1.00040 1.00040
Первая изгибная, в направлении оси y 1.00028 1.00030 1.00030
Первая крутильная 1.00025 1.00026 1.00025
Первая продольная 1.00020 1.00020 1.00020

Учет в модели температурных зависимостей для свойств материала

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

В данном онлайн справочнике вы найдёте зависимости модуля упругости от температуры для различных металлов. Их жёсткость уменьшается с ростом температуры. Для многих материалов относительное изменение жёсткости составляет порядка 10-4 K-1. Это значит, что при изменении температуры на 10 K, можно ожидать относительное изменение жёсткости материала порядка 0.1%. Данный эффект может быть гораздо больше, чем эффект, вызванный изменениями геометрических размеров, который мы рассчитывали выше.

Небольшое предупреждение: При измерениях температурных зависимостей Модуля Юнга важно знать, учитывались ли геометрические изменения, вызванные тепловым расширением. Другими словами, вы должны знать, рассчитывался ли модуль Юнга относительно исходных размеров или размеров после нагрева.

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

Коэффициент теплового расширения обычно увеличивается с ростом температуры. Относительная чувствительность обычно порядка 10-3 K-1. Казалось бы это довольно много. Однако, если посмотреть, как этот коэффициент записывается в уравнениях, становится понятно, что он почти не влияет на результат.

Для большинства материалов, которые собраны в модуле расширения Библиотека материалов (Material Library) в COMSOL Multiphysics, заданы температурные зависимости свойств. В данном примере давайте вручную добавим линейную зависимость температуры от модуля Юнга. Для этого выполним следующие операции:

  1. В настройках группы Basic (Основные свойства) выберете Temperature (Температура) в поле Model Inputs (Входные данные модели)
  2. Нажмите Add (Добавить), чтобы увидеть имя используемой переменной – T
  3. Запишите выражение, в котором задается зависимость модуля Юнга от температуры T

Вы также можете задать функцию, а затем вызвать её в модели, указав температуру T в качестве аргумента.

Скриншот окна с добавлением линейной температурной зависимости для свойств материала.
Добавление линейной зависимости свойства материала от температуры.

В настройках узла Linear Elastic Material (Линейный упругий материал) теперь станет активной вкладка Model Input (Входные данные модели). В ней вы можете задать температуру, которая будет использоваться для материала.

Скриншот, иллюстрирующий добавление температуры в материале с использованием поля Model Input.
Указание температуры материала во вкладке Model Input.

С учётом уменьшения модуля Юнга на 1·10-4 K-1, результирующий сдвиг частот становится отрицательным, а не положительным, как при постоянном модуле Юнга (показано в таблице ниже).

Тип моды Ограничения со снятыми внутренними напряжениями при постоянном E Ограничения со снятыми внутренними напряжениями при учете зависимости E от температуры Разность
Первая изгибная, в направлении оси z 1.00040 0.99990 -0.00050
Первая изгибная, в направлении оси y 1.00030 0.99980 -0.00050
Первая крутильная 1.00026 0.99976 -0.00050
Первая продольная 1.00020 0.99970 -0.00050

Сдвиг оказался таким же, как и предполагалось для всех мод – модуль Юнга уменьшился в 1·10-3 раза, а собственные частоты уменьшились пропорционально его квадратному корню. Фактически, вы можете добавить учет изменения модуля Юнга в линеаризованные уравнения для аналитической оценки сдвига частот:

\frac{f_{b,z}}{f_{b0,z}} \approx 1+ \left (-\frac{3\alpha_x}{2}+\frac{\alpha_y}{2}+\frac{3\alpha_z}{2}+\frac{\beta}{2} \right)\Delta T
\frac{f_{b,y}}{f_{b0,y}} \approx 1 + \left (-\frac{3\alpha_x}{2}+\frac{3 \alpha_y}{2}+\frac{\alpha_z}{2}+\frac{\beta}{2} \right)\Delta T
\frac{f_a}{f_{a0}} \approx 1 + \left (-\frac{\alpha_x}{2}+\frac{\alpha_y}{2}+\frac{\alpha_z}{2}+\frac{\beta}{2} \right)\Delta T
\frac{f_{t}}{f_{t0}} \approx 1 + \left (-\frac{3\alpha_x}{2}-\frac{3 \alpha_y}{10}+\frac{6\alpha_z}{5}+\frac{\beta}{2} \right)\Delta T

Предполагается, что E=E_0(1+\beta \Delta T). Значение коэффициента β обычно отрицательно, в данном случае β = -10-4 K-1.

Для общего случая изотропного теплового расширения эти выражения упрощаются до:

\frac{f_{b,z}}{f_{b0,z}} = \frac{f_{b,y}}{f_{b0,y}} = \frac{f_a}{f_{a0}} \approx 1+ \left (\frac{\alpha}{2}+\frac{\beta}{2} \right)\Delta T
\frac{f_{t}}{f_{t0}} \approx 1 + \left (-\frac{3\alpha}{5}+\frac{\beta}{2} \right)\Delta T

Комментарии по численной точности

В рассматриваемой задаче мы ищем изменения частот с точностью в одну миллионную часть (ppm). Это значит, что необходимо избегать ложных ошибок при округлениях. В этом разделе статьи речь пойдёт о том, как можно обеспечить оптимальную точность.

В настройках узла Eigenfrequency (Исследование на собственные частоты) в поле Search for eigenfrequencies around (Поиск ближайший собственных частот) запишите корректное значение с нужным порядком величины.

Скриншот окна с обновлёнными настройками исследования Eigenfrequency.
Обновлёние настроек исследования Eigenfrequency.

Также следует уменьшить Relative tolerance (Относительная погрешность) в настройках узла Eigenvalue Solver (Решатель на собственные значения).

Скриншот окна с уменьшенным значением Relative tolerance в настройках узла Eigenvalue Solver.
Уменьшение значения параметра Relative tolerance в настройках узла Eigenvalue Solver.

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

Если у вас есть основания полагать, что задача стала плохо обусловленной, как это может быть в случае с тонкими конструкциями, в настройках прямого решателя (Direct) активируйте настройку Iterative refinement (Итеррационное сгущение).

Скриншот настроек прямого (Direct) решателя.
Окно настроек прямого решателя с активированной опцией Iterative refinement.

Комментарии по используемой методике обработки неупругих деформаций в условиях геометрической нелинейности

Начиная с версии 5.3 в COMSOL Multiphysics® был изменён метод обработки неупругих деформаций в условиях геометрической нелинейности. Мультипликативное разложение тензора деформаций теперь применяется по умолчанию, вместо вычитания деформаций, которое использовалось ранее. Именно в следствие этих изменений теперь можно выполнять рассматриваемый в данной статье тип анализа с очень высокой точностью. Давайте рассмотрим (несколько искусственный) случай, где превышение температуры составляет 3·104 K и у материалов нет температурно-зависимых свойств. Это означает, что растяжения

1+\epsilon_x = 1.3
1+\epsilon_y = 1.6
1+\epsilon_z = 1.9

приводят к изменению объёма в 3.952 раза.

Тогда можно сравнить результаты расчета на собственные частоты в случае предварительного напряжения со стандартным расчётом собственных частот большой балки с параметрами L = 13 мм, a = 1.6 мм, b = 0.95 мм и плотностью в 3.952 раза меньше – ρ = 253.036 кг/м3. Естественно, это приведёт к большому увеличению собственных частот, поскольку нагреваемый объект намного больше, но плотность его меньше. Относительные изменения частоты для двух случаев показаны в таблице ниже.

Тип моды Тепловое расширение и расчет на собственные частоты с учетом преднапряжения Большая геометрия и меньшая плотность
Первая изгибная, в направлении оси z 2.2309 2.2308
Первая изгибная, в направлении оси y 1.8759 1.8759
Первая крутильная 1.6702 1.6695
Первая продольная 1.5292 1.5292

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

Резюме

В данной статье мы показали, как точно проводить расчет сдвига собственных частот, вызванного температурными изменениями с использованием COMSOL Multiphysics. Также мы обсудили учет и вклад таких эффектов как: снижение жёсткости под механическим напряжением, геометрические деформации и температурные зависимости свойств материала.

Дополнительные материалы


Комментарии (0)

Оставить комментарий
Войти | Регистрация
Загрузка...
РУБРИКАТОР БЛОГА COMSOL
РУБРИКИ
ТЕГИ