Динамический анализ зафиксированного луча

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

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

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

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

Уравнения

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

Лагранжевая формулировка уравнений равновесия может быть записана

$ρ u_{}^{¨} - \nabla \cdot (F \cdot S) = f$

где $ρ$ существенная плотность, $u$ вектор смещения, $F$ градиент деформации, $S$ второй тензор напряжения Пайола-Кирчофф, и $f$ вектор массовой силы. Это уравнение может также быть написано в следующей форме тензора:

$ρ {u_{}^{¨}}_{i} - \frac{\partial}{\partial x_{j}} ((\frac{\partial u_{i}}{\partial x_{k}} + δ_{i k}) S_{k j}) = f_{i}$

Несмотря на то, что большие отклонения составляются в этой формулировке, она принята, что деформации остаются маленькими так, чтобы линейные эластичные конститутивные отношения применялись. Кроме того, материал принят, чтобы быть изотропным. Для 2D плоского случая напряжения конститутивные отношения могут быть написаны в матричной форме:

${\begin{array}{c} S_{11} \\ S_{22} \\ S_{12} \end{array}} = [\begin{array}{ccc} C_{11} & C_{12} \\ C_{12} & C_{22} \\ 2 G_{12} \end{array}] {\begin{array}{c} E_{11} \\ E_{22} \\ E_{12} \end{array}}$

где $E_{i j}$ Зелено-лагранжев тензор деформации, заданный как

$E_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}} + \frac{\partial u_{k}}{\partial x_{i}} \frac{\partial u_{k}}{\partial x_{j}})$

Для изотропного материала

$C_{11} = C_{22} = \frac{E}{1 - ν^{2}}$

$C_{12} = \frac{E ν}{1 - ν^{2}}$

$G_{12} = \frac{E}{2 (1 + ν)}$

где $E$ модуль Янга и $ν$ отношение Пуассона.

Читатели, которым интересно более подробно о лагранжевой формулировке для нелинейной эластичности, могут консультироваться, например, со ссылкой 1.

Уравнения, представленные выше полностью, описывают геометрически нелинейную плоскую задачу напряжения. Работа, требуемая преобразовывать их в форму, приемлемую для Тулбокса УЧП, значительно упрощена при помощи Symbolic Math Toolbox MATLAB. Symbolic Math Toolbox может выполнить необходимые алгебраические манипуляции и вывести функцию MATLAB, задающую c-коэффициент, который может быть передан функциям Тулбокса УЧП. Эта функция, cCoefficientLagrangePlaneStress, показан в приложении ниже.

Setup задач

Создайте модель PDE.

N = 2; % Two PDE in plane stress elasticity
model = createpde(N);

Задайте геометрию. Во-первых, задайте длину и толщину луча,

blength = 5; % Beam length, in.
height = .1; % Thickness of the beam, in.

Рисунок зафиксированного луча показан на рисунке ниже.

Поскольку геометрия луча и загрузка симметричны о центре луча (x = blength/2), модель может быть упрощена путем рассмотрения только правильной половины луча.

l2 = blength/2;
h2 = height/2;

Создайте ребра прямоугольника, представляющего луч с этими двумя операторами.

rect = [3 4 0 l2 l2 0 -h2 -h2  h2 h2]';
g = decsg(rect,'R1',('R1')');

geometryFromEdges функция создает геометрический объект из ребер и хранит его в модели.

pg = geometryFromEdges(model,g);

Постройте геометрию и отобразите метки ребра. Метки ребра необходимы для идентификации ребра при применении граничных условий.

figure
pdegplot(g,'EdgeLabels','on');
title('Geometry With Edge Labels Displayed');

Масштабируйте график, таким образом, метки являются просматриваемыми.

axis([-.1 1.1*l2 -5*h2 5*h2]);

Выведите коэффициенты уравнения с помощью свойств материала. Для линейного случая c-матрица-коэффициентов является постоянной.

E = 3.0e7; % Young's modulus of the material, lbs/in^2
gnu = .3; % Poisson's ratio of the material
rho = .3/386; % Density of the material
G = E/(2.*(1+gnu));
mu = 2*G*gnu/(1-gnu);
c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
f = [0 0]'; % No body forces
specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0, 'f', f);

Примените граничные условия. Условием симметрии является x- смещение равный нуль на левом крае.

symBC = applyBoundaryCondition(model,'mixed','Edge',4,'u',0,'EquationIndex',1);

x- и y- смещения равняются нулю вдоль правого края.

clampedBC = applyBoundaryCondition(model,'dirichlet','Edge',2,'u',[0 0]);

Примените постоянное вертикальное напряжение вдоль верхнего края.

sigma = 2e2;
presBC = applyBoundaryCondition(model,'neumann','Edge',3,'g',[0 sigma]);

Установите нулевые начальные смещения и скорости.

setInitialConditions(model,0,0);

Сгенерируйте mesh приблизительно с восемью элементами через толщину луча.

generateMesh(model,'Hmax',height/8,'MesherVersion','R2013a');

Линейное решение

Настройте аналитический промежуток. Здесь, tfinal итоговое время в анализе.

tfinal = 3e-3;
tlist = linspace(0,tfinal,100);

Вычислите зависящее от времени решение.

result = solvepde(model,tlist);

Интерполируйте решение в центре геометрии y- компонент (2 компонента) в любом случае.

xc = 1.25;
yc = 0;
u4Linear = interpolateSolution(result,xc,yc,2,1:length(tlist));

Нелинейное решение

Задайте коэффициенты для нелинейного случая. cCoefficientLagrangePlaneStress функционируйте берет изотропные свойства материала и местоположение и структуры состояния, и возвращает c-матрицу для нелинейного плоского расчета напряжений. Приняты маленькие деформации; т.е. E и $ν$ независимы от решения. Тулбокс УЧП вызывает пользовательские коэффициентные функции с местоположением аргументов и состоянием. Функциональный cCoefficientLagrangePlaneStress ожидает аргументы E, гну, местоположение, состояние. c задан ниже как анонимная функция, чтобы обеспечить интерфейс между этими двумя функциональными подписями. (Функциональный cCoefficientLagrangePlaneStress может использоваться с любым геометрическим нелинейным плоским расчетом напряжений модели, сделанной из изотропного материала.)

c  = @(location, state) cCoefficientLagrangePlaneStress(E,gnu,location,state);

specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0 , 'f', f);

Вычислите зависящее от времени решение.

result = solvepde(model,tlist);

Как прежде, интерполируйте решение в центре геометрии y-компонента (2 компонента) в любом случае.

u4NonLinear = interpolateSolution(result,xc,yc,2,1:length(tlist));

Постройте решения

Фигура ниже показов y-отклонение в центре луча как функция времени. Нелинейный анализ вычисляет смещения, которые существенно меньше линейного анализа. Это "напряжение, укрепляющее" эффект, также отражается в более высокой частоте колебания от нелинейного анализа.

figure
plot(tlist,u4Linear(:),tlist,u4NonLinear(:));
legend('Linear','Nonlinear');
title 'Deflection at Beam Center'
xlabel 'Time, seconds'
ylabel 'Deflection, inches'
grid on

Ссылки

Малверн, Лоуренс Э., введение в механику непрерывного носителя, Prentice Hall, 1969.

Приложение - нелинейная C-коэффициентная-функция

Функциональный cCoefficientLagrangePlaneStress вычисляет c-матрицу-коэффициентов для большой формулировки напряжения плоскости функции Лагранжа смещения.

type cCoefficientLagrangePlaneStress

function c = cCoefficientLagrangePlaneStress(E, nu, loc, state)
%cCoefficientLagrangePlaneStress Calculate c-coefficient for nonlinear plane stress
% Calculate the c-coefficient for a geometrically nonlinear Lagrangian formulation 
% of plane stress elasticity. The strain measure is the Green-Lagrange strain
% tensor. The stress is the second Piola-Kirchoff stress tensor. The material
% is assumed to be isotropic with linear behavior (Hooke's law applies).
% 
% E  - Young's modulus of the linear isotropic material
% nu - Poisson's ratio for the material
% p  - matrix of point (node) locations
% t  - element connectivity matrix
% u  - current displacement vector

%    This function was generated by the Symbolic Math Toolbox version 6.0.
%    31-Jan-2014 09:50:09

% Copyright 2014-2015 The MathWorks, Inc.


ux = reshape(state.ux,2,[]);
uy = reshape(state.uy,2,[]);



dudx=ux(1,:); 
dvdx=ux(2,:);
dudy=uy(1,:);  
dvdy=uy(2,:);
 
% if(~isempty(u))
%   [ux,uy] = pdegrad(p,t,u);
%   dudx=ux(1,:); dudy=uy(1,:); dvdx=ux(2,:); dvdy=uy(2,:);
% else
%   dudx = zeros(1, size(t,2)); dudy=dudx; dvdx=dudx; dvdy=dudx;
% end

t4 = 1/(nu^2-1);
t6 = 1/(1+nu);
t7 = E*dudy.*t4*.25;
t8 = dudx+1.0;
t9 = E*dudy.*t4.*t8*.25;
t10 = dvdy+1.0;
t11 = t7+t9-E*dvdx.*t6.*t10*.25;
t12 = dvdy.*2.0;
t13 = dudx.^2;
t14 = dudy.^2;
t15 = dvdy.^2;
t16 = dvdx.^2;
t17 = E*dvdx.*t4.*(1.0./2.0);
t18 = E*dudx.*dvdx.*t4*.25;
t19 = t17+t18-E*dudy.*t6*.25-E*dudy.*dvdy.*t6.*(1.0./8.0);
t20 = E*dudy.*dvdx.*nu.*t4*.25;
t21 = t20-E*t6.*(1.0./2.0)-E*dudx.*t6*.25-E*dvdy.*t6*.25-E*dudx.*dvdy.*t6.*(1.0./8.0);
t22 = dudx.*2.0;
t23 = dvdy+2.0;
t24 = nu-1.0;
t25 = E*nu.*t4;
t26 = E*dudx.*nu.*t4.*(1.0./2.0);
t27 = E*dvdy.*nu.*t4.*(1.0./2.0);
t28 = E*dudx.*dvdy.*nu.*t4*.25;
t29 = t25+t26+t27+t28-E*dudy.*dvdx.*t6.*(1.0./8.0);
t30 = E*dudy.*t4.*t23.*(1.0./8.0);
t31 = E*dudy.*dvdy.*t4.*(1.0./8.0);
t32 = t7+t30+t31-E*dvdx.*t6.*(1.0./8.0)-E*dvdx.*t4.*t8.*t24.*(1.0./8.0);
t33 = dudy.*2.0;
t34 = dvdx.*2.0;
t35 = dudx.*dudy.*2.0;
t36 = dvdx.*dvdy;
t37 = t33+t34+t35+t36;
t38 = 1.0./t24;
t39 = E*dvdx.*t23.*t38.*(1.0./8.0);
t40 = t39-E*t6.*t37.*(1.0./8.0);
out1 = [E*t4.*(dudx.*6.0+t13.*2.0+t14+t16+4.0)*.25+E*nu.*t4.*(t12+t15)*.25;
  t11;
  t19;
  t29;
  t11;
  E*t4.*(t12+t13+t14.*2.0+t15+t22+2.0)*.25+E*nu.*t4.*(t16-2.0)*.25;
  t21;
  t32;
  t19;
  t21;
  E*t4.*(t12+t13+t15+t16.*2.0+t22+2.0)*.25+E*nu.*t4.*(t14-2.0)*.25;
  t40;
  t29;
  t32;
  t40;
  E*t4.*(dvdy.*6.0+t14+t15.*2.0+t16+4.0)*.25+E*nu.*t4.*(t13+t22)*.25];

c = -out1([1 5 6 9 10 13 14 11 15 16], :);

end

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