www.gusucode.com > pde 案例源码 matlab代码程序 > pde/clampedBeamDynamics.m
%% Dynamic Analysis of a Clamped Beam
% This example shows the Partial Differential Equation Toolbox(TM) analysis
% of the dynamic behavior of a beam clamped at both ends and loaded with a
% uniform pressure load. The pressure load is suddenly applied at time
% equal zero and the magnitude is high enough to produce deflections on the
% same order as the beam thickness.
%
% Accurately predicting this type of behavior requires a geometrically-nonlinear
% formulation of the elasticity equations. One of the main purposes
% of this example is to show how PDE Toolbox can be used to solve a problem
% in nonlinear elasticity. The analysis will be performed with both linear
% and nonlinear formulations to demonstrate the importance of the latter.
%
% This example specifies values of parameters using the imperial system of
% units. You can replace them with values specified in the metric system.
% If you do so, ensure that you specify all values throughout the example
% using the same system.
%
% Copyright 2014-2015 The MathWorks, Inc.
%% Equations
% This section describes the equations of geometrically nonlinear elasticity.
% One approach to handling the large deflections is to consider
% the elasticity equations in the deformed position. However, PDE Toolbox
% formulates the equations based on the original geometry. This motivates
% using a Lagrangian formulation of nonlinear elasticity where stresses,
% strains, and coordinates refer to the original geometry.
%
% The Lagrangian formulation of the equilibrium equations can be written
%
% $$ \rho\ddot u -\nabla \cdot (F \cdot S) = f $$
%
% where $\rho$ is the material density, $u$ is the displacement vector,
% $F$ is the deformation gradient, $S$ is the second Piola-Kirchoff
% stress tensor, and $f$ is the body force vector.
% This equation can also be written in the following tensor form:
%
% $$
% \rho\ddot u_i -\frac{\partial}{\partial x_j}\left(
% \left(\frac{\partial u_i}{\partial x_k}
% + \delta_{ik}\right) S_{kj} \right) = f_i
% $$
%
% Although large deflections are accounted for in this formulation, it is
% assumed that the strains remain small so that linear elastic constitutive
% relations apply. Also, the material is assumed to be isotropic. For the
% 2D plane stress case, the constitutive relations may be written in matrix
% form:
%
% $$
% \left\{
% \begin{array}{c}
% S_{11} \\ S_{22} \\ S_{12}
% \end{array} \right\} =
% \left[
% \begin{array}{ccc}
% C_{11} & C_{12} & \\
% C_{12} & C_{22} & \\
% & & 2 G_{12} \\
% \end{array} \right]
% \left\{
% \begin{array}{c}
% E_{11} \\ E_{22} \\ E_{12}
% \end{array} \right\}
% $$
%
% where $E_{ij}$ is the Green-Lagrange strain tensor defined as
%
% $$
% E_{ij} = \frac{1}{2} \left(
% \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} +
% \frac{\partial u_k}{\partial x_i} \frac{\partial u_k}{\partial x_j}
% \right)
% $$
%
% For an isotropic material
%
% $$ C_{11} = C_{22} = \frac{E}{1 - \nu^2} $$
%
% $$ C_{12} = \frac{E\nu}{1 - \nu^2} $$
%
% $$ G_{12} = \frac{E}{2(1 + \nu)} $$
%
% where $E$ is Young's modulus and $\nu$ is Poisson's ratio.
%
% Readers who are interested in more details about the Lagrangian
% formulation for nonlinear elasticity can consult, for example, reference 1.
%
% The equations presented above completely define the geometrically
% nonlinear plane stress problem. The work required to convert them to
% a form acceptable to PDE Toolbox is considerably simplified by using
% the MATLAB Symbolic Math Toolbox. Symbolic Math Toolbox can perform the
% necessary algebraic manipulations and output a MATLAB function defining
% the c-coefficient that can be passed to PDE Toolbox functions. This
% function, |cCoefficientLagrangePlaneStress|, is shown in the appendix below.
%% Create the PDE Model
N = 2; % Two PDE in plane stress elasticity
model = createpde(N);
%% Define the Geometry
blength = 5; % Beam length, in.
height = .1; % Thickness of the beam, in.
%%
% A drawing of the clamped beam is shown in the figure below.
%
% <<../clamped_beam_drawing.png>>
%%
% Since the beam geometry and loading are symmetric about the beam
% center(x = blength/2), the model can be simplified by considering
% only the right-half of the beam.
l2 = blength/2;
h2 = height/2;
% Create the edges of the rectangle representing the beam with these
% two statements:
rect = [3 4 0 l2 l2 0 -h2 -h2 h2 h2]';
g = decsg(rect,'R1',('R1')');
% The geometryFromEdges function creates a geometry object from the edges
% and stores it within the model.
pg = geometryFromEdges(model,g);
%%
% Plot the geometry and display the edge labels.
% The edge labels are needed for edge identification when applying
% boundary conditions.
figure
pdegplot(g,'EdgeLabels','on');
title('Geometry With Edge Labels Displayed');
axis([-.1 1.1*l2 -5*h2 5*h2]); % Scale the plot so the labels are viewable
%% Specify Equation Coefficients
% Derive the equation coefficients using the material properties.
% For the linear case, the c-coefficient matrix is constant
E = 3.0e7; % Young's modulus of the material, lbs/in^2
gnu = .3; % Poisson's ratio of the material
rho = .3/386; % Density of the material
G = E/(2.*(1+gnu));
mu = 2*G*gnu/(1-gnu);
c = [2*G+mu; 0; G; 0; G; mu; 0; G; 0; 2*G+mu];
f = [0 0]'; % No body forces
specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0, 'f', f);
%% Apply the Boundary Conditions
% Symmetry condition is x-displacement equal zero at left edge.
symBC = applyBoundaryCondition(model,'mixed','Edge',4,'u',0,'EquationIndex',1);
%%
% x- and y-displacements equal zero along right edge.
clampedBC = applyBoundaryCondition(model,'dirichlet','Edge',2,'u',[0 0]);
%%
% Apply a constant z-direction stress along the top edge.
sigma = 2e2;
presBC = applyBoundaryCondition(model,'neumann','Edge',3,'g',[0 sigma]);
%% Set the Initial Conditions
% Zero initial displacements and velocities
setInitialConditions(model,0,0);
%% Create the Mesh
% Create a mesh with approximately eight elements through the thickness
% of the beam.
generateMesh(model,'Hmax',height/8,'MesherVersion','R2013a');
%% Linear Solution
% Set up the analysis timespan
tfinal = 3e-3; % Final time in the analysis
tlist = linspace(0,tfinal,100); % Save the output at 100 time points
% Compute the time-dependent solution
result = solvepde(model,tlist);
%%
% Interpolate the solution at the geometry center for the
% y-component (component 2) at all times.
xc = 1.25;
yc = 0;
u4Linear = interpolateSolution(result,xc,yc,2,1:length(tlist));
%% Specify Equation Coefficients for Nonlinear Solution
% The function |cCoefficientLagrangePlaneStress| takes the isotropic material
% properties and location and state structures, and returns a c-matrix
% for a nonlinear plane-stress analysis. Small strains are assumed; i.e. E
% and $\nu$ are independent of the solution. PDE Toolbox calls
% user-defined coefficient functions with the arguments location and state.
% The function |cCoefficientLagrangePlaneStress| expects the arguments
% E, gnu, location, state. |c| is defined below as an anonymous function to
% provide an interface between these two function signatures. (The function
% |cCoefficientLagrangePlaneStress| can be used with any geometric nonlinear
% plane stress analysis of a model made from an isotropic material.)
c = @(location, state) cCoefficientLagrangePlaneStress(E,gnu,location,state);
specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0 , 'f', f);
%% Nonlinear Solution
% Compute the time-dependent solution.
result = solvepde(model,tlist);
%%
% As before, interpolate the solution at the geometry center for the
% y-component (component 2) at all times.
u4NonLinear = interpolateSolution(result,xc,yc,2,1:length(tlist));
%% Plot Solutions
% The figure below shows the y-deflection at the center of the beam
% as a function of time. The nonlinear analysis computes displacements
% that are substantially less than the linear analysis. This "stress
% stiffening" effect is also reflected in the higher oscillation frequency
% from the nonlinear analysis.
figure
plot(tlist,u4Linear(:),tlist,u4NonLinear(:));
legend('Linear','Nonlinear');
title 'Deflection at Beam Center'
xlabel 'Time, seconds'
ylabel 'Deflection, inches'
grid on
%% References
% # Malvern, Lawrence E., Introduction to the Mechanics of a Continuous
% Medium, Prentice Hall, 1969.
%% Appendix - Nonlinear C-Coefficient Function
%
% The function |cCoefficientLagrangePlaneStress|
% calculates the c-coefficient matrix for a large displacement Lagrangian
% plane stress formulation.
type cCoefficientLagrangePlaneStress
displayEndOfDemoMessage(mfilename)