evaluateCGradient
Оцените поток решения для УЧП
Синтаксис
Описание
[___] = evaluateCGradient( возвращает поток решения для УЧП для стационарного уравнения в 2D или 3-D точках, заданных в results,querypoints)querypoints.
[___] = evaluateCGradient(___, возвращает поток решения системы УЧП для индексов уравнения (компоненты) iU)iU. При оценке потока для системы УЧП задайте iU после входных параметров в любом из предыдущих синтаксисов.
Первая размерность cgradx, cgrady, и, в 3-D случае, cgradz соответствует точкам запроса. Второе измерение соответствует индексам уравнения iU.
[___] = evaluateCGradient(___, возвращает поток решения для УЧП для зависящего от времени уравнения или системы зависящих от времени уравнений во времена iT)iT. При оценке потока для зависящего от времени УЧП задайте iT после входных параметров в любом из предыдущих синтаксисов. Для системы зависящих от времени УЧП задайте оба индекса уравнения (компоненты) iU и индексы времени iT.
Первая размерность cgradx, cgrady, и, в 3-D случае, cgradz соответствует точкам запроса. Для одного зависящего от времени УЧП второе измерение соответствует тактам iT. Для системы зависящих от времени УЧП второе измерение соответствует индексам уравнения iU, и третья размерность соответствует тактам iT.
[ возвращает поток решения для УЧП 2D проблемы в узлах треугольной mesh. Форма выходных массивов, cgradx,cgrady]
= evaluateCGradient(results)cgradx и cgrady, зависит от количества УЧП для который results решение. Первая размерность cgradx и cgrady представляет индексы узла. Для системы стационарных или зависящих от времени УЧП второе измерение представляет индексы уравнения. Для одного зависящего от времени УЧП второе измерение представляет такты. Третья размерность представляет индексы такта для системы зависящих от времени УЧП.
[ возвращает поток решения для УЧП 3-D проблемы в узлах четырехгранной mesh. Первая размерность cgradx,cgrady,cgradz]
= evaluateCGradient(results)cgradx, cgrady, и cgradz представляет индексы узла. Второе измерение представляет индексы уравнения. Для системы стационарных или зависящих от времени УЧП второе измерение представляет индексы уравнения. Для одного зависящего от времени УЧП второе измерение представляет такты. Третья размерность представляет индексы такта для системы зависящих от времени УЧП.
Примеры
Скалярная эллиптическая проблема
Решите задачу на 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)
Проблема теплопередачи на квадрате
Решите 2D переходную задачу теплопередачи на квадратной области и вычислите тепловой поток через конвективный контур.
Создайте модель 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')
Теплопередача между двумя квадратами, сделанными из различных материалов
Решите задачу теплопередачи для следующей 2D геометрии, состоящей из квадрата и ромба, сделанного из различных материалов. Вычислите плотность потока тепла и постройте ее как векторное поле.
Создайте модель 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), это не содержит поток решения для УЧП. Чтобы вычислить поток в узловых местоположениях, вызватьevaluateCGradientне задавая местоположения. По умолчанию,evaluateCGradientиспользует узловые местоположения.
Смотрите также
evaluateGradient | interpolateSolution | PDEModel | StationaryResults | TimeDependentResults