Решите 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)

Решите 2D переходную тепловую задачу.

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

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{W}/{\mathrm{m}}^2$.

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

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

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

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

thermalIC(thermalmodel,0);

Сгенерируйте mesh.

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]);

Сгенерируйте mesh и решите задачу.

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);

Сгенерируйте mesh с линейными элементами.

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]

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

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

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

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

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

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')

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

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

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);

Сгенерируйте mesh.

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