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 решения для нелинейных моделей. Можно использовать этот аргумент с любым из предыдущих синтаксисов.

Примеры

свернуть все

Скрипт Open Live Script

Создайте модель PDE для уравнения Poisson на 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]
Скрипт Open 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 соответствуют дискретизированным версиям коэффициентов УЧП m, c, a, и f. Эти четыре матрицы также представляют область модели конечного элемента УЧП. Вместо того, чтобы задать их явным образом, можно использовать 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 в спецификации граничного условия Неймана и Дирихле. Эти четыре матрицы также представляют контур модели конечного элемента УЧП. Используйте boundary аргумент, чтобы собрать только эти матрицы.

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

Создайте модель PDE для уравнения Poisson на 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
Скрипт Open 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]
Скрипт Open Live Script

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

Создайте установившуюся тепловую модель.

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

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

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, и так далее — Собирают соответствующие матрицы. Каждая прописная буква представляет одну матрицу: KAFQGHRM, и 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', затем полями является KAFQGHR, и M.

  • Если значением является 'nullspace', затем полями является KcФК 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 граничный вектор загрузки.

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

Наложение граничных условий Дирихле

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

От '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

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

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

Примечание

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

Степени свободы (число степеней свободы)

Если количеством узлов в модели является 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. Граничные загрузки, такие как давление, поверхностные тяги, и поступательный stiffnesses, соответствуют Неймановым условиям граничного условия q и g. Когда вы задаете модель затухания при помощи Рейли, ослабляющего параметры Alpha и Beta, дискретизированный ослабляющий матричный C вычисляется при помощи большой матрицы M и матрица жесткости K как C = Alpha*M + Beta*K.

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

Смотрите также

| | | | |

Темы

Введенный в R2016a

Документация Partial Differential Equation Toolbox

Поддержка