Page 1 of 18. Copyright J.E. Akin. All rights reserved. Matlab Plane Stress Example (Draft 2, April 9, 2007) Introduction Here the Matlab closed form element matrices for the T3 element (3 node triangle, constant stress) is illustrated for a square plate, 2 by 2 inches. It is fixed at the top left corner, is restrained from horizontal (but not vertical) displacement at its bottom lect corner. It is loaded only by a horizontal body force load acting to the right. Away from the stress concentrations at the two corners it is essentially a 1-D problem. If the whole left edge were restrained againts horizontal motion it would correspond to a axial bar hanging under its own weight. Then the free end deflection is δ = W Lx /(2 A E) and the axial stress varies from σmax = W / A at the support to zero at the free end. Here Lx is the horizontal length, W the weight, E the elastic modulus, A = Lx t is the area for a given thickness (t=0.005 here). Due to Poisson’s ratio, the solution here has a vertical (y) displacement. The above plot shows the plate and a finer mesh solution for the resultant displacements, from CosmosWorks. The stress in the y-direction should be zero except at the two stress concentrations.Page 2 of 18. Copyright J.E. Akin. All rights reserved. Matlab solution Here a crude mesh is formed by using the two diagonal lines to form four elements that are three node (constant strain) triangles, CST. There are two displacement components (degrees of freedom, dof) at each node. There are three strains and thress stresses to be determined. The execution and sample plots will be shown first. The modular source script is listed at the end of this document, and as a downloadable file on the class web site. The nodal data are stored in the file msh_bc_xyz.tmp, the element type and connectivity in msh_typ_nodes.tmp, the essential (displacement component) boundary condition data in msh_ebc.tmp. If point loads or sources existed then file msh_load_pt.tmp would also be present. Those files are used for data validation plots as well as input to the stress calculations. The current example starts in Unix by invoking Matlab and running mesh plotting options to check the data. % Matlab >> addpath /net/course-a/mech517/public_html/Matlab_Plots >> mesh_shrink_plot Read 5 mesh coordinate pairs,4 elements with 3 nodes each >> bc_flags_plot This displays the packed binary code for each node that has an essential boundary condition. Since there are two displacement components here there are two digits in the packed integer. A one denotes true (a restraint exists) while a zero denotes a free displacement. They are ordered as Ux (horizontal) and Uy (vertical) components. For this application the allowed flags are 00, 10, 01, or 11.Page 3 of 18. Copyright J.E. Akin. All rights reserved. Having validated the input geometry and boundary condition flags the stress calculations are invoked: >> Modular_Plane_Sress_XY Read 5 nodes with bc_flag & 2 coordinates. (file msh_bc_xyz.tmp) 10 0 0 11 0 2 00 1 1 00 2 0 00 2 2 Read 4 elements with type number & 3 nodes each. (file msh_typ_nodes.tmp) Maximum number of element types = 1. 1 1 3 2 1 1 4 3 1 4 5 3 1 3 5 2 Read 3 EBC data with Node, DOF, Value. (file msh_ebc.tmp) 1 1 0 2 1 0 2 2 0 Application properties are: Thickness = 0.005 Body Force = 500000 0 Elastity matrix: 1.0e+10 * 1.6000 0.4000 0 0.4000 1.6000 0 0 0 0.6000Page 4 of 18. Copyright J.E. Akin. All rights reserved. Node, DOF, Resultant Load Value 1 1 1666.67 2 1 1666.67 3 1 3333.33 4 1 1666.67 5 1 1666.67 Displacement Solution: X_disp Y_disp at 5 nodes 1.0e-04 * 0 0.3403 0 0 0.5243 0.1701 0.6667 0.1667 0.6667 0.1736 Node, DOF, Reaction Force Value 1 1 -5000 2 1 -5000 2 2 -9.37916e-13 Elem, QP, X_qp, Y_qp (Location of stresses) Elem, QP, Stress_qp: xx yy xy 1 1 0.333333 1 1 1 770833 -62500 1.01644e-10 2 1 1 0.333333 2 1 500000 -5.82077e-11 62500 3 1 1.66667 1 3 1 229167 62500 -2.03288e-11 4 1 1 1.66667 4 1 500000 -2.91038e-11 -62500 End of stress analysis. Now plot various results. >> quiver_resultant_load_mesh(0.5) Using a scale of 0.5 and vector increment of 1Page 5 of 18. Copyright J.E. Akin. All rights reserved. >> quiver_disp_vec_mesh(0.25) >> quiver_reaction_vec_mesh(0.5) Using a scale of 0.5 and vector increment of 1Page 6 of 18. Copyright J.E. Akin. All rights reserved. >> color_result_surface (resultant displacement magnitude) >> result_shrink_surface (resultant displacement magnitude)Page 7 of 18. Copyright J.E. Akin. All rights reserved. >> color_result(1) (X-displacement component) Max value is 6.66667e-05 at node 4 Min value is 0 at node 1 >> color_result(2) ((X-displacement component)-displacement component) Max value is 3.40278e-05 at node 1 Min value is 0 at node 2Page 8 of 18. Copyright J.E. Akin. All rights reserved. >> contour_result_on_mesh (resultant displacement magnitude) >> smooth_result_w_bar (resultant displacement magnitude) Max value 6.88902e-05 at node 5,Min value 0 at node 2 >> quitPage 9 of 18. Copyright J.E. Akin. All rights reserved. Source listing function Modular_Plane_Stress_XY (load_pt, pre_p, pre_e) %............................................................. % Plane Stress with body and point loads, T3 triangle % XY COORDINATES CLOSED FORM INTEGRALS % Change E_e to get Plane Strain %............................................................. % pre_e = # of dummy items before el_type & connectivity % pre_p = # of dummy items before BC_flag % coordinates if ( nargin == 2 ) ; % check for optional data pre_e = 0 ; % no el # in msh_typ_nodes elseif ( nargin == 1 ) ; % check for optional data pre_p = 0 ; pre_e = 0 ; % no pt # in msh_bc_xyz elseif ( nargin == 0 ) ; % check for optional data load_pt = 0 ; pre_p = 0 ; pre_e = 0 ; % no point source data end % if from argument count % Application and element dependent controls