Зная решение с точки зрения интерфейса DoFs и модальных Число степеней свободы, восстановите решение для полной модели.

Создайте несущую модель для временного анализа.

modelT = createpde('structural','transient-solid');

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

gm = multicuboid(0.05,0.003,0.003);
modelT.Geometry = gm;

Постройте график геометрии, отобразив метки граней и ребер.

figure
pdegplot(modelT,'FaceLabels','on','FaceAlpha',0.5)
view([71 4])

figure
pdegplot(modelT,'EdgeLabels','on','FaceAlpha',0.5)
view([71 4])

Задайте модуль Юнга, отношение Пуассона и массовую плотность материала.

structuralProperties(modelT,'YoungsModulus',210E9, ...
                            'PoissonsRatio',0.3, ...
                            'MassDensity',7800);

Закрепите один конец балки.

structuralBC(modelT,'Edge',[2 8 11 12],'Constraint','fixed');

Добавить вершину в центре грани 3.

loadedVertex = addVertex(gm,'Coordinates',[0.025 0.0 0.0015]);

Сгенерируйте mesh.

generateMesh(modelT);

Применить синусоидальную концентрированную силу в z-направлении на новой вершине.

structuralBoundaryLoad(modelT,'Vertex',loadedVertex, ...
                              'Force',[0;0;10],'Frequency',6000);

Задайте нулевые начальные условия.

structuralIC(modelT,'Velocity',[0 0 0],'Displacement',[0 0 0]);

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

structuralSEInterface(modelT,'Edge',[2 8 11 12]);
structuralSEInterface(modelT,'Vertex',loadedVertex);

Уменьшите структуру, сохраняя все фиксированные режимы интерфейса до 5e5.

rom = reduce(modelT,'FrequencyRange',[-0.1,5e5]);

Затем используйте модель пониженного порядка, чтобы симулировать динамику переходного процесса. Используйте ode15s функция непосредственно для интегрирования ОДУ редуцированной системы. Работа с уменьшенной моделью требует индексации в уменьшенные системные матрицы rom.K и rom.M. Сначала создайте отображения индексов K и M к загруженным и фиксированным Число степеней свободы при помощи данных, доступных в rom.

Число степеней свободы соответствуют поступательным перемещениям. Если количество точек сетки в модели Nn, затем тулбокс присваивает идентификаторы Число степеней свободы следующим образом: первый 1 на Nn являются перемещениями X, Nn+1 на 2*Nn являются y-перемещениями и 2Nn+1 на 3*Nn являются z-перемещениями. Уменьшенный объект модели rom содержит эти идентификаторы для сохранённых чисел степеней свободы в rom.RetainedDoF.

Создайте функцию, которая возвращает идентификаторы DoF, заданные идентификаторами узла, и число узлов.

getDoF = @(x,numNodes) [x(:); x(:) + numNodes; x(:) + 2*numNodes];

Зная идентификаторы DoF для заданных идентификаторов узла, используйте intersect функция для поиска необходимых индексов.

numNodes = size(rom.Mesh.Nodes,2);
 
loadedNode = findNodes(rom.Mesh,'region','Vertex',loadedVertex);
loadDoFs = getDoF(loadedNode,numNodes);
[~,loadNodeROMIds,~] = intersect(rom.RetainedDoF,loadDoFs);

В редуцированных матрицах rom.K и rom.Mобобщенные модальные Число степеней свободы появляются после сохраненных Число степеней свободы.

fixedIntModeIds = (numel(rom.RetainedDoF) + 1:size(rom.K,1))';

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

odeDoFs = [loadNodeROMIds;fixedIntModeIds];

Соответствующие компоненты rom.K и rom.M для интегрирования во времени:

Kconstrained = rom.K(odeDoFs,odeDoFs);
Mconstrained = rom.M(odeDoFs,odeDoFs);
numODE = numel(odeDoFs);

Теперь у вас есть система второго порядка ОДУ. Как использовать ode15s, преобразуйте это в систему ОДУ первого порядка путем применения линеаризации. Такая система первого порядка в два раза превышает размер системы второго порядка.

Mode = [eye(numODE,numODE),   zeros(numODE,numODE); ...
        zeros(numODE,numODE), Mconstrained];
Kode = [zeros(numODE,numODE), -eye(numODE,numODE); ...
        Kconstrained,         zeros(numODE,numODE)];
Fode = zeros(2*numODE,1);

Заданная концентрированная силовая нагрузка в полной системе - по z-направлению, которое является третьим степенью свободы в системе ОДУ. Учет линеаризации для получения системы первого порядка дает загруженную ОДУ Степень свободы.

loadODEDoF = numODE + 3;

Задайте большую матрицу и якобиан для решателя ОДУ.

odeoptions = odeset;
odeoptions = odeset(odeoptions,'Jacobian',-Kode);
odeoptions = odeset(odeoptions,'Mass',Mode);

Задайте нулевые начальные условия.

u0 = zeros(2*numODE,1);

Решить сокращенную систему можно с помощью ode15s и вспомогательной функции CMSODEf, который определен в конце этого примера.

tlist = 0:0.00005:3E-3;
sol = ode15s(@(t,y) CMSODEf(t,y,Kode,Fode,loadODEDoF),tlist,u0,odeoptions);

Вычислите значения переменной ODE и производных по времени.

[displ,vel] = deval(sol,tlist);

Зная решение с точки зрения интерфейса DoFs и модальных Число степеней свободы, можно восстановить решение для полной модели. The reconstructSolution функция требует смещения, скорости и ускорения на всех числах степеней свободы в rom. Создайте полный вектор решения, включая нулевые значения в фиксированных числах степеней свободы.

u = zeros(size(rom.K,1),numel(tlist));
ut = zeros(size(rom.K,1),numel(tlist));
utt = zeros(size(rom.K,1),numel(tlist));
u(odeDoFs,:) = displ(1:numODE,:);
ut(odeDoFs,:) = vel(1:numODE,:);
utt(odeDoFs,:) = vel(numODE+1:2*numODE,:);

Создайте объект результатов переходного процесса, используя это решение.

RTrom = reconstructSolution(rom,u,ut,utt,tlist);

Вычислите перемещение внутри в центре балки с помощью восстановленного решения.

coordCenter = [0;0;0];
iDispRTrom = interpolateDisplacement(RTrom, coordCenter);
figure
plot(tlist,iDispRTrom.uz)
title('Z-Displacement at Geometric Center')

function f = CMSODEf(t,u,Kode,Fode,loadedVertex)
Fode(loadedVertex) = 10*sin(6000*t);
f = -Kode*u +Fode;
end