Решите 3-D стационарную тепловую проблему.

Создайте тепловую модель для этой проблемы.

thermalmodel = createpde('thermal');

Импорт и печать геометрии блока.

importGeometry(thermalmodel,'Block.stl'); 
pdegplot(thermalmodel,'FaceLabel','on','FaceAlpha',0.5)
axis equal

Назначение свойств материала.

thermalProperties(thermalmodel,'ThermalConductivity',80);

Применить постоянную температуру 100 ° C к левой стороне блока (грань 1) и постоянную температуру 300 ° C к правой стороне блока (грань 3). Все остальные грани по умолчанию изолированы.

thermalBC(thermalmodel,'Face',1,'Temperature',100);
thermalBC(thermalmodel,'Face',3,'Temperature',300);

Выполните сетку геометрии и решите проблему.

generateMesh(thermalmodel);
thermalresults = solve(thermalmodel)

thermalresults = 
  SteadyStateThermalResults with properties:

    Temperature: [12691x1 double]
     XGradients: [12691x1 double]
     YGradients: [12691x1 double]
     ZGradients: [12691x1 double]
           Mesh: [1x1 FEMesh]

Решатель находит температуры и градиенты температуры в узловых местах. Для доступа к этим значениям используйте thermalresults.Temperature, thermalresults.XGradientsи так далее. Например, постройте график температур в узловых местах.

pdeplot3D(thermalmodel,'ColorMapData',thermalresults.Temperature)

Решите 2-D нестационарную тепловую проблему.

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

thermalmodel = createpde('thermal','transient');

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

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(thermalmodel,dl);
pdegplot(thermalmodel,'EdgeLabels','on','FaceLabels','on')
xlim([-1.5 4.5])
ylim([-0.5 3.5])
axis equal

Для квадратной области присвойте следующие тепловые свойства:

thermalProperties(thermalmodel,'ThermalConductivity',10, ...
                               'MassDensity',2, ...
                               'SpecificHeat',0.1, ...
                               'Face',1);

Для алмазной области присвойте следующие тепловые свойства:

thermalProperties(thermalmodel,'ThermalConductivity',2, ...
                               'MassDensity',1, ...
                               'SpecificHeat',0.1, ...
                               'Face',2);

Предположим, что ромбовидная область является источником тепла с плотностью $4\mathrm{}{\mathrm{}}^{\mathrm{W/m2}}$.

internalHeatSource(thermalmodel,4,'Face',2);

Прикладывайте постоянную температуру $$\mathrm{}{}^{}\mathrm{0\circ C}$$ к сторонам квадратной пластины.

thermalBC(thermalmodel,'Temperature',0,'Edge',[1 2 7 8]);

Установите начальную температуру 0 ° C.

thermalIC(thermalmodel,0);

Создайте сетку.

generateMesh(thermalmodel);

Динамика для этой проблемы очень быстрая. Температура достигает стационарного состояния примерно за 0,1 секунды. Чтобы зафиксировать интересную часть динамики, установите время решения на logspace(-2,-1,10). Эта команда возвращает 10 логарифмически разнесенных значений времени решения от 0,01 до 0,1.

tlist = logspace(-2,-1,10);

Решите уравнение.

thermalresults = solve(thermalmodel,tlist)

thermalresults = 
  TransientThermalResults with properties:

      Temperature: [1481x10 double]
    SolutionTimes: [1x10 double]
       XGradients: [1481x10 double]
       YGradients: [1481x10 double]
       ZGradients: []
             Mesh: [1x1 FEMesh]

Постройте график решения с изотермическими линиями, используя контурный график.

T = thermalresults.Temperature;
pdeplot(thermalmodel,'XYData',T(:,10),'Contour','on','ColorMap','hot')

Решите статическую структурную модель, представляющую биметаллический кабель, находящийся под натяжением.

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

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

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

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

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

structuralProperties(structuralmodel,'Cell',1,'YoungsModulus',110E9, ...
                                              'PoissonsRatio',0.28);
structuralProperties(structuralmodel,'Cell',2,'YoungsModulus',210E9, ...
                                              'PoissonsRatio',0.3);

Укажите, что грани 1 и 4 являются фиксированными границами.

structuralBC(structuralmodel,'Face',[1,4],'Constraint','fixed');

Задайте тягу поверхности для граней 2 и 5.

structuralBoundaryLoad(structuralmodel,'Face',[2,5],'SurfaceTraction',[0;0;100]);

Создайте сетку и решите проблему.

generateMesh(structuralmodel);
structuralresults = solve(structuralmodel)

structuralresults = 
  StaticStructuralResults with properties:

      Displacement: [1x1 FEStruct]
            Strain: [1x1 FEStruct]
            Stress: [1x1 FEStruct]
    VonMisesStress: [22281x1 double]
              Mesh: [1x1 FEMesh]

Решатель находит значения смещения, напряжения, деформации и напряжения Мизеса в узловых расположениях. Для доступа к этим значениям используйте structuralresults.Displacement, structuralresults.Stressи так далее. Значения смещения, напряжения и деформации в узловых местоположениях возвращаются как FEStruct объекты со свойствами, представляющими их компоненты. Обратите внимание, что свойства FEStruct объекты доступны только для чтения.

structuralresults.Displacement

ans = 
  FEStruct with properties:

           ux: [22281x1 double]
           uy: [22281x1 double]
           uz: [22281x1 double]
    Magnitude: [22281x1 double]

structuralresults.Stress

ans = 
  FEStruct with properties:

    sxx: [22281x1 double]
    syy: [22281x1 double]
    szz: [22281x1 double]
    syz: [22281x1 double]
    sxz: [22281x1 double]
    sxy: [22281x1 double]

structuralresults.Strain

ans = 
  FEStruct with properties:

    exx: [22281x1 double]
    eyy: [22281x1 double]
    ezz: [22281x1 double]
    eyz: [22281x1 double]
    exz: [22281x1 double]
    exy: [22281x1 double]

Постройте график деформированной формы с z-компонентом нормального напряжения.

pdeplot3D(structuralmodel,'ColorMapData',structuralresults.Stress.szz, ...
                          'Deformation',structuralresults.Displacement)

Решение для переходного отклика тонкой 3-D пластины под гармонической нагрузкой в центре.

Создайте переходную динамическую модель для 3-D проблемы.

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

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

gm = multicuboid([5,0.05],[5,0.05],0.01);
structuralmodel.Geometry = gm;
pdegplot(structuralmodel,'FaceLabels','on','FaceAlpha',0.5)

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

figure
pdegplot(structuralmodel,'FaceLabels','on','FaceAlpha',0.25)
axis([-0.2 0.2 -0.2 0.2 -0.1 0.1])

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

structuralProperties(structuralmodel,'YoungsModulus',210E9,...
                                     'PoissonsRatio',0.3,...
                                     'MassDensity',7800);

Укажите, что все грани на периферии тонкой 3-D пластины являются фиксированными границами.

structuralBC(structuralmodel,'Constraint','fixed','Face',5:8);

Приложите синусоидальную нагрузку давления к небольшой грани в центре пластины.

structuralBoundaryLoad(structuralmodel,'Face',12,'Pressure',5E7,'Frequency',25);

Создайте сетку с линейными элементами.

generateMesh(structuralmodel,'GeometricOrder','linear','Hmax',0.2);

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

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

Решите модель.

tlist = linspace(0,1,300);
structuralresults = solve(structuralmodel,tlist)

structuralresults = 
  TransientStructuralResults with properties:

     Displacement: [1x1 FEStruct]
         Velocity: [1x1 FEStruct]
     Acceleration: [1x1 FEStruct]
    SolutionTimes: [1x300 double]
             Mesh: [1x1 FEMesh]

Решатель находит значения смещения, скорости и ускорения в узловых местоположениях. Для доступа к этим значениям используйте structuralresults.Displacement, structuralresults.Velocityи так далее. Значения смещения, скорости и ускорения возвращаются как FEStruct объекты со свойствами, представляющими их компоненты. Обратите внимание, что свойства FEStruct объекты доступны только для чтения.

structuralresults.Displacement

ans = 
  FEStruct with properties:

           ux: [1873x300 double]
           uy: [1873x300 double]
           uz: [1873x300 double]
    Magnitude: [1873x300 double]

structuralresults.Velocity

ans = 
  FEStruct with properties:

           vx: [1873x300 double]
           vy: [1873x300 double]
           vz: [1873x300 double]
    Magnitude: [1873x300 double]

structuralresults.Acceleration

ans = 
  FEStruct with properties:

           ax: [1873x300 double]
           ay: [1873x300 double]
           az: [1873x300 double]
    Magnitude: [1873x300 double]

Найдите фундаментальный (самый низкий) режим 2-D консольного луча, предполагая распространенность состояния плоского напряжения.

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

length = 5;
height = 0.1;
E = 3E7;
nu = 0.3;
rho = 0.3/386;

Создайте модель «плоскость-напряжение» модели, назначьте геометрию и создайте сетку.

structuralmodel = createpde('structural','modal-planestress');
gdm = [3;4;0;length;length;0;0;0;height;height];
g = decsg(gdm,'S1',('S1')');
geometryFromEdges(structuralmodel,g);

Определите максимальный размер элемента (пять элементов по толщине балки).

hmax = height/5;
msh=generateMesh(structuralmodel,'Hmax',hmax);

Задайте свойства несущих конструкций и граничные зависимости.

structuralProperties(structuralmodel,'YoungsModulus',E, ...
                                     'MassDensity',rho, ... 
                                     'PoissonsRatio',nu);
structuralBC(structuralmodel,'Edge',4,'Constraint','fixed');

Вычислите аналитическую основную частоту (Гц), используя теорию луча.

I = height^3/12;
analyticalOmega1 = 3.516*sqrt(E*I/(length^4*(rho*height)))/(2*pi)

analyticalOmega1 = 126.9498

Укажите диапазон частот, включающий аналитически вычисленную частоту, и решите модель.

modalresults = solve(structuralmodel,'FrequencyRange',[0,1e6])

modalresults = 
  ModalStructuralResults with properties:

    NaturalFrequencies: [32x1 double]
            ModeShapes: [1x1 FEStruct]
                  Mesh: [1x1 FEMesh]

Решатель находит собственные частоты и значения модального смещения в узловых местоположениях. Для доступа к этим значениям используйте modalresults.NaturalFrequencies и modalresults.ModeShapes.

modalresults.NaturalFrequencies/(2*pi)

ans = 32×1
105 ×

    0.0013
    0.0079
    0.0222
    0.0433
    0.0711
    0.0983
    0.1055
    0.1462
    0.1930
    0.2455
      ⋮

modalresults.ModeShapes

ans = 
  FEStruct with properties:

           ux: [6511x32 double]
           uy: [6511x32 double]
    Magnitude: [6511x32 double]

Постройте график y-составляющей решения для основной частоты.

pdeplot(structuralmodel,'XYData',modalresults.ModeShapes.uy(:,1))
title(['First Mode with Frequency ', ...
        num2str(modalresults.NaturalFrequencies(1)/(2*pi)),' Hz'])
axis equal

Выполните анализ частотной характеристики камертона.

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

model = createpde('structural','frequency-solid');

Импортируйте геометрию камертона.

importGeometry(model,'TuningFork.stl');

Задайте модуль Юнга, коэффициент Пуассона и плотность массы для моделирования поведения линейного эластичного материала. Укажите все физические свойства в согласованных единицах измерения.

structuralProperties(model,'YoungsModulus',210E9, ...
                           'PoissonsRatio',0.3, ...
                           'MassDensity',8000);

Определение граней для применения граничных зависимостей и нагрузок путем печати геометрии с метками граней.

figure('units','normalized','outerposition',[0 0 1 1])
pdegplot(model,'FaceLabels','on')
view(-50,15)
title 'Geometry with Face Labels'

Наложение достаточных граничных ограничений для предотвращения движения жесткого тела при приложенной нагрузке. Обычно вы держите камертон вручную или устанавливаете его на стол. Чтобы создать простую аппроксимацию этого граничного условия, зафиксируйте область вблизи пересечения зубцов и ручки (грани 21 и 22).

structuralBC(model,'Face',[21,22],'Constraint','fixed');

Укажите нагрузку давления на стойку (грань 11) в виде короткого прямоугольного импульса давления. В частотной области этот импульс давления представляет собой единичную нагрузку, равномерно распределенную по всем частотам.

structuralBoundaryLoad(model,'Face',11,'Pressure',1);
flist = linspace(0,4000,150);
mesh = generateMesh(model,'Hmax',0.005); 
R = solve(model,2*pi*flist)

R = 
  FrequencyStructuralResults with properties:

           Displacement: [1x1 FEStruct]
               Velocity: [1x1 FEStruct]
           Acceleration: [1x1 FEStruct]
    SolutionFrequencies: [1x150 double]
                   Mesh: [1x1 FEMesh]

Постройте график частоты колебаний наконечника стойки, который является поверхностью 12. Найдите узлы на грани наконечника и постройте график y-компонента смещения на частоте, используя один из этих узлов.

excitedTineTipNodes = findNodes(mesh,'region','Face',12);
tipDisp = R.Displacement.uy(excitedTineTipNodes(1),:);

figure 
plot(flist,abs(tipDisp)) 
xlabel('Frequency');
ylabel('|Y-Displacement|');

Найдите отклонение 3-D консольной балки под неравномерной тепловой нагрузкой. Укажите тепловую нагрузку на структурную модель, используя решение из нестационарного теплового анализа для той же геометрии и сетки.

thermalmodel = createpde('thermal','transient');

gm = multicuboid(0.5,0.1,0.05);
thermalmodel.Geometry = gm;
pdegplot(thermalmodel,'FaceLabels','on','FaceAlpha',0.5)

mesh = generateMesh(thermalmodel);

thermalProperties(thermalmodel,'ThermalConductivity',5e-3, ...
                               'MassDensity',2.7*10^(-6), ...
                               'SpecificHeat',10);

thermalBC(thermalmodel,'Face',3,'Temperature',100);
thermalBC(thermalmodel,'Face',5,'Temperature',0);

internalHeatSource(thermalmodel,10);

thermalIC(thermalmodel,0);

tlist = [0:1e-4:2e-4];
thermalresults = solve(thermalmodel,tlist)

thermalresults = 
  TransientThermalResults with properties:

      Temperature: [3870x3 double]
    SolutionTimes: [0 1.0000e-04 2.0000e-04]
       XGradients: [3870x3 double]
       YGradients: [3870x3 double]
       ZGradients: [3870x3 double]
             Mesh: [1x1 FEMesh]

for n = 1:numel(thermalresults.SolutionTimes)
    figure
    pdeplot3D(thermalmodel,'ColorMapData',thermalresults.Temperature(:,n))
    title(['Temperature at Time = ' num2str(tlist(n))])
    caxis([0 100])
end

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

structuralmodel.Geometry = gm;

structuralmodel.Mesh = mesh;

structuralProperties(structuralmodel,'YoungsModulus',1e10, ...
                                     'PoissonsRatio',0.3, ...'
                                     'CTE',11.7e-6);

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

structuralBodyLoad(structuralmodel,'Temperature',thermalresults);

structuralmodel.ReferenceTemperature = 10;

thermalstressresults = solve(structuralmodel);

pdeplot3D(structuralmodel,'ColorMapData',thermalstressresults.Displacement.Magnitude, ...
                          'Deformation',thermalstressresults.Displacement)
title(['Thermal Expansion at Solution Time = ' num2str(tlist(end))])
caxis([0 3e-3])

for n = 1:numel(thermalresults.SolutionTimes)
    structuralBodyLoad(structuralmodel,'Temperature',thermalresults,'TimeStep',n);
    thermalstressresults = solve(structuralmodel);
    figure
    pdeplot3D(structuralmodel,'ColorMapData', ...
                               thermalstressresults.Displacement.Magnitude, ...
                              'Deformation', ...
                               thermalstressresults.Displacement)
    title(['Thermal Results at Solution Time = ' num2str(tlist(n))])
    caxis([0 3e-3])
end

Решите для переходного отклика в центре 3-D пучка под гармонической нагрузкой на один из его углов.

modelM = createpde('structural','modal-solid');

gm = multicuboid(0.05,0.003,0.003);
modelM.Geometry = gm;
pdegplot(modelM,'EdgeLabels','on','VertexLabels','on');
view([95 5])

msh = generateMesh(modelM);

structuralProperties(modelM,'YoungsModulus',210E9, ...
                            'PoissonsRatio',0.3, ...
                            'MassDensity',7800);

structuralBC(modelM,'Edge',4,'Constraint','fixed');
structuralBC(modelM,'Vertex',7,'Constraint','fixed');

Rm = solve(modelM,'FrequencyRange',[-0.1,500000]);

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

modelD.Geometry = gm;
modelD.Mesh = msh;

structuralProperties(modelD,'YoungsModulus',210E9, ...
                            'PoissonsRatio',0.3, ...
                            'MassDensity',7800);

structuralBC(modelD,'Edge',4,'Constraint','fixed');
structuralBC(modelD,'Vertex',7,'Constraint','fixed');

structuralBoundaryLoad(modelD,'Vertex',5,'Force',[0,0,10],'Frequency',7600);

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

modelD.SolverOptions.RelativeTolerance = 1E-5;
modelD.SolverOptions.AbsoluteTolerance = 1E-9;

tlist = linspace(0,0.004,120);
Rdm = solve(modelD,tlist,'ModalResults',Rm)

Rdm = 
  TransientStructuralResults with properties:

     Displacement: [1x1 FEStruct]
         Velocity: [1x1 FEStruct]
     Acceleration: [1x1 FEStruct]
    SolutionTimes: [1x120 double]
             Mesh: [1x1 FEMesh]

intrpUdm = interpolateDisplacement(Rdm,0,0,0.0015);

plot(Rdm.SolutionTimes,intrpUdm.uz)
grid on
xlabel('Time');
ylabel('Center of beam displacement')

Решите электромагнитную проблему и найдите электрический потенциал и распределение поля для 2-D геометрии, представляющей пластину с отверстием.

Создание электромагнитной модели для электростатического анализа.

emagmodel = createpde('electromagnetic','electrostatic');

Импорт и печать геометрии, представляющей пластину с отверстием.

importGeometry(emagmodel,'PlateHolePlanar.stl');
pdegplot(emagmodel,'EdgeLabels','on')

Укажите диэлектрическую проницаемость вакуума в системе единиц измерения СИ.

emagmodel.VacuumPermittivity = 8.8541878128E-12;

Укажите относительную диэлектрическую проницаемость материала.

electromagneticProperties(emagmodel,'RelativePermittivity',1);

Примените граничные условия напряжения к кромкам, обрамляющим прямоугольник и окружность.

electromagneticBC(emagmodel,'Voltage',0,'Edge',1:4);
electromagneticBC(emagmodel,'Voltage',1000,'Edge',5);

Укажите плотность заряда для всей геометрии.

electromagneticSource(emagmodel,'ChargeDensity',5E-9);

Создайте сетку.

generateMesh(emagmodel);

Решите модель.

R = solve(emagmodel)

R = 
  ElectrostaticResults with properties:

      ElectricPotential: [1218x1 double]
          ElectricField: [1x1 FEStruct]
    ElectricFluxDensity: [1x1 FEStruct]
                   Mesh: [1x1 FEMesh]

Постройте график электрического потенциала и поля.

pdeplot(emagmodel,'XYData',R.ElectricPotential, ...
                  'FlowData',[R.ElectricField.Ex ...
                              R.ElectricField.Ey])
axis equal