Showing posts with label analysis. Show all posts
Showing posts with label analysis. Show all posts

Friday 3 August 2018

Computational Fluid Dynamics (CFD) Analysis of a Horizontal Axis Tidal Turbine (Update 01)

     In this post, the results of numerical simulations carried out on a three-blade horizontal-axis tidal turbine blades are made available. These simulations are a part of larger project relating to a horizontal-axis tidal turbine.
 
     The turbine under investigation had a diameter of 4 m. The aero-foils employed include the NACA 4424 (chord length 0.25 m), NACA 4420 (chord length 0.2312 m), NACA 4418 (chord length 0.2126 m), NACA 4417 (chord length 0.1938 m), NACA 4416 (chord length 0.175 m), NACA 4415 (chord length 0.1562 m), NACA 4414 (chord length 0.1376 m), NACA 4413 (chord length 0.1188 m) and the NACA 4412 (chord length 0.1 m) at a distance of 0.4 m, 0.6 m, 0.8 m, 1 m, 1.2 m, 1.4 m, 1.6 m, 1.8 m, 2 m from the blade root, respectively [1].
 
     The numerical simulations were carried out on SolidWorks Flow Simulation Premium software. 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. To predict turbulent flows, the Favre-averaged Navier-Stokes equations are used. The software employs a Cartesian mesh which is created using the immersed boundary method.
 
     The simulations were carried out using the Local rotating region(s) (Sliding) technique. The mesh had a total of 589,763 cells. 132,045 cells were at the turbine blade-rotating region boundary. While the process of mesh generation, a mesh control was employed to refine the mesh near the turbine blades and at the boundary of the rotating region and the stationery computational domain. To refine the mesh for the mesh independent test, the curvature level was increased i.e. the mesh in the areas of interest was refined by a factor of 8. As a result, the refined mesh had a total of 916,577 cells. A total of 221,978 were present around the turbine blade. The computational mesh around the turbine blade is shown in Fig. 1. The refined mesh in the areas of high gradients is also clearly visible. The computational domain had dimensions of 2D x 2D x 2.4D, where D is the diameter of the turbine. The dimensions for the computational domain resemble to those in [2].

     Fig. 2 shows computational domain around the turbine blades. In Fig. 2, the curved teal arrow represents the direction of rotation of the turbine. The green arrows represent the respective co-ordinate directions. The brown arrow represents the direction of the force of gravity. The blue arrow represents the direction of the fluid velocity. The circular disc represents the rotating region.
 

Fig. 1, The computational mesh.
    
Fig. 2, The computational domain and the orientation of the boundary conditions.
 
     The simulations were carried out for a total of 5 tip-speed ratios for the turbine ranging from 2 to 10. The fluid; water, velocity was set at 2 m/s. The results are indeed, mesh independent. The mesh independence test was conducted on the design point of the turbine i.e. at the TSR of 6. The plot between the turbine tip-speed ratio; TSR, and the co-efficient of power is shown in Fig. 3. It can be clearly seen from Fig. 3 that the results are in close agreement with the results from [1, 3]. The CFD results from both studies are lower then the BEM, Blade-Element Momentum, results because the three-dimensional effects are not considered while implementing the BEM method.
 

Fig. 3, Turbine efficiency plot.
 
     Fig. 4 shows velocity streamlines colored by velocity magnitude, both of these features are drawn relative to the rotating reference frame, around the turbine blade cross-section at various tip-speed ratios. It can be seen from Fig. 4 that the turbine stalls at TSR of 2 due to a large positive angle of attack. It can also be seen that as the TSR increases, the angle of attack on the blade decreases, it is because of this reason that the power output from the turbine increases.
 
 
Fig. 4, Row 1, L-R; TSR of 2 and 4. Row 2, L-R; TSR of 6 and 8. Row 3, TSR of 10.
 
     Thank you for reading. If you would like to collaborate with research projects, please reach out.

     [1] Binoe E. Abuan, and Robert J. Howell, "The Influence of Unsteady Flow to the Performance of a Horizontal Axis Tidal Turbine," Proceedings of the World Congress on Engineering. London, 2018.
     [2] Ibrahim, I. H., and T. H. New, "A numerical study on the effects of leading-edge modifications upon propeller flow characteristics," Proceedings of Ninth International Symposium on Turbulence and Shear Flow Phenomena. Melbourne, 2015.
     [3] Bahaj, A., Batten, W., & McCann, J. (2007). Experimental verifications of numerical predictions for the hydrodynamic performance of horizontal axis marine current turbines. Renewable Energy, 2479-2490.

Update 01:

     CFD post processing added. One more TSR simulated.

Thursday 16 June 2016

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 (u1=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 (u2=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 (u1=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 (u3=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 (u3=u1=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 (u1=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 (u4=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 (u4=u1=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 (u1=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 (u5=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 (u5=u1=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