This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

solve

Solve heat transfer or structural analysis problem

Syntax

thermalresults = solve(thermalmodel)
thermalresults = solve(thermalmodel,tlist)
structuralresults = solve(structuralmodel)
structuralresults = solve(structuralmodel,tlist)
structuralresults = solve(structuralmodel,'FrequencyRange',[omega1,omega2])

Description

example

thermalresults = solve(thermalmodel) returns the solution to the steady-state thermal model represented in thermalmodel.

example

thermalresults = solve(thermalmodel,tlist) returns the solution to the transient thermal model represented in thermalmodel at the times tlist.

example

structuralresults = solve(structuralmodel) returns the solution to the static structural analysis model represented in structuralmodel.

example

structuralresults = solve(structuralmodel,tlist) returns the solution to the transient structural dynamics model represented in structuralmodel.

example

structuralresults = solve(structuralmodel,'FrequencyRange',[omega1,omega2]) returns the solution to the modal analysis model for all modes in the frequency range [omega1,omega2]. Define omega1 as slightly smaller than the lowest expected frequency and omega2 as slightly larger than the highest expected frequency. For example, is the lowest expected frequency is zero, then use a small negative value for omega1.

Examples

collapse all

Solve a 3-D steady-state thermal problem.

Create a thermal model for this problem.

thermalmodel = createpde('thermal');

Import and plot the block geometry.

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

Assign material properties.

thermalProperties(thermalmodel,'ThermalConductivity',80);

Apply a constant temperature of to the left side of the block (face 1) and a constant temperature of to the right side of the block (face 3). All other faces are insulated by default.

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

Mesh the geometry and solve the problem.

generateMesh(thermalmodel);
thermalresults = solve(thermalmodel)
thermalresults = 
  SteadyStateThermalResults with properties:

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

The solver finds the temperatures and temperature gradients at the nodal locations. To access these values, use thermalresults.Temperature, thermalresults.XGradients, and so on. For example, plot temperatures at nodal locations.

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

Solve a 2-D transient thermal problem.

Create a transient thermal model for this problem.

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

Create the geometry and include it in the model.

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

For the square region, assign these thermal properties:

  • Thermal conductivity is .

  • Mass density is .

  • Specific heat is .

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

For the diamond region, assign these thermal properties:

  • Thermal conductivity is .

  • Mass density is .

  • Specific heat is .

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

Assume that the diamond-shaped region is a heat source with a density of .

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

Apply a constant temperature of to the sides of the square plate.

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

Set the initial temperature to .

thermalIC(thermalmodel,0);

Mesh the geometry.

generateMesh(thermalmodel);

The dynamics for this problem are very fast. The temperature reaches a steady state in about 0.1 second. To capture the interesting part of the dynamics, set the solution time to logspace(-2,-1,10). This command returns 10 logarithmically spaced solution times between 0.01 and 0.1.

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

Solve the equation.

thermalresults = solve(thermalmodel,tlist)
thermalresults = 
  TransientThermalResults with properties:

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

Plot the solution with isothermal lines by using a contour plot.

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

Solve a static structural model representing a bimetallic cable under tension.

Create a static structural model for solving a solid (3-D) problem.

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

Create the geometry and include it in the model. Plot the geometry.

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

Specify the Young's modulus and Poisson's ratio for each metal.

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

Specify that faces 1 and 4 are fixed boundaries.

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

Specify the surface traction for faces 2 and 5.

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

Generate a mesh and solve the problem.

generateMesh(structuralmodel);
structuralresults = solve(structuralmodel)
structuralresults = 
  StaticStructuralResults with properties:

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

The solver finds the values of displacement, stress, strain, and von Mises stress at the nodal locations. To access these values, use structuralresults.Displacement, structuralresults.Stress, and so on. Displacement, stress, and strain at the nodal locations are structure arrays with fields representing their components.

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]

Plot the deformed shape with the z-component of normal stress.

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

Solve for transient response of a thin 3-D plate under a harmonic load at the center.

Create a transient dynamic model for a 3-D problem.

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

Create the geometry and include it in the model. Plot the geometry.

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

Zoom in to see the face labels on the small plate in the center.

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

Specify the Young's modulus, Poisson's ratio, and mass density of the material.

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

Specify that all faces on the periphery of the thin 3-D plate are fixed boundaries.

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

Apply a sinusoidal pressure load on the small face at the center of the plate.

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

Generate a mesh with linear elements.

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

Specify the zero initial displacement and velocity.

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

Solve the model.

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]

The solver finds the values of displacement, velocity, and acceleration at the nodal locations. To access these values, use structuralresults.Displacement, structuralresults.Velocity, and so on. Displacement, velocity, and acceleration are structure arrays with fields representing their components.

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]

Find the fundamental (lowest) mode of a 2-D cantilevered beam, assuming a prevalence of the plane-stress condition.

Specify the following geometric and structural properties of the beam, along with a unit plane-stress thickness.

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

Create a model plane-stress model, assign a geometry, and generate a 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);

Define a maximum element size (five elements through the beam thickness).

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

Specify the structural properties and boundary constraints.

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

Compute the analytical fundamental frequency (Hz) using the beam theory.

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

Specify a frequency range that includes an analytically computed frequency and solve the model.

modalresults = solve(structuralmodel,'FrequencyRange',[0,1e6])
modalresults = 
  ModalStructuralResults with properties:

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

The solver finds natural frequencies and modal displacement values at nodal locations. To access these values, use modalresults.NaturalFrequencies and 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]

Plot the y-component of the solution for the fundamental frequency.

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

Find the deflection of a 3-D cantilever beam under a nonuniform thermal load. Specify the thermal load on the structural model using the solution from a transient thermal analysis on the same geometry and mesh.

Transient Thermal Model Analysis

Create a transient thermal model.

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

Create and plot the geometry.

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

Generate a mesh.

mesh = generateMesh(thermalmodel);

Specify the thermal properties of the material.

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

Specify the constant temperatures applied to the left and right ends on the beam.

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

Specify the heat source over the entire geometry.

internalHeatSource(thermalmodel,10);

Set the initial temperature.

thermalIC(thermalmodel,0);

Solve the model.

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]

Plot the temperature distribution for each time step.

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

Structural Analysis with Thermal Load

Create a static structural model.

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

Include the same geometry as for the thermal model.

structuralmodel.Geometry = gm;

Use the same mesh that you used to obtain the thermal solution.

structuralmodel.Mesh = mesh;

Specify the Young's modulus, Poisson's ratio, and coefficient of thermal expansion.

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

Apply a fixed boundary condition on face 5.

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

Apply a body load using the transient thermal model solution. By default, structuralBodyLoad uses the solution for the last time step.

structuralBodyLoad(structuralmodel,'Temperature',thermalresults);

Specify the reference temperature.

structuralmodel.ReferenceTemperature = 10;

Solve the structural model.

thermalstressresults = solve(structuralmodel);

Plot the deformed shape of the beam corresponding to the last step of the transient thermal model solution.

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

Now specify the body loads as the thermal model solutions for all time steps. For each body load, solve the structural model and plot the corresponding deformed shape of the beam.

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

Input Arguments

collapse all

Thermal model, specified as a ThermalModel object. The model contains the geometry, mesh, thermal properties of the material, internal heat source, boundary conditions, and initial conditions.

Example: thermalmodel = createpde('thermal','steadystate')

Structural model, specified as a StructuralModel object. The model contains the geometry, mesh, structural properties of the material, body loads, boundary loads, and boundary conditions.

Example: structuralmodel = createpde('structural','transient-solid')

Solution times, specified as a real vector of monotonically increasing or decreasing values.

Example: 0:20

Data Types: double

Frequency range, specified as a vector of two elements. Define omega1 as slightly smaller than the lowest expected frequency and omega2 as slightly larger than the highest expected frequency. For example, is the lowest expected frequency is zero, then use a small negative value for omega1.

Example: [-0.1,1000]

Data Types: double

Output Arguments

collapse all

Thermal results, returned as a SteadyStateThermalResults object or TransientThermalResults object. The type of thermalresults depends on whether thermalmodel represents a steady-state problem (thermalmodel.AnalysisType = 'steadystate') or a transient problem (thermalmodel.AnalysisType = 'transient').

Structural results, returned as a StaticStructuralResults, TransientStructuralResults, or ModalStructuralResults object. The type of structuralresults depends on whether structuralmodel represents a static problem, a transient problem, or a modal analysis problem. To check the type of a structural analysis problem, type structuralmodel.AnalysisType.

Introduced in R2017a