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

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

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]);

Создайте сетку.

generateMesh(modelT);

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

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

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

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

Определите интерфейсы виража, используя фиксированные и загруженные границы. В этом случае модель с уменьшенным порядком сохраняет DoF на неподвижной грани и загруженной вершине при конденсации всех других DoF в пользу модальных DoF. Для повышения производительности используйте набор ребер, ограничивающих грань 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 к загруженным и фиксированным DoF с использованием данных, доступных в rom.

DoF соответствуют поступательным перемещениям. Если количество точек сетки в модели равно Nn, то панель инструментов присваивает идентификаторы DoF следующим образом: 1 кому Nn - перемещения по оси X, Nn+1 кому 2*Nn - смещения по оси y, и 2Nn+1 кому 3*Nn - перемещения по оси Z. Уменьшенный объект модели rom содержит эти идентификаторы для сохраненных DoF в 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обобщенные модальные DoF появляются после сохраненных DoF.

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

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

odeDoFs = [loadNodeROMIds;fixedIntModeIds];

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

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

Теперь у вас есть система ODE второго порядка. Использовать 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-направления, которое является третьим ДоФ в системе ОДУ. Учет линеаризации для получения системы первого порядка дает загруженному ОДУ DoF.

loadODEDoF = numODE + 3;

Укажите массовую матрицу и якобиан для решателя ОДУ.

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

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

u0 = zeros(2*numODE,1);

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

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

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

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

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

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