В этом примере показано, как устранить степени свободы (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]); figure pdegplot(modelT,'VertexLabels','on','FaceAlpha',0.5) view([78 2.5])
Создайте сетку.
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);
Определите интерфейсы виража, используя фиксированные и загруженные границы. В этом случае модель уменьшенного порядка сохраняет степени свободы (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, которая определена в конце этого примера.
sol = ode15s(@(t,y) CMSODEf(t,y,Kode,Fode,loadODEDoF),tlist,u0,odeoptions);
Вычислите значения переменной ОДУ и производных времени.
[displ,vel] = deval(sol,tlist);
Постройте график z-смещения в нагруженной вершине и сравните его с третьим DoF в решении уменьшенной системы ОДУ.
figure plot(tlist,RT.Displacement.uz(loadedVertex,:)) hold on plot(tlist,displ(3,:),'r*') title('Z-Displacement at Loaded Vertex') legend('full model','rom')
Зная решение с точки зрения интерфейсов 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]; 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