Computational Fluid Dynamics: 1 Core VS 2 Cores VS 3 Cores VS 4 Cores. (SolidWorks Flow Simulation Premium)

Problem:
Lift and Drag calculation for a rectangular wing (0.5 m span x 0.15 m chord x uniform cross section of NACA 9410), at horizontal velocity of 10 km/h.

A very simple CFD problem (with 177250 cells) solved on a 4 core Haswell-DT, varying the number of CPU cores used each time; CPU times are as below:

Cores used            CPU time taken for solution (1 travel convergence)
1/4 core                 782 seconds (2.31x slower)
2/4 cores               468 seconds (1.38x slower)
3/4 cores               380 seconds (1.12x slower)
4/4 cores               339 seconds (1x)

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