В этом примере показано, как устранить степени свободы (число степеней свободы), которое не находится на контурах интереса при помощи техники моделирования уменьшаемого порядка Craig-Бамптона. Пример также использует меньший суперэлемент размерности, чтобы анализировать динамику системы. Для сравнения пример также выполняет прямой анализ переходных процессов исходной структуры.
Создайте структурную модель для анализа переходных процессов.
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]); figure pdegplot(modelT,'VertexLabels','on','FaceAlpha',0.5) view([78 2.5])
Сгенерируйте mesh.
generateMesh(modelT);
Прикладывайте синусоидальную сконцентрированную силу в z-направлении на новой вершине.
structuralBoundaryLoad(modelT,'Vertex',loadedVertex, ... 'Force',[0;0;10],'Frequency',6000);
Задайте нулевые начальные условия.
structuralIC(modelT,'Velocity',[0 0 0],'Displacement',[0 0 0]);
Решите модель.
tlist = 0:0.00005:3E-3; RT = solve(modelT,tlist);
Задайте интерфейсы суперэлемента с помощью фиксированных и загруженных контуров. В этом случае уменьшаемая модель порядка сохраняет степени свободы (число степеней свободы) на подкрашенном и загруженной вершине при сжатии всего другого числа степеней свободы в пользу модального числа степеней свободы. Для лучшей производительности используйте набор ребер, ограничивающих поверхность 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, который задан в конце этого примера.
sol = ode15s(@(t,y) CMSODEf(t,y,Kode,Fode,loadODEDoF),tlist,u0,odeoptions);
Вычислите значения переменной ODE и производных времени.
[displ,vel] = deval(sol,tlist);
Постройте z-смещение в загруженной вершине и сравните его с третьей степенью свободы в решении уменьшаемой системы ОДУ.
figure plot(tlist,RT.Displacement.uz(loadedVertex,:)) hold on plot(tlist,displ(3,:),'r*') title('Z-Displacement at Loaded Vertex') legend('full model','rom')
Зная решение в терминах интерфейса 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]; iDispRT = interpolateDisplacement(RT, coordCenter); iDispRTrom = interpolateDisplacement(RTrom, coordCenter); figure plot(tlist,iDispRT.uz,'k') hold on plot(tlist,iDispRTrom.uz,'g*') title('Z-Displacement at Geometric Center') legend('full model','rom')
function f = CMSODEf(t,u,Kode,Fode,loadedVertex) Fode(loadedVertex) = 10*sin(6000*t); f = -Kode*u +Fode; end