See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/315741785

FEA FOR A 2D METAL PLATE USING MATLAB

Method · March 2017

DOI: 10.13140/RG.2.2.14711.42402

CITATIONS READS

0 1,050

1 author:

Hamad Khaled Alkhareef

University College London
8 PUBLICATIONS   0 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Open for new assignment in reasearch View project

All content following this page was uploaded by Hamad Khaled Alkhareef on 02 April 2017.

The user has requested enhancement of the downloaded file.

FINITE ELEMENT ANALYSIS OF A 2D METAL
PLATE USING MATLAB
HAMAD KH ALKHAREEF - DEPARTMENT OF MECHANICAL ENGINEERING (MSC)
UNIVERSITY COLLEGE LONDON
FEBRUARY 20, 2017

Abstract
This paper investigates the accuracy of Finite Element Analysis using commercial engineering
applications, Matlab and ANSYS, to compare result. Programming language for a two-
dimensional linear elastic metal plate made of Aluminium alloy consists of two and four; 3-
noded and 6-noded triangular elements under uniform tension was used for analysing
displacement and strain. Three additional Young’s modulus were examined, 148 GPa,
200 GPa, and 300 GPa. Also, three different load conditions were applied,  = 30o,  = 60o,
and  = 110o. The aim of this study is to examine the accuracy of displacement and strain
obtained from the FEM from the code developed and validate it with ANSYS result.

1. Introduction
Structures such as bridges and beams generally consist of number of elements which is
assembled together. Typically, bridges and beams are supported at one end and free at the
other end to allow material deformation resulting from thermal expansion and high loading.
Therefore, it is essential to carry out stress analysis to define the critical points of the structure
design and restrict it with boundary conditions to avoid failure. Also, similar test is required for
the selected material to study the behaviour of it which provides the maximum allowable stress
that can be applied on the material and hence its maximum deformation. These tests can be
carried out experimentally but its required an enormous amount of effort and money
particularly for complicated designs. However, due to the development of technology, FEM
can be used through different commercial engineering applications such as ANSYS, which
allows the user to analyse all the points of interest on the model with less effort and much
cheaper comparing to experimental test. FEM was established to study the stresses in
complex structures by solving large number of simultaneous equations.
In this report, Matlab R2016 and ANSYS v15 mechanical were used as a tool to compute the
Finite Element Approach to analyse displacement and strain of a rectangular metal plate under
uniform tensile stress with a uniform thickness. The plate is simply supported at two nodes as
shown in figure 1 and consists of two and four triangular elements with 3 and 6 nodes. The
purpose of this study is to compare the displacement and strain of the plate for various
elements, material properties, and load conditions. The accuracy of the code will be validated
using ANSYS result. Model description, material properties, boundary conditions, and element
type/shape will be illustrated in the report as well as the equations used. In addition, results
obtained from the code and ANSYS will be evaluated and relationship between strain and
elastic modulus will be achieved.

1
2. Methodology:
a. Description of the model
The model used consist of a uniform rectangular plate that has a height of
0.03 m and a width of 0.05 m with a thickness of 0.002 m.

b. Material selection
Isotropic Aluminium alloy with elastic modulus of 100 GPa and Poisson’s ratio
of 0.3 was initially implemented. In addition, different elastic materials have
been tested in this study such as Copper alloy with E = 148 GPa, Low Carbon
Steel with E = 200 GPa, and Silicon Carbide with E = 300 GPa.

c. Boundary conditions and load

The plate is fixed supported (UX=UY=0) at the bottom point on the left-hand
side and roller support was used for the top point, which permit the movement
in the y-direction (UX=0). On the other hand, tensile stress of 200 MPa is
applied uniformly in the x-direction at 90o as represented in figure 1. The
direction of load has been altered three different times from 90 o to  = 30o,  =
60o, and  = 110o with keeping other variables constant.

d. Element type
3-noded and 6-noded triangular elements will be used in this study. In order to
model the plate in ANSYS, 2-D structural solid mass quad 4-node PLANE 182
with triangular element shape was used for the 3-noded element method. This
element type can be used as either plane stress/strain or axisymmetric
element. Each node has two DOF, x and y. It has plasticity, large deflection,
and strain capability. PLANE 183 quad 8 was used for modelling the 6-noded
elements because of its high order nodes. This element type has quadratic
displacement behaviour and typically used for modelling irregular meshes. Like
PLANE 182, quad 8 has two DOF at each node. Figure 1 part 2 shows the 3
and 6 noded elements order as it has been done in the code and figure 2
represents the boundary conditions and load applied in ANSYS.

e. Coding and simulation

Programming language requires from the user to define all the model variables
such as dimensions, loads, and material properties. After this the FEM can be
computed using the same language for both 3-node and 6-node elements. The
level of complexity of the code depends on the designer skills. It is
recommended to use Matlab for this study. Also, the same approach has been
applied on ANSYS to examine the validity of the code.

f. Equations and measurements

• 3-noded triangular element:
The plate coordinates are plotted as shown in figure 1. In order to calculate the
strain and displacement of the element it is important to obtain the stiffness
matrix K for each element. This can be calculated by the following equation:
[𝑘] 𝑇 [𝐷][𝐵] (1)
= 𝑡𝐴[𝐵]

2
where t is the thickness of the model and A is the area of each element which
can be defined as:
1
𝐴= ∗ 𝑏𝑎𝑠𝑒 ∗ ℎ𝑒𝑖𝑔ℎ𝑡 (2)
2

and the strain-displacement matrix B is constructed as follows:

1 𝛽1 0 𝛽2 0 𝛽3 0
[𝐵] = [ 0 𝛾1 0 𝛾2 0 𝛾3 ] (3)
2𝐴
𝛾1 𝛽1 𝛾2 𝛽2 𝛾3 𝛽3
where 𝛽 and 𝛾 are the shape function coefficients, and can be calculated as
shown in equation 4.
𝛽1 = 𝑦2 − 𝑦3 𝛾1 = 𝑥2 − 𝑥3
𝛽2 = 𝑦3 − 𝑦1 𝛾2 = 𝑥3 − 𝑥1 (4)
𝛽3 = 𝑦1 − 𝑦2 𝛾3 = 𝑥2 − 𝑥1

the stress-strain matrix D:

1 𝑣 0
𝐸 𝑣 1 0 (5)
[𝐷] = [ 1 − 2𝑣 ]
1 − 𝑣2
0 0
2

where v is the Poisson’s ratio of the material. Therefore, the displacement

was obtained from:
[𝐹] = [𝑘][𝑈] (6)

the [K] matrix has been reduced due to the boundary conditions using
Gaussian elimination and F is the force matrix which contains the total force
applied on the plate and can be calculated from this relationship:

 = 𝐹/𝐴
where (A = height x thickness) the area of the plate in the z direction. After
obtaining the total force applied, it will be divided by the number of nodes in
which the load is acting. Hence, the strain can be calculated as follows:

[𝜀 ] = [𝐵 ][𝑈] (7)

• 6-noded triangular element:

The stiffness matrix of the quadratic triangular element requires the calculation
of three [B] and [K] matrices. The method used to calculate the displacement
and strain in the 6-noded element will be as follows:
1 1−𝐿1
[𝐾]𝑒 = ∫ [𝐵]𝑇 [𝐷][𝐵]𝑑𝑣 = ∫ ∫ [𝐵]𝑇 [𝐷][𝐵]𝑡|𝐽|𝑑𝐿1 𝑑𝐿2 (8)
𝑉 0 −0

3
or [𝐾]𝑒 = 𝐻. 𝑡. [𝐵]𝑇 [𝐷][𝐵]det[𝐽], where H is the weighting factor = 1/6

where [D] can be calculated from equation 5, and the strain-displacement

matrix [B]3x18 is calculated from:
[𝐵]3𝑥18 = [𝐵1 ]3𝑥4 [𝐵2 ]4𝑥4 [𝐵3 ]4𝑥18
(9)

where [B1], [B2], and [B3] are obtained using the derivative shape functions
matrix and Jacobian matrix:

1 0 0 0
[𝐵1 ] = [0 0 0 1]
0 1 1 0

𝐽𝐽22 𝐽𝐽12
− 0 0
|𝐽| |𝐽|
𝐽𝐽21 𝐽𝐽11
− 0 0
|𝐽| |𝐽|
[𝐵2 ] =
𝐽𝐽22 𝐽𝐽12
0 0 −
|𝐽| |𝐽|
𝐽𝐽21 𝐽𝐽11
0 0 −
[ |𝐽| |𝐽| ]

𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6

0 0 0 0 0 0
𝜕𝐿1 𝜕𝐿1 𝜕𝐿1 𝜕𝐿1 𝜕𝐿1 𝜕𝐿1
𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6
0 0 0 0 0 0
𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2
[𝐵3 ] =
𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6
0 0 0 0 0 0
𝜕𝐿1 𝜕𝐿1 𝜕𝐿1 𝜕𝐿1 𝜕𝐿1 𝜕𝐿1
𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6
0 0 0 0 0 0
[ 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 𝜕𝐿2 ]

And the quadratic shape functions in natural coordinates are given by:
𝑁4 = 4𝐿1 𝐿2
𝑁1 = 𝐿1 (2𝐿1 − 1)
𝑁5 = 4𝐿2 (1 − 𝐿1 − 𝐿2 )
𝑁2 = 𝐿2 (2𝐿2 − 1)
𝑁6 = 4𝐿1 (1 − 𝐿1 − 𝐿2 )
𝑁3 = (1 − 𝐿1 − 𝐿2 )(1 − 2𝐿1 − 2𝐿2 )

4
 derivatives of shape functions are given by:

𝜕𝑁1 𝜕𝑁1
= 4𝐿1 − 1 =0
𝜕𝐿1 𝜕𝐿2
𝜕𝑁2 𝜕𝑁2
=0 = 4𝐿2 − 1
𝜕𝐿1 𝜕𝐿2
𝜕𝑁3 𝜕𝑁3
= −3 + 4𝐿1 + 4𝐿2 = −3 + 4𝐿1 + 4𝐿2
𝜕𝐿1 𝜕𝐿2
𝜕𝑁4 𝜕𝑁4
= 4𝐿2 = 4𝐿1
𝜕𝐿1 𝜕𝐿2
𝜕𝑁5 𝜕𝑁5
= −4𝐿2 = 4 − 4𝐿1 − 8𝐿2
𝜕𝐿1 𝜕𝐿2
𝜕𝑁6 𝜕𝑁6
= 4 − 4𝐿2 − 8𝐿1 = −4𝐿1
𝜕𝐿1 𝜕𝐿2

The Jacobian for the isoparametric element is:

6 6
𝜕𝑁𝑖 𝜕𝑁𝑖
∑ 𝑥 ∑ 𝑦
𝜕𝐿1 𝑖 𝜕𝐿1 𝑖
𝑖=1 𝑖=1
[𝐽] = 6 6
𝜕𝑁𝑖 𝜕𝑁𝑖
∑ 𝑥𝑖 ∑ 𝑦
[ 𝑖=1 𝜕𝐿2 𝜕𝐿2 𝑖 ]
𝑖=1

Therefore, the displacement and strain can be calculated from:

𝑒
{𝜀}3𝑥1 = [𝐵]3𝑥18 {𝛿}18𝑥1 (10)

And the displacement is found from the relationship in equation 6. However,

the forces will be distributed on three nodes and will be, 1/6 the total force for
the top and bottom nodes and 2/3 for the middle node.

g. Validity of the program

The code was validated by comparing the displacement and strain obtained
with ANSYS result. The plate was modelled in ANSYS and the boundary
conditions and load from the program were applied in a similar manner on
ANSYS. The result was fairly satisfying; however, the order of the nodes is
different in ANSYS than the program as well as the decimal places.

3. Results
The stiffness matrices obtained from the program for one element for both, 3-noded and 6-
noded triangular, are shown in figure 3.

The code has been extended to solve the displacement and strain of two triangular elements.
Table 1 represents the comparison of displacement and strain obtained from the code with
ANSYS result for the two and four 3-noded triangular elements. In addition, table 2 compares
the displacement and strain obtained from the code with ANSYS for the two and four 6-noded
triangular elements. The deformation shape of the plate for both cases are shown in Figure 4.

5
Moreover, the code has been tested for different materials such as Copper alloy (E=148 GPa),
Low Carbon steel (E=200 GPa), and Silicon Carbide (E=300 GPa). The calculated
displacement and strain of these materials based on the four 6-noded triangular elements are
summarised in figure 5. Furthermore, displacement and strain have been computed when the
angle of the load applied is altered i.e.  = 30o,  = 60o, and  = 110o, the effect of orientated
load is shown in figure 5.

4. Discussion
Figure 3 part (a) represents the stiffness matrix for a 3-noded triangular element and part (b)
for the 6-noded element. These matrices were expanded to the global coordinate in order to
find the displacements, forces and strains of the element.
The result listed in table 1 in accordance to the global coordinate arrangement in figure 1
and2, proves the validity of the code for the two and four 3-noded elements. However, ANSYS
result were more accurate than the code, for instance, ANSYS provides more decimals in the
solution while the code gives zeros. It can be noticed that the 4 elements solution has more
nodes, hence additional degree of freedom which gives more accurate result than the 2
elements problem. In addition, the two and four 6-noded elements result for displacement and
strain are summarised in table 2. The values of displacement obtained from the code are
identical to ANSYS, which validates the accuracy of the code. Nevertheless, the strain values
have a high percentage of error, average of 110%, which dictates that the code has some
limitations when solving quadratic (6-noded) problems. Also, it is hard to find analytical form
of integrals, equation 8, and the jacobian calculation depends on the selection of shape
functions, which depends on the choice of elements. If the values of nodes added up for each
element, it will give a fairly close value comparing to ANSYS.
Furthermore, the four 6-noded triangular elements code was tested for 3 additional elastic
moduli as mentioned in the result section. The relationship between the elastic moduli of the
material and the strain is found to be inverse, low Young’s modulus has small strain and small
displacement, as shown in figure 5 (a). While the relationship between stress and strain is
proportional as shown in figure 5 (c). Moreover, the angle of the load applied was altered three
different times, and the effect of it was recognised and plotted in figure 5 (b). The maximum
force can be achieved on the X-direction when the load is applied at  = 90o, and hence
maximum strain. While the minimum can be obtained at  = 0o and  = 180o. On the other
hand, the maximum strain and displacement can be achieved at   0o to 15o for the Y-direction
forces and the minimum is at  = 180o.
Generally, the four 6-noded triangular elements have more accurate result than the 3-noded
(linear) problem due to the high number of simultaneous equations which allows the model to
displace in a flexible manner as shown in figure 4 (c) and (d).

5. Conclusion
The two and four 3-noded elements (linear) problem has less accurate result than the 6-noded
one. The code has some limitations when solving high order quadratic problems as shown in
table 2 for the strain. The strain is directly affected by the angle of the applied load as well as
the elastic moduli of the material used. High order elements provide accurate result and allow
flexibility when measuring the displacement or strain in a particular model.

6
References
[1] P. Seshu, Textbook of finite element analysis, 1st ed. New Delhi: Prentice-Hall of India, 2010.

[2] "4.182 PLANE182 2-D Structural Solid (UP19980821)", Ansys.stuba.sk. [Online]. Available:
http://www.ansys.stuba.sk/html/elem_55/chapter4/ES4-182.htm. [Accessed: 04- Feb- 2017].

[3] "PLANE183", Sharcnet.ca. [Online]. Available: https://www.sharcnet.ca/Software/Ansys/16.2.3/en-

us/help/ans_elem/Hlp_E_PLANE183.html. [Accessed: 04- Feb- 2017].

[4] Materials Data Book, 1st ed. Cambridge, UK: Cambridge University Engineering Department, 2003,
p. 11

7
Figures and tables

(a1) (c1)

0m 0.05 m 0m 0.05 m
(b1) (d1)
3 2 1
3 3 2 3 3 2
2
2 4
1
1 3
3
1 2 1 1 1
1 2 2
(a2)
(c2)
4 4 4
2 1 2 1 2 1
3 3 3
2 2 4
5 6 5 6 5 6
6 6 5 6 5
5
1 1 3
1 3 1 3 1 3
4 2 4 2 4 2
(b2) (d2)
Figure 1.Part 1 - (a1) shows the global nodes ordering for the 3-noded 2 elements triangular and (c1) is the 4 elements. (b1)
represents the global nodes of the 6-noded 2 elements triangular and (d1) is the 4 elements. Part 2 – (a2) shows the local
nodes arrangements of the 3-noded 2 elements and (c2) represents the 4 elements. (b2) and (d2) shows the 6-noded 2 and
4 elements local nodes arrangements.

8
 (a) (b)

(c) (d)
Figure 2. the boundary conditions and loads applied using ANSYS. (a) and (b) show the 3-noded 2 and 4 elements problem.
While (c) and (d) represent the 6-noded 2 and 4 elements problem. The number of nodes represent the global coordinate for
ANSYS solution.

(a)

(b)

Figure 3. (a) represents the stiffness matrix obtained for the 3-noded element and (b) is the [K] for the 6-noded element
model.

9
 3-noded Displacement 3-noded Strain
2 elements 4 elements 2 elements 4 elements
Node Matlab ANSYS Matlab ANSYS Element Matlab ANSYS Matlab ANSYS
u 0 0 0 0 εx 0.002 2.00E-03 0.002 2.00E-03
1
v 0 0 0 0 1 εy -0.0006 -6.00E-04 -0.0006 -6.00E-04
u 0 0 0 0 γxy 0.0000 5.96E-19 -0.0000 -1.35E-18
2
v -1.80E-05 -1.80E-05 -1.80E-05 -1.80E-05 εx 0.002 2.00E-03 0.002 2.00E-03
u 1.00E-04 1.00E-04 1.00E-04 5.00E-05 2 εy -0.0006 -6.00E-04 -0.0006 -6.00E-04
3
v -0.0000 -1.80E-05 -0.0000 -1.80E-05 γxy -0.0000 -2.10E-19 -0.0000 -2.10E-19
u 1.00E-04 1.00E-04 1.00E-04 5.00E-05 εx 0.002 2.00E-03
4
v -1.80E-05 5.08E-20 -1.80E-05 -4.27E-20 3 εy -0.0006 -6.00E-04
u 5.00E-05 1.00E-04 γxy -0.0000 -1.31E-18
5
v -0.0000 -1.80E-05 εx 0.002 2.00E-03
u 5.00E-05 1.00E-04 4 εy -0.0006 -6.00E-04
6
v -1.80E-05 -8.81E-20 γxy -0.0000 -1.00E-18

Table 1. comparison between Matlab result and ANSYS for displacement and strain obtained from the 3-noded; 2 and 4
elements.

6-noded Displacement 6-noded Strain

2 elements 4 elements
2 elements 4 elements
Node Matlab ANSYS Matlab ANSYS
u 0 0 0 0 Element Matlab ANSYS Matlab ANSYS
1 εx 0.005 3.37E-03 0.0055 5.11E-03
v 0 0 0 0
u 0 0 0 0 1 εy -0.0018 -1.98E-03 -0.0027 -3.59E-03
2
v -3.43E-05 -3.43E-05 -5.81E-05 -5.81E-05 γxy 0.0008 4.60E-03 0.0025 4.90E-03
u 1.26E-04 4.01E-05 1.40E-04 6.07E-05
3 εx 0.0065 2.02E-03 0.0067 2.23E-03
v -1.11E-05 -2.40E-05 -2.52E-05 -4.14E-05
u 1.26E-04 1.26E-04 1.40E-04 8.69E-05 2 εy -0.0018 -5.99E-04 -0.0016 -6.20E-04
4
v -2.38E-05 -2.38E-05 -4.14E-05 -3.90E-05 γxy -0.0008 -2.75E-05 -0.0025 -7.60E-04
u 7.37E-05 6.95E-05 5.43E-05 4.78E-05
5 εx 0.0057 1.95E-03
v -1.21E-05 -2.74E-05 -2.79E-05 -4.41E-05
u 1.27E-04 1.26E-04 4.78E-05 8.93E-05 3 εy -0.0012 -5.87E-04
6
v -1.62E-05 -1.11E-05 -4.41E-05 -3.17E-05 γxy 0.0003 1.69E-04
u 6.95E-05 1.27E-04 8.93E-05 9.30E-05
7 εx 0.0061 1.29E-03
v -2.74E-05 -1.62E-05 -3.17E-05 -3.24E-05
u 4.01E-05 7.37E-05 8.69E-05 5.43E-05 4 εy -0.0017 -7.14E-04
8
v -2.34E-05 -1.21E-05 -3.90E-05 -2.79E-05 γxy -0.0003 -3.69E-03
u 7.78E-05 7.78E-05 1.14E-04 7.30E-05
9
v -1.73E-05 -1.73E-05 -2.58E-05 -3.37E-05
u 1.14E-04 1.40E-04
10
v -4.14E-05 -4.14E-05
u 6.07E-05 1.14E-04
11
v -4.14E-05 -4.14E-05
u 7.30E-05 1.40E-04
12
v -3.36E-05 -2.52E-05
u 9.30E-05 1.14E-04
13
v -3.24E-05 -2.58E-05
u 1.16E-04 1.41E-04
14
v -3.29E-05 -3.29E-05
u 1.41E-04 1.16E-04
15
v -3.29E-05 -3.29E-05

Table 2. comparison between Matlab and ANSYS result of displacement and strain obtained for the 6-noded; 2 and 4
elements

10
 1

(a)

4
1

2
3

(b)

(c)

(d)
Figure 4. (a) shows the deformation shape of the plate from ANSYS and Matlab for the 3-noded; 2 elements and (b) for the
4 elements. while (c) and (d) show the deformation in the 6-noded; 2 elements and 4 elements respectively.

11
 The Relationship between Young's modulus and strain at (200 MPa) load
4.50E-03
4.00E-03
3.50E-03
Strain (unitless)

3.00E-03
2.50E-03
2.00E-03
1.50E-03
1.00E-03
5.00E-04
0.00E+00
0 50 100 150 200 250 300 350 400 450
Modulus of Elasticity (GPa)

The affect of theta of the applied force on the strain

0.0004
[top & bottom nodes] 0.0004
0.0003
0.0003 0.0002
0.0001
strain

0.0002 0
-0.0001
0.0001 -0.0002
-0.0003
0 -0.0004
0 15 30 45 60 75 90 105 120 135 150 165 180
strain x strain y angle (deg)

Relationship between stress-strain

2.50E-12

2.00E-12

1.50E-12
strain

1.00E-12

5.00E-13

0.00E+00
0 0.05 0.1 0.15 0.2 0.25
stress (GPa)

Figure 5. (a) represents the relationship between the elastic modulus and strain based on the dimension given. And (b)
shows the impact of the orientation of the load applied on the strain. And (c) represents the stress-strain relationship
(graph).

12
 k 2 = t * A* B2' * D* B2
Appendix k 2_gl obal = zer os ( 12, 12);
% % % % P ART 2 %%%( i ) 4- ELE ME NT S - 3- NODE D
k 2_gl obal ([ 3, 4, 9, 10, 11, 12] ,[ 3, 4, 9, 10, 11, 12] ) = ...
cl c;
k2 ([ 5, 6, 1, 2, 3, 4] ,[ 5, 6, 1, 2, 3, 4] )
E = 100e9; %Y OU N G' S MOD UL US
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
v = 0. 3; %P OI SS ON' S RATI O
%% % %% % %% % %% % %% %
t = 0. 002; %T HI CK NE S S ( m)
% %E L E ME NT 3
P = 200e6; %P RE S S URE APPLI ED ( Pa)
x 13=0. 025;
A = 0. 5* 0. 025* 0. 03; %A RE A ( m^ 2)
y 13=0;
f = ( P*t * 0. 03)/ 2; %F OR CE APPLI ED ( N)
x 23=0. 05;
D = ( E/ ( 1-( v^ 2)))*[ 1, v, 0 ; v, 1, 0 ; 0, 0, (( 1-
y 23=0;
v)/ 2)]; %ST RE S S- ST RAI N
x 33=0. 025;
%REL ATI ONS HI P
y 33=0. 03;
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
%% % %% % %% % %% % %% %
% %S HA P E F UNCTI ON COEF FI CI ENT S
% %E L E ME NT 1
ai 3 = ( x23* y 33) - ( x33* y 23);
x 11 =0;
aj 3 = ( x33* y 13) - ( x13* y 33);
y 11=0;
ak 3 = ( x13* y23) - ( x23* y13);
x 21=0. 025;
y 21=0;
bi 3 = y23- y33;
x 31=0;
bj 3 = y33- y13;
y 31=0. 03;
bk 3 = y13- y23;
% %S HA P E F UNCTI ON COEF FI CI ENT S
ci 3 = x33- x23;
ai 1 = ( x21* y 31) - ( x31* y 21);
cj 3 = x13- x33;
aj 1 = ( x31* y 11) - ( x11* y 31);
ck 3 = x23- x13;
ak 1 = ( x11* y21) - ( x21* y11);

% %[ J A C OBI AN] MAT RI X

bi 1 = y21- y31;
j 3 = [ cj 3, bj 3; ci 3, bi 3];
bj 1 = y31- y11;
J 3 = det (j 3);
bk 1 = y11- y21;

% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S

ci 1 = x31- x21;
MA T RI X
cj 1 = x11- x31;
B3 = ( 1/ J3) *[ bi 3, 0, bj 3, 0, bk3, 0; 0, ci 3, 0, cj 3, 0, ck3; ...
ck 1 = x21- x11;
ci 3, bi 3, cj 3, bj 3, ck3, bk 3]
% %[ J A C OBI AN] MAT RI X
k 3 = t * A* B3' * D* B3
j 1 = [ cj 1, bj 1; ci 1, bi 1];
k 3_gl obal = zer os ( 12, 12);
J 1 = det (j 1);
k 3_gl obal ([ 5, 6, 9, 10, 11, 12] ,[ 5, 6, 9, 10, 11, 12] ) = ...
k3 ([ 3, 4, 1, 2, 5, 6] , [ 3, 4, 1, 2, 5, 6] )
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
MA T RI X
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
B1 = ( 1/ J1) *[ bi 1, 0, bj 1, 0, bk1, 0; 0, ci 1, 0, cj 1, 0, ck1; ...
%% % %% % %% % %% % %% %
ci 1, bi 1, cj 1, bj 1, ck1, bk 1]
% %E L E ME NT 4
x 14=0. 05;
k 1 = t * A* B1' * D* B1
y 14=0;
k 1_gl obal = zer os ( 12, 12);
x 24=0. 05;
k 1_gl obal ([ 1, 2, 3, 4, 9, 10] ,[ 1, 2, 3, 4, 9, 10] ) = ...
y 24=0. 03;
k1([ 1, 2, 5, 6, 3, 4] ,[ 1, 2, 5, 6, 3, 4] )
x 34=0. 025;
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
y 34=0. 03;
%% % %% % %% % %% % %% %
%E L E ME NT 2
% %S HA P E F UNCTI ON COEF FI CI ENT S
x 12=0. 025;
ai 4 = ( x24* y 34) - ( x34* y 24);
y 12=0;
aj 4 = ( x34* y 14) - ( x13* y 34);
x 22=0. 025;
ak 4 = ( x14* y24) - ( x24* y14);
y 22=0. 03;
x 32=0;
bi 4 = y24- y34;
y 32=0. 03;
bj 4 = y34- y14;
bk 4 = y14- y24;
% %S HA P E F UNCTI ON COEF FI CI ENT S
ai 1 = ( x22* y 32) - ( x32* y 22);
ci 4 = x34- x24;
aj 1 = ( x32* y 12) - ( x12* y 32);
cj 4 = x14- x34;
ak 1 = ( x12* y22) - ( x22* y12);
ck 4 = x24- x14;
bi 2 = y22- y32;
% %[ J A C OBI AN] MAT RI X
bj 2 = y32- y12;
j 4 = [ cj 4, bj 4; ci 4, bi 4];
bk 2 = y12- y22;
J 4 = det (j 4);
ci 2 = x32- x22;
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
cj 2 = x12- x32;
MA T RI X
ck 2 = x22- x12;
B4 = ( 1/ J4) *[ bi 4, 0, bj 4, 0, bk4, 0; 0, ci 4, 0, cj 4, 0, ck4; ...
ci 4, bi 4, cj 4, bj 4, ck4, bk 4]
% %[ J A C OBI AN] MAT RI X
k 4 = t * A* B4' * D* B4
j 2 = [ cj 2, bj 2; ci 2, bi 2];
k 4_gl obal = zer os ( 12, 12);
J 2 = det (j 2);
k 4_gl obal ([ 5, 6, 7, 8, 11, 12] ,[ 5, 6, 7, 8, 11, 12] ) = ...
k4 ([ 1, 2, 3, 4, 5, 6] , [ 1, 2, 3, 4, 5, 6] )
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
MA T RI X
%% % %% % %% % %% % %% %
B2 = ( 1/ J2) *[ bi 2, 0, bj 2, 0, bk2, 0; 0, ci 2, 0, cj 2, 0, ck2; ...
% %S U M OF GL OB AL [ K]
ci 2, bi 2, cj 2, bj 2, ck2, bk 2]
K = k1_gl obal + k2_gl obal + k3_gl obal + k4_gl obal

13
 xi _yi 1 = [ x11 y11; x21 y21; x31 y31; x41 y41; x51 y51; x61
y 61] ;

% %A P PL Y B OUN DA RY CONDI TI ONS - [ u1=v 1=u2=0] % %L 1 and L2 VAL UE S

di sp(' u1=v 1=u2=0 ( due t o const r ai ns) THE REF ORE: ' ) L1 = 0. 5;
K_e = K ( 4: 12, 4: 12) L2 = 0. 5;

% %DI SPL A CE ME NT [ U] & ST RAI N [ S] CAL CUL ATI OON - % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
[ KQ= F] MA T RI X 1
F = [ 0; f; 0; f; 0; 0; 0; 0; 0]; %%F ORCE MAT RI X B = el e ment 1b_6 N( xi _yi 1, L1, L2);
U = i nv( K_e) * F; %%DI SPL A CE ME NT B11 = B;
u1 = 0; k = el e ment 1k _6 N( xi _yi 1, H, t, D, L1, L2);
v 1 = 0; k 11 = k;
u2 = 0;
v 2 = U( 1, 1); % %L 1 and L2 VAL UE S
u3 = U( 2, 1); L1 = 0. 5;
v 3 = U( 3, 1); L2 = 0;
u4 = U( 4, 1); % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
v 4 = U( 5, 1); MA T RI X 2
u5 = U( 6, 1); B = el e ment 1b_6 N( xi _yi 1, L1, L2);
v 5 = U( 7, 1); B21 = B;
u6 = U( 8, 1); k = el e ment 1k _6 N( xi _yi 1, H, t, D, L1, L2);
v 6 = U( 9, 1); k 21 = k;

di sp (' di spl ace me nt ' ); % %L 1 and L2 VAL UE S

di sp (' x y' ); L1 = 0;
U_ = [ u1 v1; u2 v2; u3 v3; u4 v4; u5 v5; u6 v6] L2 = 0. 5;
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
%% % %% % %% % %% % %% % MA T RI X 3
di sp (' ELE ME NT 1 : DI SPL A CE ME NT & F ORCE S' ) B = el e ment 1b_6 N( xi _yi 1, L1, L2);
U_ 1 = [ u1; v1; u3; v3; u2; v2] B31 = B;
f 1 = k1 * U_1 k = el e ment 1k _6 N( xi _yi 1, H, t, D, L1, L2);
di sp (' ELE ME NT 2 : DI SPL A CE ME NT & F ORCE S' ) k 31 = k;
U_ 2 = [ u3; v3; u4; v4; u2; v2]
f 2 = k2 * U_2 B1 = B11 + B21 + B31; %%[ B] MAT RI X EL E ME NT 1
di sp (' ELE ME NT 3 : DI SPL A CE ME NT & F ORCE S' ) K1 = k11 + k21 + k31; %%L OCAL [ K] EL E ME NT 1
U_ 3 = [ u3; v3; u5; v5; u4; v4] % %L OCAL T O GL OB AL [ K] MAT RI X 1
f 3 = k3 * U_3 K1_gl obal = zer os( 30);
di sp (' ELE ME NT 4 : DI SPL A CE ME NT & F ORCE S' ) K1_gl obal
U_ 4 = [ u5; v5; u6; v6; u4; v4] ([ 1, 2, 3, 4, 9, 10, 13, 14, 21, 22, 23, 24] ,[ 1, 2, 3, 4, 9, 10, 13, 14, 21, 22, 23,
f 4 = k4 * U_4 24] ) = ...
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % K1([ 1, 2, 5, 6, 7, 8, 3, 4, 11, 12, 9, 10] ,
%% % %% % %% % %% % %% % [ 1, 2, 5, 6, 7, 8, 3, 4, 11, 12, 9, 10] );
di sp (' ELE ME NT 1 : ST RAI N' ) %% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
S1 = B1* U_1 %% % %% % %% % %% % %% %
di sp (' ELE ME NT 2 : ST RAI N' ) % %E L E ME NT 2
S2 = B2* U_2 x 12 = 0. 025;
di sp (' ELE ME NT 3 : ST RAI N' ) y 12 = 0. 03;
S3 = B3* U_3 x 22 = 0;
di sp (' ELE ME NT 4 : ST RAI N' ) y 22 = 0. 03;
S4 = B4* U_4 x 32 = 0. 025;
% % % % % % % % P ART 2 %%%( i i ) 4- el e me nt s - 6- NODE D y 32 = 0;
cl c; x 42 = 0. 0125;
E = 100* 10^ 9; %Y OUN G' S MOD UL US y 42 = 0. 03;
v = 0. 3; %P OI SS ON' S RATI O x 52 = 0. 0125;
t = 0. 002; %T HI CK NE S S ( m) y 52 = 0. 015;
A = 0. 5* 0. 025* 0. 03; %A RE A ( m^ 2) x 62 = 0. 025;
P = 200* 10^ 6; %P RE S S URE APPLI ED ( Pa) y 62 = 0. 015;
f = ( P*t * 0. 03); %F ORCE APPLI ED ( N) % %[ C OOR DI NAT E S] MAT RI X E2
H = 1/ 6; %WEI GHTI NG F UNCTI ON xi _yi 2 = [ x12 y12; x22 y22; x32 y32; x42 y42; x52 y52; x62
D = ( E/ ( 1-( v^ 2)))*[ 1, v, 0 ; v, 1, 0 ; 0, 0, (( 1- y 62] ;
v)/ 2)]; %ST RE S S- ST RAI N
%REL ATI ONS HI P % %L 1 and L2 VAL UE S
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % L1 = 0. 5;
%% % %% % %% % %% % %% % L2 = 0. 5;
% %E L E ME NT 1 % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
x 11 = 0; MA T RI X 1
y 11 = 0; B = el e ment 2b_6 N( xi _yi 2, L1, L2);
x 21 = 0. 025; B12 = B;
y 21 = 0; k = el e ment 2k _6 N( xi _yi 2, H, t, D, L1, L2);
x 31 = 0; k 12 = k;
y 31 = 0. 03;
x 41 = 0. 0125; % %L 1 and L2 VAL UE S
y 41 = 0; L1 = 0. 5;
x 51 = 0. 0125; L2 = 0;
y 51 = 0. 015; % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
x 61 = 0; MA T RI X 2
y 61 = 0. 015; B = el e ment 2b_6 N( xi _yi 2, L1, L2);
% %[ C OOR DI NAT E S] MAT RI X E1 B22 = B;
k = el e ment 2k _6 N( xi _yi 2, H, t, D, L1, L2);

14
k 22 = k; %% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
%% % %% % %% % %% % %% %
% %L 1 and L2 VAL UE S % %E L E ME NT 4
L1 = 0; x 14 = 0. 05;
L2 = 0. 5; y 14 = 0. 03;
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S x 24 = 0. 025;
MA T RI X 3 y 24 = 0. 03;
B = el e ment 2b_6 N( xi _yi 2, L1, L2); x 34 = 0. 05;
B32 = B; y 34 = 0;
k = el e ment 2k _6 N( xi _yi 2, H, t, D, L1, L2); x 44 = 3* 0. 0125;
k 32 = k; y 44 = 0. 03;
x 54 = 3* 0. 0125;
B2 = B12 + B22 + B32; %%[ B] MAT RI X EL E ME NT 2 y 54 = 0. 015;
K2 = k12 + k22 + k32; %%L OCAL [ K] EL E ME NT 2 x 64 = 0. 05;
% %L OCAL T O GL OB AL [ K] MAT RI X 2 y 64 = 0. 015;
K2_gl obal = zer os( 30); % %[ C OOR DI NAT E S] MAT RI X E4
K2_gl obal xi _yi 4 = [ x14 y14; x24 y24; x34 y34; x44 y44; x54 y54; x64
([ 3, 4, 11, 12, 13, 14, 15, 16, 23, 24, 25, 26] ,[ 3, 4, 11, 12, 13, 14, 15, 16, 2 y 64] ;
3, 24, 25, 26] ) = ...
K2([ 3, 4, 7, 8, 5, 6, 1, 2, 9, 10, 11, 12] , % %L 1 and L2 VAL UE S
[ 3, 4, 7, 8, 5, 6, 1, 2, 9, 10, 11, 12] ); L1 = 0. 5;
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % L2 = 0. 5;
%% % %% % %% % %% % %% % % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
% %E L E ME NT 3 MA T RI X 1
x 13 = 0. 025; B = el e ment 4b_6 N( xi _yi 4, L1, L2);
y 13 = 0; B14 = B;
x 23 = 0. 05; k = el e ment 4k _6 N( xi _yi 4, H, t, D, L1, L2);
y 23 = 0; k 14 = k;
x 33 = 0. 025;
y 33 = 0. 03; % %L 1 and L2 VAL UE S
x 43 = 3* 0. 0125; L1 = 0. 5;
y 43 = 0; L2 = 0;
x 53 = 3* 0. 0125; % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
y 53 = 0. 015; MA T RI X 2
x 63 = 0. 025; B = el e ment 4b_6 N( xi _yi 4, L1, L2);
y 63 = 0. 015; B24 = B;
% %[ C OOR DI NAT E S] MAT RI X E3 k = el e ment 4k _6 N( xi _yi 4, H, t, D, L1, L2);
xi _yi 3 = [ x13 y13; x23 y23; x33 y33; x43 y43; x53 y53; x63 k 24 = k;
y 63] ;
% %L 1 and L2 VAL UE S
% %L 1 and L2 VAL UE S L1 = 0;
L1 = 0. 5; L2 = 0. 5;
L2 = 0. 5; % %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S MA T RI X 3
MA T RI X 1 B = el e ment 4b_6 N( xi _yi 4, L1, L2);
B = el e ment 3b_6 N( xi _yi 3, L1, L2); B34 = B;
B13 = B; k = el e ment 4k _6 N( xi _yi 4, H, t, D, L1, L2);
k = el e ment 3k _6 N( xi _yi 3, H, t, D, L1, L2); k 34 = k;
k 13 = k;
B4 = B14 + B24 + B34; %%[ B] MAT RI X EL E ME NT 4
% %L 1 and L2 VAL UE S K4 = k14 + k24 + k34; %%L OCAL [ K] EL E ME NT 4
L1 = 0. 5; % %L OCAL T O GL OB AL [ K] MAT RI X 4
L2 = 0; K4_gl obal = zer os( 30);
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S K4_gl obal
MA T RI X 2 ([ 5, 6, 7, 8, 15, 16, 19, 20, 27, 28, 29, 30] ,[ 5, 6, 7, 8, 15, 16, 19, 20, 27, 28,
B = el e ment 3b_6 N( xi _yi 3, L1, L2); 29, 30] ) = ...
B23 = B; K4([ 5, 6, 1, 2, 3, 4, 7, 8, 9, 10, 11, 12] ,
k = el e ment 3k _6 N( xi _yi 3, H, t, D, L1, L2); [ 5, 6, 1, 2, 3, 4, 7, 8, 9, 10, 11, 12] );
k 23 = k; %% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
%% % %% % %% % %% % %% %
% %L 1 and L2 VAL UE S % %S U M OF GL OB AL [ K]
L1 = 0; K_gl obal = K1_gl obal + K2_gl obal + K3_gl obal + K4_gl obal ;
L2 = 0. 5; % %A P PL Y B OUN DA RY CONDI TI ONS - [ u1=v 1=u2=0]
% %S T RAI N- DI SPL A CE ME NT REL ATI ON & STI FF NE S S di sp(' u1=v 1=u2=0 ( due t o const r ai ns) THE REF ORE: ' );
MA T RI X 3 K_e = K_gl obal ( 4: 30, 4: 30);
B = el e ment 3b_6 N( xi _yi 3, L1, L2); % %DI SPL A CE ME NT [ U] & ST RAI N [ S] CAL CUL ATI OON -
B33 = B; [ KQ= F]
k = el e ment 3k _6 N( xi _yi 3, H, t, D, L1, L2); F = [ 0; f/ 6; 0; f/ 6; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0;
k 33 = k; ...
0; 0; 0; (( 2/ 3) *f ); 0] ; %%F ORCE MAT RI X
B3 = B13 + B23 + B33; %%[ B] MAT RI X EL E ME NT 3 U = i nv( K_e) * F; %%DI SPL A CE ME NT
K3 = k13 + k23 + k33; %%L OCAL [ K] EL E ME NT 3
% %L OCAL T O GL OB AL [ K] MAT RI X 3 u1 = 0;
K3_gl obal = zer os( 30); v1 = 0;
K3_gl obal u2 = 0;
([ 5, 6, 13, 14, 15, 16, 17, 18, 25, 26, 27, 28] ,[ 5, 6, 13, 14, 15, 16, 17, 18, 2 v2 = U( 1, 1);
5, 26, 27, 28] ) = ... u3 = U( 2, 1);
K3([ 3, 4, 1, 2, 5, 6, 7, 8, 11, 12, 9, 10] , v3 = U( 3, 1);
[ 3, 4, 1, 2, 5, 6, 7, 8, 11, 12, 9, 10] ); u4 = U( 4, 1);
v4 = U( 5, 1);

15
 u5 = U( 6, 1); dN1dL2, 0, dN2dL2, 0, dN3dL2, 0, dN4dL2, 0, dN5dL2, 0,
v 5 = U( 7, 1); d N6dL2, 0;
u6 = U( 8, 1); 0, dN1dL1, 0, dN2dL1, 0, dN3dL1, 0, dN4dL1, 0, dN5dL1,
v 6 = U( 9, 1); 0, dN6dL1;
u7 = U( 10, 1); 0, dN1dL2, 0, dN2dL2, 0, dN3dL2, 0, dN4dL2, 0, dN5dL2,
v 7 = U( 11, 1); 0, dN6dL2] ;
u8 = U( 12, 1); %% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
v 8 = U( 13, 1); %% % %% % %
u9 = U( 14, 1); B = B1* B2* B3;
v 9 = U( 15, 1); end
u10 = U( 16, 1); %% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
v 10 = U( 17, 1); %% % %% % %% % %% % %% %
u11 = U( 18, 1); f uncti on k = el e ment 1k _6 N ( xi _yi 1, H, t, D, L1, L2)
v 11 = U( 19, 1); % % %Deri vati ves of t he Shape Funct i on
u12 = U( 20, 1); d N1dL1 = ( 4* L1) - 1;
v 12 = U( 21, 1); d N1dL2 = 0;
u13 = U( 22, 1); d N2dL1 = 0;
v 13 = U( 23, 1); d N2dL2 = ( 4* L2) - 1;
u14 = U( 24, 1); d N3dL1 = - 3 + ( 4* L1) + ( 4* L2);
v 14 = U( 25, 1); d N3dL2 = - 3 + ( 4* L1) + ( 4* L2);
u15 = U( 26, 1); d N4dL1 = ( 4* L2);
v 15 = U( 27, 1); d N4dL2 = ( 4* L1);
di sp (' DI SPL A CE ME NT' ); d N5dL1 = (- 4* L2);
di sp (' x y' ); d N5dL2 = 4 - ( 4* L1) - ( 8* L2);
U_ = [ u1 v1; u2 v2; u3 v3; u4 v4; u5 v5; u6 v6; u7 v7; u8 v8; d N6dL1 = 4 - ( 4* L2) - ( 8* L1);
u9 v9; ... d N6dL2 = (- 4* L1);
u10 v10; u11 v11; u12 v12; u13 v13; u14 v14; u15 v15] % % %S h a pe Funct i on Mat ri x
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % j = [ dN1dL1, dN2dL1, dN3dL1, dN4dL1, dN5dL1, dN6dL1;
%% % %% % %% % %% % %% % dN1dL2, dN2dL2, dN3dL2, dN4dL2, dN5dL2, dN6dL2] ;
di sp (' ELE ME NT 1 : DI SPL A CE ME NT & F ORCE S' ) % % %J ac obi an mat ri x ([ xi _yi ] mat ri x f r o m mai n scri pt )
U_ 1 = [ u1; v1; u7; v7; u2; v2; u5; v5; u12; v12; u11; v11] J = j * ( xi _yi 1);
f 1 = K1 * U_1 % % %St r ai n- Di spl ace me nt Mat ri x
di sp (' ELE ME NT 2 : DI SPL A CE ME NT & F ORCE S' ) B1 = [ 1, 0, 0, 0; 0, 0, 0, 1; 0, 1, 1, 0] ;
U_ 2 = [ u8; v8; u2; v2; u7; v7; u6; v6; u12; v12; u13; v13] B2 = [ J( 2, 2)/ det ( J), - J( 1, 2)/ det ( J), 0, 0;
f 2 = K2 * U_2 - J( 2, 1)/ det ( J), J( 1, 1)/ det ( J), 0, 0;
di sp (' ELE ME NT 3 : DI SPL A CE ME NT & F ORCE S' ) 0, 0, J( 2, 2)/ det ( J), - J( 1, 2)/ det ( J);
U_ 3 = [ u7; v7; u3; v3; u8; v8; u9; v9; u14; v14; u13; v13] 0, 0, - J( 2, 1)/ det ( J), J( 1, 1)/ det ( J)];
f 3 = K3 * U_3 B3 = [ dN1dL1, 0, dN2dL1, 0, dN3dL1, 0, dN4dL1, 0, dN5dL1,
di sp (' ELE ME NT 4 : DI SPL A CE ME NT & F ORCE S' ) 0, dN6dL1, 0;
U_ 4 = [ u4; v4; u8; v8; u3; v3; u10; v10; u14; v14; u15; v15] dN1dL2, 0, dN2dL2, 0, dN3dL2, 0, dN4dL2, 0, dN5dL2, 0,
f 4 = K4 * U_4 d N6dL2, 0;
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % 0, dN1dL1, 0, dN2dL1, 0, dN3dL1, 0, dN4dL1, 0, dN5dL1,
%% % %% % %% % %% % %% % 0, dN6dL1;
di sp (' ELE ME NT 1 : ST RAI N' ) 0, dN1dL2, 0, dN2dL2, 0, dN3dL2, 0, dN4dL2, 0, dN5dL2,
S1 = B1* U_1 0, dN6dL2] ;
di sp (' ELE ME NT 2 : ST RAI N' ) %% % %% % %% % %% % %% % %% % %% % %% % %% % %% % %
S2 = B2* U_2 %% % %% % %
di sp (' ELE ME NT 3 : ST RAI N' ) B = B1* B2* B3;
S3 = B3* U_3 k = H*t * B' * D* B* det ( J); %St i ff ness Mat ri x
di sp (' ELE ME NT 4 : ST RAI N' ) end
S4 = B4* U_4
%% % %% % %% % %% % %% % %% % %% % %% % %% % %% % % % % % %P L E A S E NOT E T HAT T HE S E T WO F UNCTI ONS
%% % %% % %% % %% % %% % HA V E BEE N USE D F OR EA CH EL E ME NT T O CAL CUL AT E
f uncti on B = el e ment 1b_6 N ( xi _yi 1, L1 , L2) [ B] AND [ K]
% % %Deri vati ves of t he Shape Funct i on
d N1dL1 = ( 4* L1) - 1;
d N1dL2 = 0;
d N2dL1 = 0;
d N2dL2 = ( 4* L2) - 1;
d N3dL1 = - 3 + ( 4* L1) + ( 4* L2);
d N3dL2 = - 3 + ( 4* L1) + ( 4* L2);
d N4dL1 = ( 4* L2);
d N4dL2 = ( 4* L1);
d N5dL1 = (- 4* L2);
d N5dL2 = 4 - ( 4* L1) - ( 8* L2);
d N6dL1 = 4 - ( 4* L2) - ( 8* L1);
d N6dL2 = (- 4* L1);
% % %S h a pe Funct i on Mat ri x
j = [ dN1dL1, dN2dL1, dN3dL1, dN4dL1, dN5dL1, dN6dL1;
dN1dL2, dN2dL2, dN3dL2, dN4dL2, dN5dL2, dN6dL2] ;
% % %J ac obi an mat ri x ([ xi _yi ] mat ri x f r o m mai n scri pt )
J = j * ( xi _yi 1);
% % %St r ai n- Di spl ace me nt Mat ri x
B1 = [ 1, 0, 0, 0; 0, 0, 0, 1; 0, 1, 1, 0] ;
B2 = [ J( 2, 2)/ det ( J), - J( 1, 2)/ det ( J), 0, 0;
- J( 2, 1)/ det ( J), J( 1, 1)/ det ( J), 0, 0;
0, 0, J( 2, 2)/ det ( J), - J( 1, 2)/ det ( J);
0, 0, - J( 2, 1)/ det ( J), J( 1, 1)/ det ( J)];
B3 = [ dN1dL1, 0, dN2dL1, 0, dN3dL1, 0, dN4dL1, 0, dN5dL1,
0, dN6dL1, 0;

16

View publication stats