Computational Fluid Dynamics Analysis of a Drone

     This post is about the computational fluid dynamics analysis of a small drone. The drone features a blended body-wing design with various cross sections at different span-wise locations. The drone has a wing-span and length of 6 ft. and 4.92 ft., respectively. The root (center) portion of the drone is relatively thicker and symmetrical in cross section for increased mechanical strength while the the mid-section and wing tips are thinner and utilize more cambered aero-foils. This is purely a concept design and as of now, no physical model of this drone exists.

     The numerical simulations for the present study were carried out using SolidWorks Flow Simulation Premium© code. The code employs κ-ε model with Two-Scales Wall Functions approach as the turbulence model. The numerical algorithm implemented is the SIMPLE-R, modified. The second-order upwind discretization scheme is used to approximate the convective fluxes while the diffusive terms are approximated using the central differencing scheme. The time derivatives are approximated with an implicit first order Euler scheme. The SolidWorks Flow Simulation© solves the Navier-Stokes equations, equations 1-3, which are formulations of mass, momentum and energy conservation laws for fluid flows. Turbulent flows are predicted using the Favre-averaged Navier-Stokes equations.

     The mesh independence test was carried out starting with 348,679 fluid cells. The mesh density was then increased up to 2,360,514 cells. The results of mesh independence study are mentioned below.

                          Mesh Name             Cells            Lift [N]         Drag [N]      Lift/Drag

                          M1                          348,679       382.41 -48.14      7.95

                          M2                          1,032,665    466.08      -48.73      9.57

                          M3                          1,559,516    473.48 -47.89      9.89

                          M4                          1,990,010    486.38 -48.08      10.12

                          M5                          2,360,514    491.07 -48.32      10.16

     It can be seen that as the mesh density increased, the difference in the critical parameters between two successive meshes also reduced. The mesh independence test was stopped as the difference between all of the critical parameters was less than one percent for the meshes M4 and M5.

     The pressure and velocity plots at various span-wise locations are shown in Fig. 1-2. It can be clearly seen that there is a negligible change in the velocity and pressure distributions around the drone between meshes M4 and M5. It can also be seen that as the mesh becomes finer, the resolution of both the pressure and velocity plots also increases.

Fig. 1, Velocity contours of various meshes.

Fig. 2, Pressure contours of various meshes.

     Aero-acoustics around the drone were also examined, as shown Fig. 3.

Fig. 3, Sound level contours of various meshes.

     A zoomed in view of the computational mesh is shown in Fig. 4. The refined mesh at the drone walls as a result of the mesh controls employed is clearly visible. The hump near the root of the drone is also visible, it was added in order to prevent the span-wise flow.

Fig. 4, Mesh level M4.

     The boundary conditions and the computational domain are shown in Fig. 5. The large red arrows represents inlet velocity boundary condition and the large blue arrows represents the atmospheric pressure outlet boundary condition. The red squares represents real wall boundary condition (slip) applied to the computational domain walls so that the boundary layer from the walls does not effect the flow around the drone.

An animation of an aileron roll can be seen here.

Thank you for reading. If you would like to collaborate on research projects, please reach out.

Automotive Computational Fluid Dynamics (CFD) Analysis

     This post is about the numerical analysis of an Ahmed body. It was a new experience because of the area between the car's floor and the road which is different from most of the numerical analysis performed in open channel aeronautics and turbo-machinery.

     The numerical analysis was performed using the commercial software, SolidWorks Flow Simulation. The software employs κ − ε turbulence model with damping functions, SIMPLE-R (modified), as the numerical algorithm and second-order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The time derivatives are approximated with an implicit first-order Euler scheme. Flow simulation solves the Navier–Stokes equations, which are formulations of mass, momentum, and energy conservation laws for fluid flows. To predict turbulent flows, the Favre-averaged Navier–Stokes equations are used.

     The software generates Cartesian mesh using immersed boundary method. The mesh had a cell size of 0.035 m in the far field regions within the computational domain. A fine mesh was need between the road the the car's floor to make sure the interaction of the car's floor with the road was captured accurately. Therefore the mesh between the car's floor and the road was refined to have a cell size of 0.00875 m. Another mesh control was applied around the body to refine the mesh with a cells size of 0.0175 m to capture the trailing vortices. The resulting mesh had 209,580 total cells, among those cells, 31,783 cells were at the solid fluid boundary. The computational domain size was ~1L x 1.12L x 3L where L being the vehicle's length. The computational domain along with the computational mesh is shown in Fig. 1.

Fig. 1 Mesh, computational domain and the boundary conditions.

     The red arrows within the Fig. 1 represents the inlet boundary condition of ambient (free-stream) velocity and the blue arrows represent the outlet boundary condition of the ambient pressure. The green arrows represents the co-ordinates axes direction.

     The results from the numerical analysis were compared with [1-3]. The results are within 10% of the experimental results. The velocity (superimposed by the velocity streamlines) and pressure profiles around the car body at various free-stream velocities is shown in Fig. 2.

Fig. 3 Velocity and pressure plots. From the top, Row 1, L-R; ambient velocity of 30 and 40 m/s. Row 2, L-R; ambient velocity of 60 and 80 m/s. Row 3, ambient velocity of 105 m/s.

     It was a good experience learning about automotive CFD after spending a long time in aeronautic/turbo-machinery CFD. Thank you for reading. Please share my work. If you would like to collaborate on a project please reach out.

[1] F.J.Bello-Millán, T.Mäkelä, L.Parras, C.delPino, C.Ferrera, "Experimental study on Ahmed's body drag coefficient for different yaw angles", Journal of Wind Engineering and Industrial Aerodynamics, Volume 157, October 2016, Pages 140-144.

[2] Guilmineau E., Deng G.B., Queutey P., Visonneau M. (2018) Assessment of Hybrid LES Formulations for Flow Simulation Around the Ahmed Body. In: Deville M. et al. (eds) Turbulence and Interactions. TI 2015. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol 135. Springer, Cham.

[3] A. Thacker, S.Aubrun, A.Leroy, P.Devinant, "Effects of suppressing the 3D separation on the rear slant on the flow structures around an Ahmed body", Journal of Wind Engineering and Industrial Aerodynamics, Volumes 107–108, August–September 2012, Pages 237-243.

Bond Graph Representation of the Coaxial Contra-Rotating Propeller

Abstract of my research in to the bond graph modeling for a coaxial contra-rotating propeller is given below.

A system of propellers formed when two propellers rotate about the same axis but in opposite direction is called coaxial contra-rotating propeller (Noriyuki Sasaki, 1998). A well-designed contra-rotating propeller has no rotational air flow, a maximum amount of air is pushed uniformly through the propeller disk resulting in an increase in performance and low induced energy losses. This study focuses on obtaining the system of equations needed to represent a contra-rotating propeller to be designed using seven spur gears mounted on five shafts. This study was completed using a bond graph model with a desire of obtaining a clearer understanding of this multi-domain , nonlinear and nonhomogeneous  system. A bond graph is a graphical representation of a physical dynamic system. The study resulted in a total of thirteen differential equations; needed to represent the system.

Finite Element Method: 1D Elasticity (MATLAB)

Code

%code solves 1D elasticity problems with up to 4 elements (5 nodes)

clc;

n=input('Press one for one, two for two, three for three and four for four elements in the system '); %elements in the system

clc;

if n==1

    A1=input('Enter the area of the element in mm^2 '); %area input

    E1=input('Enter the Young''s modulus for the element in MPa '); %Young's modulus input

    L1=input('Enter the length of the element in mm '); %length input

    DT=input('Enter the change in temperature in degree Celsius ');

    TA1=input('Enter the coefficient of thermal expansion for the element in per degree Celsius ');

    Px1i=input('Enter the pressure acting on the face at node i '); Px1j=input('Enter the pressure acting on the face at node j '); %pressure input

    X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

    P1i=input('Enter the nodal force on node i '); P1j=input('Enter the nodal force on node j '); %nodal force input

    clc;

    I=[1 -1;-1 1]; kg=((A1*E1)/L1)*I; %global stiffness matrix calculation

    FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; %thermal force vector calculation

    FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp1j=A1*Px1j*FPj.'; %pressure force vector

    FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; %body force vector

    FN1=[P1i P1j].'; %nodal force vector

    fg=ft1+fp1i+fp1j+bf1+FN1; %global force vector calculation

    u=input('Press one for u1=0, two for u2=0 '); %boundary conditions

    clc;

    if u==1

        fgu=zeros(1,1); fgu(1,1)=fg(2,1); %global force vector when u1=0

        kgu=zeros(1,1); kgu(1,1)=kg(2,2); %global stiffness matrix when u1=0

        U=kgu\fgu; %displacement calculation

        e1=(U(1,1)-0)/L1; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; %strain due to forces

        S1=ef1*E1; %stress calculation

        R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction force calculation

        uU=['Displacement at node j is ',num2str(U),' mm']; disp(uU); %displacement output

        rR=['Reaction force at node i is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element is ',num2str(ef1)]; disp(eE); %strain output

        sS=['Stress in the element is ',num2str(S1),' MPa']; disp(sS); %stress output

    end

    if u==2

        fgu=zeros(1,1); fgu(1,1)=fg(1,1); %global force vector when u1=0

        kgu=zeros(1,1); kgu(1,1)=kg(1,1); %global stiffness matrix when u1=0

        U=kgu\fgu; %displacement calculation

        e1=(0-U(1,1))/L1; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; %strain due to forces

        S1=ef1*E1; %stress calculation

        R1=(kg(2,1)*U(1,1))-fg(2,1); %reaction force calculation

        uU=['Displacement at node i is ',num2str(U),' mm']; disp(uU); %displacement output

        rR=['Reaction force at node j is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element is ',num2str(ef1)]; disp(eE); %strain output

        sS=['Stress in the element is ',num2str(S1),' MPa']; disp(sS); %stress output

    end

end

if n==2

    A1=input('Enter the area of element number one in mm^2 '); A2=input('Enter the area of element number two in mm^2 '); %area input

    E1=input('Enter the Young''s modulus for element number one in MPa '); E2=input('Enter the Young''s modulus for element number two in MPa '); %Young's modulus input

    L1=input('Enter the length of element one in mm '); L2=input('Enter the length of element two in mm '); %length input

    DT=input('Enter the change in temperature in degree Celsius ');

    TA1=input('Enter the coefficient of thermal expansion for element number one in per degree Celsius '); TA2=input('Enter the coefficient of thermal expansion for element number two in per degree Celsius '); %thermal expansion coefficient input

    Px1i=input('Enter the pressure acting on the face at node i for element number one '); Px2j=input('Enter the pressure acting on the face at node j for element number two '); %pressure input

    X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

    P1i=input('Enter the nodal force on node i for element number one '); P1j=input('Enter the nodal force on node j for element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for element number two '); %nodal force input

    clc;

    I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I; ke2=((A2*E2)/L2)*I; %element stiffness matrix calculation

    kg=zeros(n+1,n+1); %initialize global stiffness matrix

    kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2); kg(1,3)=0; kg(2,1)=ke1(2,1); kg(2,2)=ke1(1,1)+ke2(2,2); kg(2,3)=ke2(2,1); kg(3,1)=0; kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2); %global stiffness matrix formulation

    FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; ft2=E2*A2*TA2*DT*FT.'; %thermal force vector

    FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp2j=A2*Px2j*FPj.'; %pressure force vector

    FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; bf2=(X*A2*L2/2)*FB.'; %body force vector

    FN1=[P1i P1j].'; FN2=[P2i P2j].'; %nodal force vector

    f1=ft1+fp1i+bf1+FN1; f2=ft2+fp2j+bf2+FN2; %force vector calculation

    fg=zeros(3,1); %global force vector initialize

    fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1); fg(3,1)=f2(2,1); %global force vector calculation

    u=input('Press one for u1=0, two for u3=0, three for u1=u3=0 '); %boundary conditions

    clc;

    if u==1

        fgu=zeros(2,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); %global force vector when u1=0

        kgu=zeros(2,2); kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); %global stiffness matrix when u1=0

        U=kgu\fgu; %displacement calculation

        e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; %stress calculation

        R1=(kg(2,1)*U(1,1))-fg(1,1); %reaction force calculation

        uU=['Displacement at node two and three is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one and two is ',num2str(ef1),' and ',num2str(ef2)]; disp(eE); %strain output

        sS=['Stress in the element one and two is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output

    end

    if u==2

        fgu=zeros(2,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); %global force vector when u3=0

        kgu=zeros(2,2); kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); %global stiffness matrix when u3=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(3,2)*U(2,1))-fg(3,1); %reaction force calculation

        e1=(U(2,1)-U(1,1))/L1; e2=(0-U(2,1))/L2; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; %stress calculation

        uU=['Displacement at node one and two is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element two, node j is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one and two is ',num2str(ef1),' and ',num2str(ef2)]; disp(eE); %strain output

        sS=['Stress in the element one and two is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output

    end

    if u==3

        fgu=zeros(1,1); fgu(1,1)=fg(2,1); %global force vector when u3=0

        kgu=zeros(1,1); kgu(1,1)=kg(2,2); %global stiffness matrix when u1=u3=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(2,1)*U(1,1))-fg(1,1); %reaction force calculation

        R3=(kg(3,2)*U(1,1))-fg(3,1); %reaction force calculation

        e1=(U(1,1)-0)/L1; e2=(0-U(1,1))/L2; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2;%stress calculation

        uU=['Displacement at node two is ',num2str(U(1,1)),' mm']; disp(uU); %displacement output

        rR1=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR1); %reaction force output

        rR3=['Reaction force at element two, node j is ',num2str(R3),' N']; disp(rR3); %reaction force output

        eE=['Strain in the element one and two is ',num2str(ef1), ' and ',num2str(ef2)]; disp(eE); %strain output

        sS=['Stress in the element one and two is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output

    end

end

if n==3

    A1=input('Enter the area of element number one in mm^2 '); A2=input('Enter the area of element number two in mm^2 '); A3=input('Enter the area of element number three in mm^2 '); %area input

    E1=input('Enter the Young''s modulus for element number one in MPa '); E2=input('Enter the Young''s modulus for element number two in MPa '); E3=input('Enter the Young''s modulus for element number three in MPa '); %Young's modulus input

    L1=input('Enter the length of element one in mm '); L2=input('Enter the length of element two in mm '); L3=input('Enter the length of element three in mm '); %length input

    DT=input('Enter the change in temperature in degree Celsius ');

    TA1=input('Enter the coefficient of thermal expansion for element number one in per degree Celsius '); TA2=input('Enter the coefficient of thermal expansion for element number two in per degree Celsius '); TA3=input('Enter the coefficient of thermal expansion for element number three in per degree Celsius '); %thermal expansion coefficient input

    Px1i=input('Enter the pressure acting on the face at node i for element number one '); Px3j=input('Enter the pressure acting on the face at node j for element number three '); %pressure input

    X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

    P1i=input('Enter the nodal force on node i for element number one '); P1j=input('Enter the nodal force on node j for element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for element number two '); P3i=P2j; P3j=input('Enter the nodal force on node j for element number three '); %nodal force input

    clc;

    I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I; ke2=((A2*E2)/L2)*I; ke3=((A3*E3)/L3)*I; %element stiffness matrix calculation

    kg=zeros(n+1,n+1); %initialize global stiffness matrix

    kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2); kg(2,1)=ke1(2,1); kg(2,2)=ke1(2,2)+ke2(1,1); kg(2,3)=ke2(1,2); kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2)+ke3(1,1); kg(3,4)=ke3(1,2); kg(4,3)=ke3(2,1); kg(4,4)=ke3(2,2); %global stiffness matrix calculation

    FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; ft2=E2*A2*TA2*DT*FT.'; ft3=E3*A3*TA3*DT*FT.'; %thermal force vector

    FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp3j=A3*Px3j*FPj.'; %pressure force vector

    FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; bf2=(X*A2*L2/2)*FB.'; bf3=(X*A3*L3/2)*FB.'; %body force vector

    FN1=[P1i P1j].'; FN2=[P2i P2j].';  FN3=[P3i P3j].'; %nodal force vector

    f1=ft1+fp1i+bf1+FN1; f2=ft2+bf2+FN2; f3=ft3+fp3j+bf3+FN3; %force vector calculation

    fg=zeros(4,1); %global force vector initialize

    fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1); fg(3,1)=f2(2,1)+f3(1,1); fg(4,1)=f3(2,1); %global force vector formulation

    u=input('Press one for u1=0, two for u4=0, three for u1=u4=0 '); %boundary conditions

    clc;

    if u==1

        fgu=zeros(3,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); fgu(3,1)=fg(4,1); %global force vector when u1=0

        kgu=zeros(3,3); kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); %global stiffness matrix when u1=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction force calculation

        e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; %stress calculation

        uU=['Displacement at node two, three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm and ',num2str(U(3,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one, two and three is ',num2str(ef1),', ',num2str(ef2),' and ',num2str(ef3)]; disp(eE); %strain output

        sS=['Stress in the element one, two and three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa']; disp(sS); %stress output

    end

    if u==2

        fgu=zeros(3,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); %global force vector when u4=0

        kgu=zeros(3,3); kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); kgu(2,3)=kg(2,3); kgu(3,2)=kg(3,2); kgu(3,3)=kg(3,3); %global stiffness matrix when u4=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(4,3)*U(3,1))-fg(4,1); %reaction force calculation

        e1=(0-U(3,1))/L1; e2=(U(3,1)-U(2,1))/L2; e3=(U(2,1)-U(1,1))/L3; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; %stress calculation

        uU=['Displacement at node one, two and three is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm and ',num2str(U(3,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element three, node j is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one, two and three is ',num2str(ef1),', ',num2str(ef2),' and ',num2str(ef3)]; disp(eE); %strain output

        sS=['Stress in the element one, two and three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa']; disp(sS); %stress output

    end

    if u==3

        fgu=zeros(2,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); %global force vector when u1=u4=0

        kgu=zeros(2,2); kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); %global stiffness matrix when u4=u1=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(1,2)*U(1,1))-fg(1,1); R2=(kg(4,3)*U(2,1))-fg(4,1); %reaction force calculation

        e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(0-U(2,1))/L3; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; %stress calculation

        uU=['Displacement at node two and three is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element one, node i and element three, node j is ',num2str(R1),' N and ',num2str(R2),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one, two and three is ',num2str(ef1),' , ',num2str(ef2),' and ',num2str(ef3)]; disp(eE); %strain output

        sS=['Stress in the element one, two and three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa']; disp(sS); %stress output

    end

end

if n==4

    A1=input('Enter the area of element number one in mm^2 '); A2=input('Enter the area of element number two in mm^2 '); A3=input('Enter the area of element number three in mm^2 '); A4=input('Enter the area of element number four in mm^2 '); %area input

    E1=input('Enter the Young''s modulus for element number one in MPa '); E2=input('Enter the Young''s modulus for element number two in MPa '); E3=input('Enter the Young''s modulus for element number three in MPa '); E4=input('Enter the Young''s modulus for element number four in MPa '); %Young's modulus input

    L1=input('Enter the length of element one in mm '); L2=input('Enter the length of element two in mm '); L3=input('Enter the length of element three in mm '); L4=input('Enter the length of element four in mm '); %length input

    DT=input('Enter the change in temperature in degree Celsius ');

    TA1=input('Enter the coefficient of thermal expansion for element number one in per degree Celsius '); TA2=input('Enter the coefficient of thermal expansion for element number two in per degree Celsius '); TA3=input('Enter the coefficient of thermal expansion for element number three in per degree Celsius '); TA4=input('Enter the coefficient of thermal expansion for element number four in per degree Celsius '); %thermal expansion coefficient input

    Px1i=input('Enter the pressure acting on the face at node i for element number one '); Px4j=input('Enter the pressure acting on the face at node j for element number four '); %pressure input

    X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

    P1i=input('Enter the nodal force on node i for element number one '); P1j=input('Enter the nodal force on node j for element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for element number two '); P3i=P2j; P3j=input('Enter the nodal force on node j for element number three '); P4i=P3j; P4j=input('Enter the nodal force on node j for element number four '); %nodal force input

    clc;

    I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I; ke2=((A2*E2)/L2)*I; ke3=((A3*E3)/L3)*I; ke4=((A4*E4)/L4)*I; %element stiffness matrix calculation

    kg=zeros(n+1,n+1); %initialize global stiffness matrix

    kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2); kg(2,1)=ke1(2,1); kg(2,2)=ke1(2,2)+ke2(1,1); kg(2,3)=ke2(1,2); kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2)+ke3(1,1); kg(3,4)=ke3(1,2); kg(4,3)=ke3(2,1); kg(4,4)=ke3(2,2)+ke4(1,1); kg(4,5)=ke4(1,2); kg(5,4)=ke4(2,1); kg(5,5)=ke4(2,2); %global stiffness matrix calculation

    FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; ft2=E2*A2*TA2*DT*FT.'; ft3=E3*A3*TA3*DT*FT.'; ft4=E4*A4*TA4*DT*FT.'; %thermal force vector

    FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp4j=A4*Px4j*FPj.'; %pressure force vector

    FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; bf2=(X*A2*L2/2)*FB.'; bf3=(X*A3*L3/2)*FB.'; bf4=(X*A4*L4/2)*FB.'; %body force vector

    FN1=[P1i P1j].'; FN2=[P2i P2j].';  FN3=[P3i P3j].'; FN4=[P4i P4j].';%nodal force vector

    f1=ft1+fp1i+bf1+FN1; f2=ft2+bf2+FN2; f3=ft3+bf3+FN3; f4=ft4+fp4j+bf4+FN4; %force vector calculation

    fg=zeros(5,1); %global force vector initialize

    fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1); fg(3,1)=f2(2,1)+f3(1,1); fg(4,1)=f3(2,1)+f4(1,1); fg(5,1)=f4(2,1); %global force vector formulation

    u=input('Press one for u1=0, two for u5=0, three for u1=u5=0 '); %boundary conditions0

    clc;

    if u==1

        fgu=zeros(4,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); fgu(3,1)=fg(4,1); fgu(4,1)=fg(5,1); %global force vector when u1=0

        kgu=zeros(4,4);

        kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); kgu(3,4)=kg(4,5); kgu(4,3)=kg(5,4); kgu(4,4)=kg(5,5);%global stiffness matrix when u1=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction force calculation

        e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; e4=(U(4,1)-U(3,1))/L4; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; S4=ef4*E4; %stress calculation

        uU=['Displacement at node two, three, four and five is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),' mm and ',num2str(U(4,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one, two, three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and ',num2str(ef4)]; disp(eE); %strain output

        sS=['Stress in the element one, two, three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa and ',num2str(S4),' MPa']; disp(sS); %stress output

    end

    if u==2

        fgu=zeros(4,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); fgu(4,1)=fg(4,1); %global force vector when u5=0

        kgu=zeros(4,4);

        kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); kgu(2,3)=kg(2,3); kgu(3,2)=kg(3,2); kgu(3,3)=kg(3,3); kgu(3,4)=kg(3,4); kgu(4,3)=kg(4,3); kgu(4,4)=kg(4,4); %global stiffness matrix when u5=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(5,4)*U(4,1))-fg(5,1); %reaction force calculation

        e1=(U(2,1)-U(1,1))/L1; e2=(U(3,1)-U(2,1))/L2; e3=(U(4,1)-U(3,1))/L3; e4=(0-U(4,1))/L4; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; S4=ef4*E4; %stress calculation

        uU=['Displacement at node one, two, three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),' mm and ',num2str(U(4,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element four, node j is ',num2str(R1),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one, two, three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and ',num2str(ef4)]; disp(eE); %strain output

        sS=['Stress in the element one, two, three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa and ',num2str(S4),' MPa']; disp(sS); %stress output

    end

    if u==3

        fgu=zeros(3,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); %global force vector when u1=u5=0

        kgu=zeros(3,3);

        kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); %global stiffness matrix when u5=0

        U=kgu\fgu; %displacement calculation

        R1=(kg(5,4)*U(3,1))-fg(5,1); R2=(kg(1,2)*U(1,1))-fg(5,1); %reaction force calculation

        e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; e4=(0-U(3,1))/L4; %strain calculation

        e0=TA1*DT; %thermal strain

        ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0; %strain due to forces

        S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; S4=ef4*E4; %stress calculation

        uU=['Displacement at node two, three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),' mm']; disp(uU); %displacement output

        rR=['Reaction force at element one, node i and element four, node j is ',num2str(R1),' N and ',num2str(R2),' N']; disp(rR); %reaction force output

        eE=['Strain in the element one, two, three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and ',num2str(ef4)]; disp(eE); %strain output

        sS=['Stress in the element one, two, three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa and ',num2str(S4),' MPa']; disp(sS); %stress output

    end

end

Sample calculations for problem with one element

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 1

Enter the area of the element in mm^2 20

Enter the Young's modulus for the element in MPa 200000

Enter the length of the element in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for the element in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i 0

Enter the pressure acting on the face at node j 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i 0

Enter the nodal force on node j 10000

Press one for u1=0, two for u2=0 1

Result

Displacement at node j is 0.26 mm

Reaction force at node i is -10000 N

Strain in the element is 0.0025

Stress in the element is 500 MPa

With second boundary condition (u₂=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 1

Enter the area of the element in mm^2 20

Enter the Young's modulus for the element in MPa 200000

Enter the length of the element in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for the element in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i 0

Enter the pressure acting on the face at node j 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i -10000

Enter the nodal force on node j for element number one 0

Press one for u1=0, two for u2=0 2

Result

Displacement at node i is -0.26 mm

Reaction force at node j is 10000 N

Strain in the element is 0.0025

Stress in the element is 500 MPa

Sample calculations for problem with two elements

Problem 5.1a, Fagan pg. 94

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 10000

Press one for u1=0, two for u3=0, three for u1=u3=0 1

Result

Displacement at node two and three is 0.26 mm and 0.77 mm

Reaction force at element one, node i is -10000 N

Strain in the element one and two is 0.0025 and 0.005

Stress in the element one and two is 500 MPa and 1000 MPa

With second boundary condition (u₃=0)

Problem 5.1a, Fagan pg. 94 mirrored

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 10

Enter the area of element number two in mm^2 20

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one -10000

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Press one for u1=0, two for u3=0, three for u1=u3=0 2

Result

Displacement at node one and two is -0.77 mm and -0.26 mm

Reaction force at element two, node j is 10000 N

Strain in the element one and two is 0.005 and 0.0025

Stress in the element one and two is 1000 MPa and 500 MPa

With third boundary condition (u₃=u₁=0)

Problem 5.1b, Fagan pg. 94

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Press one for u1=0, two for u3=0, three for u1=u3=0 3

Result

Displacement at node two is 0.0033333 mm

Reaction force at element one, node i is 266.6667 N

Reaction force at element two, node j is -266.6667 N

Strain in the element one and two is -6.6667e-05 and -0.00013333

Stress in the element one and two is -13.3333 MPa and -26.6667 MPa

Sample calculations for problem with three elements

Problem 5.1f, Fagan pg. 95

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 3

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 20

Enter the area of element number three in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the change in temperature in degree Celsius 0

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number three 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 20000

Enter the nodal force on node j for element number two 0

Enter the nodal force on node j for element number three 10000

Press one for u1=0, two for u4=0, three for u1=u4=0 1

Result

Displacement at node two, three and four is 1.25 mm, 1.5 mm and 2 mm

Reaction force at element one, node i is -50000 N

Strain in the element one, two and three is 0.0125, 0.0025 and 0.005

Stress in the element one, two and three is 2500 MPa, 500 MPa and 1000 MPa

With second boundary condition (u₄=0)

Problem 5.1f, Fagan pg. 95 mirrored

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 3

Enter the area of element number one in mm^2 10

Enter the area of element number two in mm^2 20

Enter the area of element number three in mm^2 20

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the change in temperature in degree Celsius 0

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number three 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one -10000

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two -20000

Enter the nodal force on node j for element number three 0

Press one for u1=0, two for u4=0, three for u1=u4=0 2

Result

Displacement at node one, two and three is -2 mm, -1.5 mm and -1.25 mm

Reaction force at element three, node j is 50000 N

Strain in the element one, two and three is 0.0125, 0.0025 and 0.005

Stress in the element one, two and three is 2500 MPa, 500 MPa and 1000 MPa

With third boundary condition (u₄=u₁=0)

Problem 5.1f, Fagan pg. 95 with temperature load only, no nodal forces

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 3

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 20

Enter the area of element number three in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number three 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Enter the nodal force on node j for element number three 0

Press one for u1=0, two for u4=0, three for u1=u4=0 3

Result

Displacement at node two and three is 0.0025 mm and 0.005 mm

Reaction force at element one, node i and element three, node j is 300 N and -300 N

Strain in the element one, two and three is -7.5e-05 , -7.5e-05 and -0.00015

Stress in the element one, two and three is -15 MPa, -15 MPa and -30 MPa

Sample calculations for problem with four elements

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 4

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 15

Enter the area of element number three in mm^2 10

Enter the area of element number four in mm^2 5

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the Young's modulus for element number four in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the length of element four in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number four in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number four 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 400

Enter the nodal force on node j for element number one 200

Enter the nodal force on node j for element number two 400

Enter the nodal force on node j for element number three 200

Enter the nodal force on node j for element number four 400

Press one for u1=0, two for u5=0, three for u1=u5=0 1

Results

Displacement at node two, three, four and five is 0.06 mm, 0.12333 mm, 0.17333 mm and 0.22333 mm

Reaction force at element one, node i is -2400 N

Strain in the element one, two, three and four is 0.0005, 0.00053333, 0.0004 and 0.0004

Stress in the element one, two, three and four is 100 MPa, 106.6667 MPa 80 MPa and 80 MPa

With second boundary condition (u₅=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 4

Enter the area of element number one in mm^2 5

Enter the area of element number two in mm^2 10

Enter the area of element number three in mm^2 15

Enter the area of element number four in mm^2 20

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the Young's modulus for element number four in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the length of element four in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number four in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number four 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one -400

Enter the nodal force on node j for element number one -200

Enter the nodal force on node j for element number two -400

Enter the nodal force on node j for element number three -200

Enter the nodal force on node j for element number four -400

Press one for u1=0, two for u5=0, three for u1=u5=0 2

Results

Displacement at node one, two, three and four is -0.22333 mm, -0.17333 mm, -0.12333 mm and -0.06 mm

Reaction force at element four, node j is 2400 N

Strain in the element one, two, three and four is 0.0004, 0.0004, 0.00053333 and 0.0005

Stress in the element one, two, three and four is 80 MPa, 80 MPa 106.6667 MPa and 100 MPa

With third boundary condition (u₅=u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 4

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 15

Enter the area of element number three in mm^2 10

Enter the area of element number four in mm^2 5

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the Young's modulus for element number four in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the length of element four in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number four in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number four 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Enter the nodal force on node j for element number three 0

Enter the nodal force on node j for element number four 0

Press one for u1=0, two for u5=0, three for u1=u5=0 3

Results

Displacement at node two, three and four is -0.0058 mm, -0.0002 mm, 0.0032 mm

Reaction force at element one, node i and element four, node j is -132 N and 132 N

Strain in the element one, two, three and four is -0.000158, -4.4e-05, -6.6e-05 and -0.000132

Stress in the element one, two, three and four is -31.6 MPa, -8.8 MPa -13.2 MPa and -26.4 MPa

Sample Calculations for problems involving variable cross-sectional area

Problem 5.1d, Fagan pg. 95

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 25

Enter the area of element number two in mm^2 15

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 10000

Press one for u1=0, two for u3=0, three for u1=u3=0 1

Results

Displacement at node two and three is 0.21 mm and 0.55333 mm

Reaction force at element one, node i is -10000 N

Strain in the element one and two is 0.002 and 0.0033333

Stress in the element one and two is 400 MPa and 666.6667 MPa