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