Решите 3-D установившуюся тепловую задачу.

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

thermalmodel = createpde('thermal');

Импортируйте и постройте геометрию блока.

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

Присвойте свойства материала.

thermalProperties(thermalmodel,'ThermalConductivity',80);

Примените постоянную температуру $$\mathrm{}\mathrm{}100{}^{\circ }C$$ к левой стороне блока (стоят 1), и постоянная температура $$\mathrm{}\mathrm{}300{}^{\circ }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);

Примите, что ромбовидная область является источником тепла с плотностью $$4W/m^3$$.

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

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

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

Установите начальную температуру на $$0{}^{\circ }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]);

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

generateMesh(structuralmodel);
structuralresults = solve(structuralmodel)

structuralresults = 
  StaticStructuralResults with properties:

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

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

structuralresults.Displacement

ans = struct with fields:
           ux: [22281x1 double]
           uy: [22281x1 double]
           uz: [22281x1 double]
    Magnitude: [22281x1 double]

structuralresults.Stress

ans = struct with fields:
    sxx: [22281x1 double]
    syy: [22281x1 double]
    szz: [22281x1 double]
    syz: [22281x1 double]
    sxz: [22281x1 double]
    sxy: [22281x1 double]

structuralresults.Strain

ans = struct with fields:
    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 struct]
         Velocity: [1x1 struct]
     Acceleration: [1x1 struct]
    SolutionTimes: [1x300 double]
             Mesh: [1x1 FEMesh]

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

structuralresults.Displacement

ans = struct with fields:
           ux: [1873x300 double]
           uy: [1873x300 double]
           uz: [1873x300 double]
    Magnitude: [1873x300 double]

structuralresults.Velocity

ans = struct with fields:
           vx: [1873x300 double]
           vy: [1873x300 double]
           vz: [1873x300 double]
    Magnitude: [1873x300 double]

structuralresults.Acceleration

ans = struct with fields:
           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 struct]
                  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 = struct with fields:
    ux: [6511x32 double]
    uy: [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 struct]
               Velocity: [1x1 struct]
           Acceleration: [1x1 struct]
    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 struct]
         Velocity: [1x1 struct]
     Acceleration: [1x1 struct]
    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')