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