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

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

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

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

Решите электромагнитную задачу и найдите электрический потенциал и распределение поля для 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);

Сгенерируйте 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