Exponenta Banner

assembleFEMatrices

Собрать матрицы конечных элементов

Свернуть все на странице

Синтаксис

FEM = assembleFEMatrices(model)
FEM = assembleFEMatrices(model,matrices)
FEM = assembleFEMatrices(model,bcmethod)
FEM = assembleFEMatrices(___,state)

Описание

пример

FEM = assembleFEMatrices(model) возвращает структурный массив, содержащий все матрицы конечных элементов для УЧП задачи, заданной как model.

пример

FEM = assembleFEMatrices(model,matrices) возвращает структурный массив, содержащий только указанные матрицы конечных элементов.

пример

FEM = assembleFEMatrices(model,bcmethod) собирает матрицы конечных элементов и накладывает граничные условия с помощью метода, заданного bcmethod.

пример

FEM = assembleFEMatrices(___,state) собирает матрицы конечных элементов с помощью входного времени или решения, заданного в state массив структур. Функция использует time поле структуры для зависящих от времени моделей и поле решения u для нелинейных моделей. Можно использовать этот аргумент с любым из предыдущих синтаксисов.

Примеры

свернуть все

Открыть Live Script

Создайте модель PDE для уравнения Пуассона на L-образной мембране с нулевыми граничными условиями Дирихле.

model = createpde(1);
geometryFromEdges(model,@lshapeg);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1);
applyBoundaryCondition(model,'Edge',1:model.Geometry.NumEdges,'u',0);

Сгенерируйте mesh и получите матрицы конечных элементов по умолчанию для задачи и mesh.

generateMesh(model,'Hmax',0.2);
FEM = assembleFEMatrices(model)
FEM = struct with fields:
    K: [401x401 double]
    A: [401x401 double]
    F: [401x1 double]
    Q: [401x401 double]
    G: [401x1 double]
    H: [80x401 double]
    R: [80x1 double]
    M: [401x401 double]
    T: [401x401 double]
Открыть Live Script

Сделайте расчеты быстрее, задав, какие матрицы конечных элементов собирать.

Создайте переходную тепловую модель и включите геометрию встроенной функции squareg.

thermalmodel = createpde('thermal','steadystate');
geometryFromEdges(thermalmodel,@squareg);

Постройте график геометрии с метками ребер.

pdegplot(thermalmodel,'EdgeLabels','on')
xlim([-1.1 1.1])
ylim([-1.1 1.1])

Figure contains an axes. The axes contains 5 objects of type line, text.

Задайте теплопроводность материала и внутреннего источника тепла.

thermalProperties(thermalmodel,'ThermalConductivity',0.2);
internalHeatSource(thermalmodel,10);

Установите граничные условия.

thermalBC(thermalmodel,'Edge',[1,3],'Temperature',100);

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

generateMesh(thermalmodel);

Соберите жесткость и большие матрицы.

FEM_KM = assembleFEMatrices(thermalmodel,'KM')
FEM_KM = struct with fields:
    K: [1541x1541 double]
    M: [1541x1541 double]

Теперь соберите матрицы конечных элементов M, K, A и F.

FEM_MKAF = assembleFEMatrices(thermalmodel,'MKAF')
FEM_MKAF = struct with fields:
    M: [1541x1541 double]
    K: [1541x1541 double]
    A: [1541x1541 double]
    F: [1541x1 double]

Четыре матрицы M, K, A и F соответствуют дискретизированным версиям коэффициентов PDE m, c, a и f. Эти четыре матрицы также представляют область конечноэлементной модели PDE. Вместо того, чтобы указывать их явным образом, можно использовать domain аргумент.

FEMd = assembleFEMatrices(thermalmodel,'domain')
FEMd = struct with fields:
    M: [1541x1541 double]
    K: [1541x1541 double]
    A: [1541x1541 double]
    F: [1541x1 double]

Четыре матрицы Q, G, H и R соответствуют дискретизированным версиям q, g, h и r в спецификации граничных условий Неймана и Дирихле. Эти четыре матрицы также представляют контур конечноэлементной модели PDE. Используйте boundary аргумент для сборки только этих матриц.

FEMb = assembleFEMatrices(thermalmodel,'boundary')
FEMb = struct with fields:
    H: [74x1541 double]
    R: [74x1 double]
    G: [1541x1 double]
    Q: [1541x1541 double]
Открыть Live Script

Создайте модель PDE для уравнения Пуассона на L-образной мембране с нулевыми граничными условиями Дирихле. 

model = createpde(1);
geometryFromEdges(model,@lshapeg);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1);
applyBoundaryCondition(model,'Edge',1:model.Geometry.NumEdges,'u',0);

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

generateMesh(model,'Hmax',0.2);

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

FEMn = assembleFEMatrices(model,'nullspace')
FEMn = struct with fields:
    Kc: [321x321 double]
    Fc: [321x1 double]
     B: [401x321 double]
    ud: [401x1 double]
     M: [321x321 double]

Получите решение для УЧП с помощью nullspace матрицы конечных элементов.

un = FEMn.B*(FEMn.Kc\FEMn.Fc) + FEMn.ud;

Сравните этот результат с решением, заданным solvepde. Эти два решения идентичны.

u1 = solvepde(model);
norm(un - u1.NodalSolution)
ans = 0

Получите матрицы конечных элементов после наложения граничного условия с помощью жесткого пружинного подхода. Этот подход сохраняет степени свободы Дирихле, но налагает на них большое наказание.

FEMs = assembleFEMatrices(model,'stiff-spring')
FEMs = struct with fields:
    Ks: [401x401 double]
    Fs: [401x1 double]
     M: [401x401 double]

Получите решение УЧП с помощью жестко-пружинных конечных элементарных матриц. Этот метод дает менее точное решение.

us = FEMs.Ks\FEMs.Fs;
norm(us - u1.NodalSolution)
ans = 0.0098
Открыть Live Script

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

Создайте переходную несущую модель для решения твердотельной (3-D) задачи.

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

Создайте геометрию и включите ее в модель. Постройте график геометрии.

gm = multicylinder(0.01,0.05);
addVertex(gm,'Coordinates',[0,0,0.05]);
structuralmodel.Geometry = gm;
pdegplot(structuralmodel,'FaceLabels','on','FaceAlpha',0.5)

Figure contains an axes. The axes contains 3 objects of type quiver, patch, line.

Задайте модуль Юнга и коэффициент Пуассона.

structuralProperties(structuralmodel,'Cell',1,'YoungsModulus',201E9, ...
                                              'PoissonsRatio',0.3, ...
                                              'MassDensity',7800);

Задайте, что нижняя часть цилиндра является фиксированным контуром.

structuralBC(structuralmodel,'Face',1,'Constraint','fixed');

Задайте гармоническое давление в верхней части цилиндра.

structuralBoundaryLoad(structuralmodel,'Face',2,'Pressure',5E7,'Frequency',50);

Задайте нулевое начальное перемещение и скорость.

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

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

generateMesh(structuralmodel,'GeometricOrder','linear');
tlist = linspace(0,1,300);

Соберите матрицы конечных элементов для начального временного шага.

state.time = tlist(1);
FEM_domain = assembleFEMatrices(structuralmodel,state)
FEM_domain = struct with fields:
    K: [6609x6609 double]
    A: [6609x6609 double]
    F: [6609x1 double]
    Q: [6609x6609 double]
    G: [6609x1 double]
    H: [252x6609 double]
    R: [252x1 double]
    M: [6609x6609 double]

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

state.time = tlist(1);
FEM_boundary_init = assembleFEMatrices(structuralmodel,'G',state)
FEM_boundary_init = struct with fields:
    G: [6609x1 double]


state.time = tlist(floor(length(tlist)/2));
FEM_boundary_half = assembleFEMatrices(structuralmodel,'G',state)
FEM_boundary_half = struct with fields:
    G: [6609x1 double]


state.time = tlist(end);
FEM_boundary_final = assembleFEMatrices(structuralmodel,'G',state)
FEM_boundary_final = struct with fields:
    G: [6609x1 double]
Открыть Live Script

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

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

thermalmodelS = createpde('thermal','steadystate');

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

r1 = [3 4 -.5 .5 .5 -.5  -.8 -.8 .8 .8];
r2 = [3 4 -.05 .05 .05 -.05  -.4 -.4 .4 .4];
gdm = [r1; r2]';

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

g = decsg(gdm,'R1-R2',['R1'; 'R2']');

Преобразуйте decsg формат в объект геометрии. Включите геометрию в модель и постройте график геометрии.

geometryFromEdges(thermalmodelS,g);
figure
pdegplot(thermalmodelS,'EdgeLabels','on'); 
axis([-.9 .9 -.9 .9]);

Figure contains an axes. The axes contains 9 objects of type line, text.

Установите температуру на левый край равной 100 степеней. Установите тепловой поток из блока на правом ребре равным -10. Верхние и нижние ребра и ребра внутри полости все изолированы: нет теплопередачи через эти ребра.

thermalBC(thermalmodelS,'Edge',6,'Temperature',100);
thermalBC(thermalmodelS,'Edge',1,'HeatFlux',-10);

Задайте теплопроводность материала как простую линейную функцию от температуры u.

k = @(~,state) 0.7+0.003*state.u;
thermalProperties(thermalmodelS,'ThermalConductivity',k);

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

generateMesh(thermalmodelS);

Вычислите установившееся решение.

Rnonlin = solve(thermalmodelS);

Поскольку теплопроводность нелинейна (зависит от температуры), вычислите системные матрицы, соответствующие сходящейся температуре. Присвойте распределение температуры u поле state массив структур. Потому что u поле должно содержать вектор-строку, транспонировать распределение температуры.

state.u = Rnonlin.Temperature.';

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

FEM = assembleFEMatrices(thermalmodelS,'nullspace',state)
FEM = struct with fields:
    Kc: [1277x1277 double]
    Fc: [1277x1 double]
     B: [1320x1277 double]
    ud: [1320x1 double]
     M: [1277x1277 double]

Вычислите решение с помощью системных матриц, чтобы убедиться, что они дают ту же температуру, что и Rnonlin.

u = FEM.B*(FEM.Kc\FEM.Fc) + FEM.ud;

Сравните этот результат с решением, заданным solve.

norm(u - Rnonlin.Temperature)
ans = 9.8994e-05

Входные параметры

свернуть все

Объект модели, заданный как PDEModel объект, ThermalModel объект, StructuralModel объект, или ElectroMagneticModel объект.

Пример: model = createpde(1)

Пример: thermalmodel = createpde('thermal','steadystate')

Пример: structuralmodel = createpde('structural','static-solid')

Пример: emagmodel = createpde('electromagnetic','electrostatic')

Метод включения граничных условий, заданный как 'none', 'nullspace', или 'stiff-spring'. Для получения дополнительной информации см. «Алгоритмы».

Пример: FEM = assembleFEMatrices(model,'nullspace')

Типы данных: char | string

Матрицы для сборки, заданные как:

  • Матричные идентификаторы, такие как 'F', 'MKF', 'K', и так далее - Собрать соответствующие матрицы. Каждая заглавная буква представляет одну матрицу: K, A, F, Q, G, H, R, M, и T. Можно объединить несколько букв в один вектор символов или строку, например 'MKF'.

  • 'boundary' - Собрать все матрицы, относящиеся к контурам геометрии.

  • 'domain' - Собрать все связанные с областью матрицы.

Пример: FEM = assembleFEMatrices(model,'KAF')

Типы данных: char | string

Время для зависящих от времени моделей и решение для нелинейных моделей, заданное в массиве структур. Поля массива представляют следующие значения:

  • state.time содержит неотрицательное число, задающее значение времени для зависящих от времени моделей.

  • state.u содержит матрицу решений размера N -by - Np, которая может использоваться для сборки матриц в нелинейной настройке задачи, где коэффициенты являются функциями state.u. Здесь N количество уравнений в системе, а Np число узлов в mesh.

Пример: state.time = tlist(end); FEM = assembleFEMatrices(model,'boundary',state)

Выходные аргументы

свернуть все

Матрицы конечных элементов, возвращенные как структурный массив. Используйте bcmethod и matrices аргументы для определения конечных матриц элемента, которые необходимо собрать.

Поля в структурном массиве зависят от bcmethod:

  • Если значение 'none', тогда поля будут K, A, F, Q, G, H, R, и M.

  • Если значение 'nullspace', тогда поля будут Kc, Fc, B, ud, и M.

  • Если значение 'stiff-spring', тогда поля будут Ks, Fs, и M.

Поля в структурном массиве также зависят от matrices:

  • Если значение boundary, тогда все поля являются матрицами, связанными с контурами геометрии.

  • Если значение domain, тогда все поля являются связанными с областью матрицами.

  • Если значение является матричным идентификатором или идентификаторами, такими как 'F', 'MKF', 'K'и так далее, тогда поля являются соответствующими матрицами.

Для получения дополнительной информации см. «Алгоритмы».

Алгоритмы

свернуть все

Partial Differential Equation Toolbox™ решает уравнения вида

m∂2u∂t2+d∂u∂t−∇·(c⊗∇u)+au=f

и собственные уравнения вида

−∇·(c⊗∇u)+au=λduили−∇·(c⊗∇u)+au=λ2mu

при условиях контура Дирихле hu = r и контур Неймана, n · (c⊗∇u)+qu=g.

assembleFEMatrices возвращает следующие полные матрицы и векторы конечных элементов, которые представляют соответствующую задачу УЧП:

  • K - матрица жесткости, интеграл дискретизированной версии c коэффициент.

  • M - большая матрица, интеграл дискретизированной версии m или d коэффициенты. M является ненулевым для зависящих от времени и собственных значений задач.

  • A является неотъемлемой частью дискретизированной версии a коэффициент.

  • F является неотъемлемой частью дискретизированной версии f коэффициент. Для тепловых, электромагнитных и структурных проблем, F является источником или вектором нагрузки тела.

  • Q является неотъемлемой частью дискретизированной версии q термин в граничном условии Неймана.

  • G является неотъемлемой частью дискретизированной версии g термин в граничном условии Неймана. Для структурных проблем, G является краевым вектором нагрузки.

  • The H и R матрицы происходят непосредственно из условий Дирихле и mesh.

Навязывание условий Контура дирихлет

The 'nullspace' метод устраняет условия Дирихле из задачи с помощью линейной алгебры подхода. Он генерирует объединенные конечноэлементные матрицы Kc, Fc, B, и векторные ud соответствует уменьшенной системе Kc*u = Fc, где Kc = B'*(K + A + Q)*B, и Fc = B'*(F + G). The B матрица охватывает ядро пространства из столбцов H (матрица условия Дирихле, представляющая h*ud = r). The R вектор представляет условия Дирихле в H*ud = R. The ud вектор имеет размер вектора решения. Его элементы являются нулями повсюду, кроме мест Дирихле степеней свободы (DoFs), где они содержат предписанные значения.

Из 'nullspace' матрицы, можно вычислить решение u как

u = B*(Kc\Fc) + ud.

Если вы собрали конкретный набор матриц, например G и M, можно наложить граничные условия на G и M следующим образом. Сначала вычислите nullspace столбцов H.

[B,Or] = pdenullorth(H);
ud = Or*((H*Or\R)); % Vector with known value of the constraint DoF.

Затем используйте B матрица следующим образом. Чтобы исключить степени свободы Дирихле из вектора нагрузки G, использовать:

GwithBC = B'*G

Чтобы исключить степени свободы Дирихле из большой матрицы, используйте:

M = B'*M*B

Можно исключить степени свободы Дирихле из других векторов и матриц с помощью того же метода.

The 'stiff-spring' метод преобразует граничные условия Дирихле в граничные условия Неймана с помощью жестко-пружинного приближения. Оно возвращает матрицу Ks и вектор Fs который вместе представляет другой тип комбинированных матриц конечных элементов. Приблизительное решение u = Ks\Fs. По сравнению с 'nullspace' метод, 'stiff-spring' метод генерирует матрицы быстрее, но обычно дает менее точные решения.

Примечание

Внутренне тулбокс использует 'nullspace' подход к наложению граничных условий Дирихле при вычислении решения с помощью solvepde и solve.

Степени свободы (DoFs)

Если число узлов в модели NumNodes, и количество уравнений N, затем длина векторов-столбцов u и ud является N*NumNodes. Тулбокс присваивает идентификаторы степеням свободы в u и ud:

  • Значения от 1 до NumNodes соответствуют первому уравнению.

  • Записи из NumNodes+1 на 2*NumNodes соответствуют второму уравнению.

  • Записи из 2*NumNodes+1 на 3*NumNodes соответствуют третьему уравнению.

Тот же подход применяется ко всем остальным записям, вплоть до N*NumNodes.

Для примера в 3-D структурной модели длина вектора решения u является 3*NumNodes. Первый NumNodes записи соответствуют x-размещение в каждом узле, следующая NumNodes записи соответствуют y-размещение, и следующая NumNodes записи соответствуют z-размещение.

Тепловой, структурный и электромагнитный анализ

В тепловом анализе m и a коэффициенты являются нулями. Теплопроводность преобразуется в c коэффициент. Продукта массовой плотности и удельного тепла сопоставляются с d коэффициент. Внутренний источник тепла преобразуется в f коэффициент. Температура на контуре соответствует члену граничного условия Дирихле r с h = 1. Различные формы граничного теплового потока, такие как сам тепловой поток, излучательная способность и коэффициент конвекции, сопоставляются с терминами граничного условия Неймана q и g.

В структурном анализе a коэффициент равен нулю. Модуль Юнга и отношение Пуассона к c коэффициент. Массовая плотность преобразуется в m коэффициент. Карта нагрузок тела на f коэффициент. Смещения, ограничения и компоненты перемещения вдоль осей сопоставляются с условиями граничного условия Дирихле h и r. Краевые нагрузки, такие как давление, поверхностные тяги и поступательные жесткости, соответствуют терминам граничных условий Неймана q и g. Когда вы задаете модель демпфирования при помощи параметров Релея Alpha и Betaдискретизированная матрица демпфирования C вычисляется с помощью большой матрицы M и матрица жесткости K как C = Alpha*M + Beta*K.

В электростатическом и магнитостатическом анализах m, a, и d коэффициенты являются нулями. Карта относительной диэлектрической проницаемости и относительной проницаемости по c коэффициент. Плотность заряда и плотность тока сопоставлены с f коэффициент. Напряжение и магнитный потенциал на контуре соответствуют члену граничного условия Дирихле r с h = 1.

См. также

| | | | |

Темы

Введенный в R2016a

Partial Differential Equation Toolbox производными

Поддержка