При знании решения в терминах интерфейса 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.

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

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

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

При знании идентификаторов степени свободы для данных идентификаторов узла используйте 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))';

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

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 и модального числа степеней свободы, можно восстановить решение для полной модели. 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