Оцените поток решения PDE
[___] = evaluateCGradient( возвращает поток решения УЧП для стационарного уравнения в 2-D или 3-D точках, заданных в results,querypoints)querypoints.
[___] = evaluateCGradient(___, возвращает поток решения системы PDE для индексов уравнений (компонентов) iU)iU. При оценке потока для системы PDE задайте iU после входных параметров в любом из предыдущих синтаксисов.
Первая размерность cgradx, cgrady, и, в 3-D случае, cgradz соответствует точкам запроса. Второе измерение соответствует индексам уравнений iU.
[___] = evaluateCGradient(___, возвращает поток решения УЧП для зависящего от времени уравнения или системы зависящих от времени уравнений во времени iT)iT. При оценке потока для зависящего от времени УЧП задайте iT после входных параметров в любом из предыдущих синтаксисов. Для системы зависящих от времени PDE задайте оба индексов уравнения (компоненты) iU и временные индексы iT.
Первая размерность cgradx, cgrady, и, в 3-D случае, cgradz соответствует точкам запроса. Для одного зависящего от времени УЧП второе измерение соответствует временным шагам iT. Для системы зависящих от времени PDE, второе измерение соответствует индексам уравнений iU, и третья размерность соответствует временным шагам iT.
[ возвращает поток УЧП решения 2-D задачи в узловых точках треугольного mesh. Форма выхода массивов, cgradx,cgrady]
= evaluateCGradient(results)cgradx и cgrady, зависит от количества PDE, для которых results является решением. Первая размерность cgradx и cgrady представляет индексы узлов. Для системы стационарных или зависящих от времени PDE, второе измерение представляет индексы уравнений. Для одного зависящего от времени УЧП, второе измерение представляет временные шаги. Третья размерность представляет индексы временного шага для системы зависящих от времени PDE.
[ возвращает поток УЧП решения 3-D задачи в узловых точках четырехгранного mesh. Первая размерность cgradx,cgrady,cgradz]
= evaluateCGradient(results)cgradx, cgrady, и cgradz представляет индексы узлов. Второе измерение представляет индексы уравнений. Для системы стационарных или зависящих от времени PDE, второе измерение представляет индексы уравнений. Для одного зависящего от времени УЧП, второе измерение представляет временные шаги. Третья размерность представляет индексы временного шага для системы зависящих от времени PDE.
Решите задачу на L-образной мембране с нулевыми граничными условиями Дирихле. Оцените тензорный продукт c- коэффициент и градиенты решения скаляра эллиптической задачи в узловых и произвольных местоположениях. Постройте график результатов.
Создайте модель PDE и геометрию для этой задачи.
model = createpde; geometryFromEdges(model,@lshapeg); pdegplot(model,'FaceLabels','on')
Задайте граничные условия и коэффициенты.
applyBoundaryCondition(model,'dirichlet','Edge',1:model.Geometry.NumEdges,'u',0); specifyCoefficients(model,'m',0,'d',0,'c',10,'a',0,'f',1,'Face',1); specifyCoefficients(model,'m',0,'d',0,'c',5,'a',0,'f',1,'Face',2); specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1,'Face',3);
Создайте сетку геометрии и решите проблему.
generateMesh(model,'Hmax',0.05);
results = solvepde(model);
u = results.NodalSolution;Вычислите поток решения и постройте график результатов.
[cgradx,cgrady] = evaluateCGradient(results); figure pdeplot(model,'XYData',u,'Contour','on','FlowData',[cgradx,cgrady])
Вычислите поток решения на сетке от -1 до 1 в каждом направлении с помощью матрицы точек запроса.
v = linspace(-1,1,37); [X,Y] = meshgrid(v); querypoints = [X(:),Y(:)]'; [cgradxq,cgradyq] = evaluateCGradient(results,querypoints);
Кроме того, можно задать точки запроса следующим X,Y вместо того, чтобы задавать их как матрицу.
[cgradxq,cgradyq] = evaluateCGradient(results,X,Y);
Постройте график результата с помощью quiver функция построения графика.
figure quiver(X(:),Y(:),cgradxq,cgradyq) xlabel('x') ylabel('y')
Вычислите напряжения в консольной балке, подверженной сдвигу, загрузке на свободном конце.
Создайте модель PDE и геометрию для этой задачи.
N = 3; model = createpde(N); importGeometry(model,'SquareBeam.STL'); pdegplot(model,'FaceLabels','on')
Задайте коэффициенты и примените граничные условия.
E = 2.1e11; nu = 0.3; c = elasticityC3D(E, nu); a = 0; f = [0;0;0]; specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a','f',f); applyBoundaryCondition(model,'dirichlet','Face',6,'u',[0 0 0]); applyBoundaryCondition(model,'neumann','Face',5,'g',[0,0,-3e3]);
Создайте сетку геометрии и решите проблему.
generateMesh(model,'Hmax',25,'GeometricOrder','quadratic'); results = solvepde(model);
Вычислите напряжение, то есть продукт c- коэффициент и градиенты перемещения.
[sig_xx,sig_yy,sig_zz] = evaluateCGradient(results);
Постройте график нормального компонента напряжения вдоль x-направление. Верхний фрагмент балки испытывает натяжение, а нижний фрагмент - сжатие.
figure
pdeplot3D(model,'ColorMapData',sig_xx(:,1))
Задайте линию поперек балки снизу вверх в середине и середине ширины. Вычислите напряжения вдоль линии.
zg = linspace(0, 100, 10); xg = 250*ones(size(zg)); yg = 50*ones(size(zg)); [sig_xx,sig_xy,sig_xz] = evaluateCGradient(results,xg,yg,zg,1);
Постройте график нормального напряжения вдоль x-направление.
figure plot(sig_xx,zg) grid on xlabel('\sigma_{xx}') ylabel('z')
Вычислите напряжения в идеализированной 3-D механической части под приложенной нагрузкой. Во-первых, создайте модель PDE для этой задачи.
N = 3; model = createpde(N);
Импортируйте геометрию и постройте график.
importGeometry(model,'BracketWithHole.stl'); figure pdegplot(model,'FaceLabels','on') view(30,30) title('Bracket with Face Labels')
figure pdegplot(model,'FaceLabels','on') view(-134,-32) title('Bracket with Face Labels, Rear View')
Задайте коэффициенты и примените граничные условия.
E = 200e9; % elastic modulus of steel in Pascals nu = 0.3; % Poisson's ratio c = elasticityC3D(E,nu); a = 0; f = [0;0;0]; % Assume all body forces are zero specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a,'f',f); applyBoundaryCondition(model,'dirichlet','Face',4,'u',[0,0,0]); distributedLoad = 1e4; % Applied load in Pascals applyBoundaryCondition(model,'neumann','Face',8,'g',[0,0,-distributedLoad]);
Создайте сетку геометрии и решите проблему.
bracketThickness = 1e-2; % Thickness of horizontal plate with hole, meters hmax = bracketThickness; % Maximum element length for a moderately fine mesh generateMesh(model,'Hmax',hmax,'GeometricOrder','quadratic'); result = solvepde(model);
Создайте сетку. Для этой сетки вычислите тензор напряжения, который является продуктом c- коэффициент и градиенты перемещения.
v = linspace(0,0.2,21); [xq,yq,zq] = meshgrid(v); [cgradx,cgrady,cgradz] = evaluateCGradient(result);
Извлечение отдельных компонентов напряжений.
sxx = cgradx(:,1); sxy = cgradx(:,2); sxz = cgradx(:,3); syx = cgrady(:,1); syy = cgrady(:,2); syz = cgrady(:,3); szx = cgradz(:,1); szy = cgradz(:,2); szz = cgradz(:,3);
Вычислите напряжение фон Мизеса.
sVonMises = sqrt( 0.5*( (sxx-syy).^2 + (syy -szz).^2 +...
(szz-sxx).^2) + 3*(sxy.^2 + syz.^2 + szx.^2));Постройте график напряжения фон Мизеса. Максимальное напряжение возникает на самом слабом участке. Этот раздел имеет наименьший материал для поддержки приложенной нагрузки.
pdeplot3D(model,'colormapdata',sVonMises)
Решите 2-D задачу переходной теплопередачи в квадратной области и вычислите тепловой поток через конвективный контур.
Создайте модель PDE для этой задачи.
numberOfPDE = 1; model = createpde(numberOfPDE);
Создайте геометрию.
g = @squareg; geometryFromEdges(model,g); pdegplot(model,'EdgeLabels','on') xlim([-1.2,1.2]) ylim([-1.2,1.2]) axis equal
Задайте свойства материала и условия окружающей среды.
rho = 7800; cp = 500; k = 100; Text = 25; hext = 5000;
Задайте коэффициенты. Примените изолированные граничные условия к три ребер и условие свободного контура конвекции к правому ребру.
specifyCoefficients(model,'m',0,'d',rho*cp,'c',k,'a',0,'f',0); applyBoundaryCondition(model,'neumann','Edge', [1,3,4],'q',0,'g',0); applyBoundaryCondition(model,'neumann','Edge', 2,'q',hext,'g',Text*hext);
Установите начальные условия: равномерную комнатную температуру в области и более высокую температуру на левом крае.
setInitialConditions(model,25);
setInitialConditions(model, 100, 'Edge', 4);Сгенерируйте mesh и решите задачу, используя 0:1000:200000 как вектор времени.
generateMesh(model); tlist = 0:1000:200000; results = solvepde(model,tlist);
Задайте линию на контуре конвекции, чтобы вычислить тепловой поток через нее.
yg = -1:0.1:1; xg = ones(size(yg));
Оцените продукт коэффициента c и пространственных градиентов в (xg,yg).
[qx,qy] = evaluateCGradient(results,xg,yg,1:length(tlist));
Пространственно интегрируйте градиенты, чтобы получить тепловой поток для каждого временного шага.
HeatFlowX(1:length(tlist)) = -trapz(yg,qx(:,1:length(tlist)));
Постройте конвективный тепловой поток с течением времени.
figure plot(tlist,HeatFlowX) title('Heat flow across convection boundary') xlabel('Time') ylabel('Heat flow')
Решите задачу теплопередачи для следующей геометрии 2-D, состоящей из квадрата и алмаза, изготовленного из различных материалов. Вычислите плотность теплового потока и постройте график как векторное поле.
Создайте модель PDE для этой задачи.
numberOfPDE = 1; model = createpde(numberOfPDE);
Создайте геометрию, которая состоит из квадрата с встраиваемым алмазом.
SQ1 = [3; 4; 0; 3; 3; 0; 0; 0; 3; 3]; D1 = [2; 4; 0.5; 1.5; 2.5; 1.5; 1.5; 0.5; 1.5; 2.5]; gd = [SQ1,D1]; sf = 'SQ1+D1'; ns = char('SQ1','D1'); ns = ns'; dl = decsg(gd,sf,ns); geometryFromEdges(model,dl); pdegplot(model,'EdgeLabels','on','FaceLabels','on') xlim([-1.5,4.5]) ylim([-0.5,3.5]) axis equal
Установите параметры для квадратной области.
rho_sq = 2; C_sq = 0.1; k_sq = 10; Q_sq = 0; h_sq = 0;
Установите параметры для алмазной области.
rho_d = 1; C_d = 0.1; k_d = 2; Q_d = 4; h_d = 0;
Задайте коэффициенты для обоих поддоменов. Примените граничные и начальные условия.
specifyCoefficients(model,'m',0,'d',rho_sq*C_sq,'c',k_sq,'a',h_sq,'f',Q_sq,'Face',1); specifyCoefficients(model,'m',0,'d',rho_d*C_d,'c',k_d,'a',h_d,'f',Q_d,'Face',2); applyBoundaryCondition(model,'dirichlet','Edge',[1,2,7,8],'h',1,'r',0); setInitialConditions(model,0);
Создайте сетку геометрии и решите проблему. Чтобы захватить самую динамическую часть процесса распределения тепла, решите задачу с помощью logspace(-2,-1,10) как вектор времени.
generateMesh(model); tlist = logspace(-2,-1,10); results = solvepde(model,tlist); u = results.NodalSolution;
Вычислите плотность теплового потока. Постройте график решения с помощью изотермических линий с помощью контурного графика и постройте график векторного поля теплового потока со стрелами. Направление теплового потока (от более высоких до более низких температур) противоположно направлению . Поэтому используйте -cgradx и -cgrady чтобы показать тепловой поток.
[cgradx,cgrady] = evaluateCGradient(results); figure pdeplot(model,'XYData',u(:,10),'Contour','on','FlowData',[-cgradx(:,10),-cgrady(:,10)],'ColorMap','hot')
В то время как results объект содержит решение и его градиент (оба рассчитаны в узловых точках треугольного или четырехгранного mesh), он не содержит потока решения PDE. Чтобы вычислить поток в узловых местоположениях, вызовите evaluateCGradient без указания местоположений. По умолчанию, evaluateCGradient использует узловые локации.
evaluateGradient | interpolateSolution | PDEModel | StationaryResults | TimeDependentResults