7 Laminate Analysis -Part I 7.1 Basic Equations Fiber-reinforced materials consist usually of multiple layers of material to form a laminate.Each layer is thin and may have a different fiber orientation-see Fig.7.1. Two laminates may have the same number of layers and the same fiber angles but the two laminates may be different because of the arrangement of the layers. In this chapter,we will evaluate the influence of fiber directions,stacking arrange- ments and material properties on laminate and structural response.We will study a simplified theory called classical lamination theory for this purpose(see [1]). Figure 7.2 shows a global Cartesian coordinate system and a general laminate consisting of N layers.The laminate thickness is denoted by HI and the thickness of an individual layer by h.Not all layers necessarily have the same thickness,so the thickness of the kth layer is denoted by hk. The origin of the through-thickness coordinate,designated z,is located at the laminate geometric midplane.The geometric midplane may be within a particular layer or at an interface between layers.We consider the +z axis to be downward and the laminate extends in the z direction from-H/2 to +H/2.We refer to the layer at the most negative location as layer 1,the next layer in as layer 2,the layer at an arbitrary location as layer k,and the layer at the most positive z position as layer N.The locations of the layer interfaces are denoted by a subscripted z;the 0000000000000000 0000000000000000 Fig.7.1.Schematic illustration of a laminate with four layers
7 Laminate Analysis – Part I 7.1 Basic Equations Fiber-reinforced materials consist usually of multiple layers of material to form a laminate. Each layer is thin and may have a different fiber orientation – see Fig. 7.1. Two laminates may have the same number of layers and the same fiber angles but the two laminates may be different because of the arrangement of the layers. In this chapter, we will evaluate the influence of fiber directions, stacking arrangements and material properties on laminate and structural response. We will study a simplified theory called classical lamination theory for this purpose (see [1]). Figure 7.2 shows a global Cartesian coordinate system and a general laminate consisting of N layers. The laminate thickness is denoted by H and the thickness of an individual layer by h. Not all layers necessarily have the same thickness, so the thickness of the kth layer is denoted by hk. The origin of the through-thickness coordinate, designated z, is located at the laminate geometric midplane. The geometric midplane may be within a particular layer or at an interface between layers. We consider the +z axis to be downward and the laminate extends in the z direction from −H/2 to +H/2. We refer to the layer at the most negative location as layer 1, the next layer in as layer 2, the layer at an arbitrary location as layer k, and the layer at the most positive z position as layer N. The locations of the layer interfaces are denoted by a subscripted z; the Fig. 7.1. Schematic illustration of a laminate with four layers
116 7 Laminate Analysis-Part I y N (a) layer N 5 C (b) Fig.7.2.Schematic illustration showing a cross-section and a plan view first layer is bounded by locations zo and z1,the second layer by z1 and z2,the kth layer by zk-1 and zk,and the Nth layer by zN-1 and zN [1]. Let us examine the deformation of an r-z cross-section [1].Figure 7.3 shows in detail the deformation of a cross-section,and in particular the displacements of point P,a point located at an arbitrary distance z below point po,a point on the reference surface,with points P and po being on line AA'.The superscript 0 will be reserved to denote the kinematics of point po on the reference surface.In particular,the horizontal translation of point po in the x direction will be denoted by u.The vertical translation will be denoted by w.The rotation of the reference surface about the y axis at point po is ow/Or.An important part of the Kirchhoff hypothesis is the assumption that line AA'remains perpendicular to the reference surface.Because of this,the rotation of line AA'is the same as the rotation of the reference surface,and thus the rotation of line AA',as viewed in the r-z plane,is 0w/Or.It is assumed that [1]: 0w0 Ox <1 (7.1) By less than unity is meant that sines and tangents of angles of rotation are replaced by the rotations themselves,and cosines of the angles of rotation are replace by 1.With this approximation,then,the rotation of point po causes point P to translate horizontally in the minus z direction by an amount equal to: 8w (7.2)
116 7 Laminate Analysis – Part I Fig. 7.2. Schematic illustration showing a cross-section and a plan view first layer is bounded by locations z0 and z1, the second layer by z1 and z2, the kth layer by zk−1 and zk, and the Nth layer by zN−1 and zN [1]. Let us examine the deformation of an x-z cross-section [1]. Figure 7.3 shows in detail the deformation of a cross-section, and in particular the displacements of point P, a point located at an arbitrary distance z below point P0, a point on the reference surface, with points P and P0 being on line AA . The superscript 0 will be reserved to denote the kinematics of point P0 on the reference surface. In particular, the horizontal translation of point P0 in the x direction will be denoted by u0. The vertical translation will be denoted by w0. The rotation of the reference surface about the y axis at point P0 is ∂w0/∂x. An important part of the Kirchhoff hypothesis is the assumption that line AA remains perpendicular to the reference surface. Because of this, the rotation of line AA is the same as the rotation of the reference surface, and thus the rotation of line AA , as viewed in the x-z plane, is ∂w0/∂x. It is assumed that [1]: ∂w0 ∂x < 1 (7.1) By less than unity is meant that sines and tangents of angles of rotation are replaced by the rotations themselves, and cosines of the angles of rotation are replace by 1. With this approximation, then, the rotation of point P0 causes point P to translate horizontally in the minus x direction by an amount equal to: z = ∂w0 ∂x (7.2)
7.1 Basic Equations 117 e r Fig.7.3.Schematic illustration showing the kinematics of deformation of a laminate Therefore,the horizontal translation of a point P with coordinates (z,y,z)in the direction of the r-axis is then given by: u,,习=ue,0-00e, (7.3) Also,the vertical translation of point P in the direction of the z-axis is given by: w(x,y,z)=w(,y) (7.4) The horizontal translation of point P in the direction of the y-axis is similar to that in the direction of the z-axis and is given by: u(z,y,z)=v°(x,)-z9 0u°(x,y） (7.5) oy Therefore,we now have the following relations: u(x,y,z)=u°(x,)）-z 0w°(x, (7.6a)
7.1 Basic Equations 117 Fig. 7.3. Schematic illustration showing the kinematics of deformation of a laminate Therefore, the horizontal translation of a point P with coordinates (x, y, z) in the direction of the x-axis is then given by: u(x, y, z) = u0 (x, y) − z ∂w0(x, y) ∂x (7.3) Also, the vertical translation of point P in the direction of the z-axis is given by: w(x, y, z) = w0 (x, y) (7.4) The horizontal translation of point P in the direction of the y-axis is similar to that in the direction of the x-axis and is given by: v(x, y, z) = v0 (x, y) − z ∂w0(x, y) ∂y (7.5) Therefore, we now have the following relations: u(x, y, z) = u0 (x, y) − z ∂w0(x, y) ∂x (7.6a)
118 7 Laminate Analysis-Part I v(e,y,z）=v(x,y)-z 0u°(t,y) (7.6b) w(x,y,z)=w°(x,) (7.6c) Next,we investigate the strains that result from the displacements according to the Kirchhoff hypothesis.This can be done by using the strain-displacement relations from the theory of elasticity.Using these relations and(7.6a,b,c),we can compute the strains at any point within the laminate,and by using these laminate strains in the stress-strain relations,we can compute the stresses at any point within the laminate. From the strain-displacement relations and (7.6a),the extensional strain in the z direction,Er,is given by: ez,,动=0u红y,=0ue,边-0uz,边 (7.7) Ox Ox 0x2 Equation (7.7)may be re-written as follows: er(x,,z)=e(x,)+z9(x,) (7.8) where the following notation is used: (红，)=0ue (7.9a) (红，0= 82w°(x,y） (7.9b) 0x2 The quantity is referred to as the extensional strain of the reference surface in the r direction,and is referred to as the curvature of the reference surface in the z direction.The other five strain components are given by: E(x,y,z)三 0u(x,2=8(红，）)+z(x,) (7.10a) 则 e:e,)=0咖=0ug业=0 (7.10b) :l,斯，=加e8+m色 0z -0mE边_0wz边=0 (7.10c) ay oy x:(x,y,2）≡ 0(,y,+0u(x,,】 0x 8z 0m(E,_u(x,型=0 (7.10d) Ox Ox Yxy(x,,2) 0u(c,,+0ue,42=,+z, (7.10e) Ox 8y where the following notation is used: (z,到= 0u(x,) (7.11a) by
118 7 Laminate Analysis – Part I v(x, y, z) = v0 (x, y) − z ∂w0(x, y) ∂y (7.6b) w(x, y, z) = w0 (x, y) (7.6c) Next, we investigate the strains that result from the displacements according to the Kirchhoff hypothesis. This can be done by using the strain-displacement relations from the theory of elasticity. Using these relations and (7.6a,b,c), we can compute the strains at any point within the laminate, and by using these laminate strains in the stress-strain relations, we can compute the stresses at any point within the laminate. From the strain-displacement relations and (7.6a), the extensional strain in the x direction, εx, is given by: εx(x, y, z) ≡ ∂u(x, y, z) ∂x = ∂u0(x, y) ∂x − z ∂2w0(x, y) ∂x2 (7.7) Equation (7.7) may be re-written as follows: εx(x, y, z) = ε 0 x(x, y) + zκ0 x(x, y) (7.8) where the following notation is used: ε 0 x(x, y) = ∂u0(x, y) ∂x (7.9a) κ0 x(x, y) = −∂2w0(x, y) ∂x2 (7.9b) The quantity ε0 x is referred to as the extensional strain of the reference surface in the x direction, and κ0 x is referred to as the curvature of the reference surface in the x direction. The other five strain components are given by: εy(x, y, z) ≡ ∂v(x, y, z) ∂y = ε 0 y(x, y) + zκ0 y(x, y) (7.10a) εz(x, y, z) ≡ ∂w(x, y, z) ∂z = ∂w0(x, y) ∂z = 0 (7.10b) γyz(x, y, z) ≡ ∂w(x, y, z) ∂y + ∂v(x, y, z) ∂z = ∂w0(x, y) ∂y − ∂w0(x, y) ∂y = 0 (7.10c) γxz(x, y, z) ≡ ∂w(x, y, z) ∂x + ∂u(x, y, z) ∂z = ∂w0(x, y) ∂x − ∂w0(x, y) ∂x = 0 (7.10d) γxy(x, y, z) ≡ ∂v(x, y, z) ∂x + ∂u(x, y, z) ∂y = γ0 xy + zκ0 xy (7.10e) where the following notation is used: ε 0 y(x, y) = ∂v0(x, y) ∂y (7.11a)
7.2 MATLAB Functions Used 119 9(x,)=- 82w°(x, 0y2 (7.11b) 0,e,=0m色+ 0u°(x, (7.11c) Ox by k(x,)=-2 0u°(x, 8xoy (7.11d) The quantitiesandare referred to as the reference surface extensional strain in the y direction,the reference surface curvature in the y direction,the reference surface inplane shear strain,and the reference surface twisting curvature, respectively. The second important assumption of classical lamination theory is that each point within the volume of a laminate is in a state of plane stress.Therefore,we can compute the stresses if we know the strains and curvatures of the reference surface. Accordingly,using the strains that result from the Kirchhoff hypothesis,(7.8)and (7.10a,e),we find that the stress-strain relations for a laminate become: 012 Q16 (7.12) Finally the force and moment resultants in the laminate can be computed using the stresses as follows: H/2 Nr= ordz (7.13a) -H/2 H/2 Ny= (7.13b) -H/2 H/2 (7.13c -H/2 H/2 M:= 0x2d2 (7.13d) -H/2 H/2 My= (7.13e) -H/2 H/2 My= Tryzdz (7.13f) -H/2 7.2 MATLAB Functions Used The only MATLAB function used in this chapter to calculate the strains is:
7.2 MATLAB Functions Used 119 κ0 y(x, y) = −∂2w0(x, y) ∂y2 (7.11b) γ0 xy(x, y) = ∂v0(x, y) ∂x + ∂u0(x, y) ∂y (7.11c) κ0 xy(x, y) = −2∂2w0(x, y) ∂x∂y (7.11d) The quantities ε0 y, κ0 y, γ0 xy, and κ0 xy are referred to as the reference surface extensional strain in the y direction, the reference surface curvature in the y direction, the reference surface inplane shear strain, and the reference surface twisting curvature, respectively. The second important assumption of classical lamination theory is that each point within the volume of a laminate is in a state of plane stress. Therefore, we can compute the stresses if we know the strains and curvatures of the reference surface. Accordingly, using the strains that result from the Kirchhoff hypothesis, (7.8) and (7.10a, e), we find that the stress-strain relations for a laminate become: ⎧ ⎨ ⎩ σx σy τxy ⎫ ⎬ ⎭ = ⎡ ⎢ ⎣ Q¯11 Q¯12 Q¯16 Q¯12 Q¯22 Q¯26 Q¯16 Q¯26 Q¯66 ⎤ ⎥ ⎦ ⎧ ⎪⎨ ⎪⎩ ε0 x + zκ0 x ε0 y + zκ0 y γ0 xy + zκ0 xy ⎫ ⎪⎬ ⎪⎭ (7.12) Finally the force and moment resultants in the laminate can be computed using the stresses as follows: Nx = H/ 2 −H/2 σxdz (7.13a) Ny = H/ 2 −H/2 σydz (7.13b) Nxy = H/ 2 −H/2 τxydz (7.13c) Mx = H/ 2 −H/2 σxzdz (7.13d) My = H/ 2 −H/2 σyzdz (7.13e) Mxy = H/ 2 −H/2 τxyzdz (7.13f) 7.2 MATLAB Functions Used The only MATLAB function used in this chapter to calculate the strains is:
120 7 Laminate Analysis-Part I Strains(eps_ro,eps-yo,gam_ryo,kap_ro,kap-yo,kap_ryo,z)-This function calcu- lates the three strains er,Ey,and yry at any point P on the normal line given the three strainsand the three curvaturesat point po,and the distance z between P and po.There are seven input arguments to this function. The function returns the 3 x 1 strain vector. The following is a listing of the MATLAB source code for this function: function y Strains(eps_xo,eps_yo,gam_xyo,kap_xo,kap_yo,kap_xyo,z) %Strains This function returns the strain vector at any point P along the normal line at distance z from point Po which % lies on the reference surface.There are seven input arguments for this function-namely the three strains and three curvatures at point Po and the distance z. 名 The size of the strain vector is 3 x 1. epsilonx eps_xo +z kap_xo; epsilony eps_yo z kap_yo; gammaxy gam_xyo +z kap_xyo; y [epsilonx epsilony gammaxy]; MATLAB Example 7.1 Consider a graphite-reinforced polymer composite laminate with the elastic con- stants as given in Example 2.2.The laminate has total thickness of 0.500 mm and is stacked as a [0/90s laminate.The four layers are of equal thickness.It is deformed so that at a point (r,y)on the reference surface,we have the following strains and curvatures: €2=400×10-6 9=2g=k2=K9=k2y=0 0 Use MATLAB to determine the following: (a)the three components of strain at the interface locations. (b)the three components of stress in each layer.Plot the stress distribution along the depth of the laminate for each component. (c)the force and moment resultants in the laminate. (d)the three components of strain at the interface locations with respect to the principal material system. (e)the three components of stress in each layer with respect to the principal material system. Solution This example is solved using MATLAB.First the strains are calculated at the five interface locations using the MATLAB function Strains as follows:
120 7 Laminate Analysis – Part I Strains(eps xo, eps yo, gam xyo, kap xo, kap yo, kap xyo, z) – This function calculates the three strains εx, εy, and γxy at any point P on the normal line given the three strains ε0 x, ε0 y, γ0 xy and the three curvatures κ0 x, κ0 y, κ0 xy at point P0, and the distance z between P and P0. There are seven input arguments to this function. The function returns the 3 × 1 strain vector. The following is a listing of the MATLAB source code for this function: function y = Strains(eps_xo,eps_yo,gam_xyo,kap_xo,kap_yo,kap_xyo,z) %Strains This function returns the strain vector at any point P % along the normal line at distance z from point Po which % lies on the reference surface. There are seven input % arguments for this function - namely the three strains % and three curvatures at point Po and the distance z. % The size of the strain vector is 3 x 1. epsilonx = eps_xo+z* kap_xo; epsilony = eps_yo+z* kap_yo; gammaxy = gam_xyo+z* kap_xyo; y = [epsilonx ; epsilony ; gammaxy]; MATLAB Example 7.1 Consider a graphite-reinforced polymer composite laminate with the elastic constants as given in Example 2.2. The laminate has total thickness of 0.500 mm and is stacked as a [0/90]S laminate. The four layers are of equal thickness. It is deformed so that at a point (x, y) on the reference surface, we have the following strains and curvatures: ε 0 x = 400 × 10−6 ε 0 y = γ0 xy = κ0 x = κ0 y = κ0 xy = 0 Use MATLAB to determine the following: (a) the three components of strain at the interface locations. (b) the three components of stress in each layer. Plot the stress distribution along the depth of the laminate for each component. (c) the force and moment resultants in the laminate. (d) the three components of strain at the interface locations with respect to the principal material system. (e) the three components of stress in each layer with respect to the principal material system. Solution This example is solved using MATLAB. First the strains are calculated at the five interface locations using the MATLAB function Strains as follows:
7.2 MATLAB Functions Used 121 >>eps11on1=Strains(400e-6,0,0,0,0,0,-0.250e-3) epsilon1 1.0e-003* 0.4000 0 0 >>epsi1on2=Strains(400e-6,0,0,0,0,0,-0.125e-3) epsilon2 1.0e-003* 0.4000 0 0 >>epsi1on3=Strains(400e-6,0,0,0,0,0,0) epsilon3 1.0e-003* 0.4000 0 0 >>epsi1on4=Strains(400e-6,0,0,0,0,0,0.125e-3) epsilon4 1.0e-003* 0.4000 0 0 >>epsi1on5=Strains(400e-6,0,0,0,0,0,0.250e-3) epsilon5 1.0e-003* 0.4000 0 0
7.2 MATLAB Functions Used 121 >> epsilon1 = Strains(400e-6,0,0,0,0,0,-0.250e-3) epsilon1 = 1.0e-003 * 0.4000 0 0 >> epsilon2 = Strains(400e-6,0,0,0,0,0,-0.125e-3) epsilon2 = 1.0e-003 * 0.4000 0 0 >> epsilon3 = Strains(400e-6,0,0,0,0,0,0) epsilon3 = 1.0e-003 * 0.4000 0 0 >> epsilon4 = Strains(400e-6,0,0,0,0,0,0.125e-3) epsilon4 = 1.0e-003 * 0.4000 0 0 >> epsilon5 = Strains(400e-6,0,0,0,0,0,0.250e-3) epsilon5 = 1.0e-003 * 0.4000 0 0
122 7 Laminate Analysis-Part I Next,the reduced stiffness [Q}in GPa is calculated for this material using the MATLAB function ReducedStiffness as follows: >>Q=ReducedStiffness(155.0,12.10,0.248,4.40) Q= 155.7478 3.0153 0 3.0153 12.1584 0 0 0 4.4000 The transformed reduced stiffnesses [Q]in GPa for the four layers are now cal- culated using the MATLAB function Qbar as follows: >Qbar1 =Qbar(Q,0) Qbar1 155.7478 3.0153 0 3.0153 12.1584 0 0 0 4.4000 >Qbar2 Qbar(Q,90) Qbar2 12.1584 3.0153 -0.0000 3.0153 155.7478 0.0000 -0.0000 0.0000 4.4000 >Qbar3 Qbar(Q,90) Qbar3 12.1584 3.0153 -0.0000 3.0153 155.7478 0.0000 -0.0000 0.0000 4.4000 >Qbar4 Qbar(Q,0) Qbar4 155.7478 3.0153 0 3.0153 12.1584 0 0 0 4.4000 Next,the stresses in each layer are calculated in MPa.Note that the stress vector is calculated twice for each layer-once at the top of the layer and once at the bottom of the layer
122 7 Laminate Analysis – Part I Next, the reduced stiffness [Q} in GPa is calculated for this material using the MATLAB function ReducedStiffness as follows: >> Q = ReducedStiffness(155.0, 12.10, 0.248, 4.40) Q = 155.7478 3.0153 0 3.0153 12.1584 0 0 0 4.4000 The transformed reduced stiffnesses [Q¯] in GPa for the four layers are now calculated using the MATLAB function Qbar as follows: >> Qbar1 = Qbar(Q,0) Qbar1 = 155.7478 3.0153 0 3.0153 12.1584 0 0 0 4.4000 >> Qbar2 = Qbar(Q,90) Qbar2 = 12.1584 3.0153 -0.0000 3.0153 155.7478 0.0000 -0.0000 0.0000 4.4000 >> Qbar3 = Qbar(Q,90) Qbar3 = 12.1584 3.0153 -0.0000 3.0153 155.7478 0.0000 -0.0000 0.0000 4.4000 >> Qbar4 = Qbar(Q,0) Qbar4 = 155.7478 3.0153 0 3.0153 12.1584 0 0 0 4.4000 Next, the stresses in each layer are calculated in MPa. Note that the stress vector is calculated twice for each layer – once at the top of the layer and once at the bottom of the layer
7.2 MATLAB Functions Used 123 >sigmala Qbar1*epsilon1*1e3 sigmala 62.2991 1.2061 0 >sigmaib Qbar1*epsilon2*1e3 sigmalb 62.2991 1.2061 0 >sigma2a =Qbar2*epsilon2*1e3 sigma2a 4.8634 1.2061 -0.0000 >sigma2b Qbar2*epsilon3*1e3 sigma2b 4.8634 1.2061 -0.0000 >sigma3a Qbar3*epsilon3*1e3 sigma3a 4.8634 1.2061 -0.0000 >sigma3b =Qbar3*epsilon4*1e3 sigma3b 4.8634 1.2061 -0.0000
7.2 MATLAB Functions Used 123 >> sigma1a = Qbar1*epsilon1*1e3 sigma1a = 62.2991 1.2061 0 >> sigma1b = Qbar1*epsilon2*1e3 sigma1b = 62.2991 1.2061 0 >> sigma2a = Qbar2*epsilon2*1e3 sigma2a = 4.8634 1.2061 -0.0000 >> sigma2b = Qbar2*epsilon3*1e3 sigma2b = 4.8634 1.2061 -0.0000 >> sigma3a = Qbar3*epsilon3*1e3 sigma3a = 4.8634 1.2061 -0.0000 >> sigma3b = Qbar3*epsilon4*1e3 sigma3b = 4.8634 1.2061 -0.0000
124 7 Laminate Analysis-Part I >sigma4a =Qbar4*epsilon4*1e3 sigma4a 62.2991 1.2061 0 >sigma4b Qbar4*epsilon5*1e3 sigma4b 62.2991 1.2061 0 Next,we setup the y-axis for the three plots: >y=[0.2500.1250.12500-0.125-0.125-0.250] y= 0.2500 0.1250 0.125000-0.1250 -0.1250 -0.2500 The distribution of the stress o along the depth of the laminate is now plotted as follows (see Fig.7.4): >>x [sigma4b(1)sigma4a(1)sigma3b(1)sigma3a(1)sigma2b(1) sigma2a(1)sigma1b(1)sigmala(1)] X= 62.2991 62.29914.8634 4.8634 4.8634 4.8634 62.2991 62.2991 >plot(x,y) >xlabel('\sigma_x (MPa)') >ylabel(‘z(mm)') The distribution of the stress ou along the depth of the laminate is now plotted as follows (see Fig.7.5): >x [sigma4b(2)sigma4a(2)sigma3b(2)sigma3a(2)sigma2b(2) sigma2a(2)sigma1b(2)sigmala(2)] X= 1.2061 1.20611.20611.20611.20611.20611.2061 1.2061 >plot(x,y) >ylabel(‘z(mm)') >xlabel('\sigma_y (MPa)')
124 7 Laminate Analysis – Part I >> sigma4a = Qbar4*epsilon4*1e3 sigma4a = 62.2991 1.2061 0 >> sigma4b = Qbar4*epsilon5*1e3 sigma4b = 62.2991 1.2061 0 Next, we setup the y-axis for the three plots: >> y = [0.250 0.125 0.125 0 0 -0.125 -0.125 -0.250] y = 0.2500 0.1250 0.1250 0 0 -0.1250 -0.1250 -0.2500 The distribution of the stress σx along the depth of the laminate is now plotted as follows (see Fig. 7.4): >> x = [sigma4b(1) sigma4a(1) sigma3b(1) sigma3a(1) sigma2b(1) sigma2a(1) sigma1b(1) sigma1a(1)] x = 62.2991 62.2991 4.8634 4.8634 4.8634 4.8634 62.2991 62.2991 >> plot(x,y) >> xlabel(‘\sigma_x (MPa)’) >> ylabel(‘z (mm)’) The distribution of the stress σy along the depth of the laminate is now plotted as follows (see Fig. 7.5): >> x = [sigma4b(2) sigma4a(2) sigma3b(2) sigma3a(2) sigma2b(2) sigma2a(2) sigma1b(2) sigma1a(2)] x = 1.2061 1.2061 1.2061 1.2061 1.2061 1.2061 1.2061 1.2061 >> plot(x,y) >> ylabel(‘z (mm)’) >> xlabel(‘\sigma_y (MPa)’)