1D Super-sonic nozzle flow simulation using Macormack Method

IN THIS PROJECT, I HAVE WORKED ON THE ONE DIMENSIONAL SUPERSONIC NOZZLE FLOW. IN THIS FLOW, I HAVE ASSUMED THE FLOW TO BE STEADY,ISENTROPIC. IN THIS NOZZLE, I HAVE CONSIDERED THE FLOW AT INLET SECTION COMES FROM RESERVOIR WHERE THE PRESSURE, TEMPERATURE ARE DENOTED AS P_0 AND T_0 RESPECTIVELY. THE CROSS-SECTIONAL AREA OF RESERVOIR IS INFINITE DUE TO WHICH THE VELOCITY BECOMES ZERO THERE.THE NOZZLE IS BASICALLY DIVIDED INTO TWO SECTIONS i.e CONVERGENT SECTION AND DIVERGENT SECTION. AS A RESULT, THE FLOW EXPANDS ISENTROPICALLY AND ATTAINS SUPERSONIC(M>1) SPEED AT EXIT OR DIVERGENT SECTION OF THE NOZZLE WHERE IT BECOMES SUBSONIC(M<1) AT THE CONVERGENT SECTION AND AT A POINT IT'S SPEED  BECOMES EQUAL TO THE SONIC SPEED(M=1) WHICH IS KNOWN AS THROTTLE POINT(CENTRE POINT OF CONVERGENT AND DIVERGENT SECTION).

IT IS ALSO KNOWN AS THE QUASI-ONE DIMENSIONAL NOZZLE FLOW BECAUSE IN THE  FLOW, THE FLOW PROPERTIES ARE CONSIDERED AS THE FUNCTION OF DIRECTION OF THE FLOW.

THE GOVENING EQUATIONS USED IN THIS FLOW ARE:-

CONTINUITY EQUATION:-

rho_1*A_1+ rho_1*A_1*(v_1)^2+int_(A_1)^(A_2)rho*dA=rho_2*A_2+rho_2*A_2*(v_2)^2

MOMENTUM EQUATION :-

h_1+(v_1)^2/2=h_2+(v_2)^2/2

THE TECHNIQUE USED IN THIS NOZZLE FLOW FOR NUMERICAL SIMULATION SO AS TO FIND NATURE OF FLOW IS :-

MacCormack method:-

IN COMPUTATIONAL FLUID DYNAMICS, IT IS WIDELY USED IN DISCRETIZATION SCHEME FOR NUMERICAL SIMULATION OF HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS.IT IS A SECOND ORDER DIFFERENCE METHOD.

The MacCormack method is a variation of the two-step Lax–Wendroff scheme but is much simpler in application. To illustrate the algorithm, consider the following first order hyperbolic equation

The application of MacCormack method to the above equation proceeds in two steps; a predictor step which is followed by a corrector step.

Predictor step: In the predictor step, a "provisional" value of  is estimated as follows

The above equation is obtained by replacing the spatial and temporal derivatives in the previous first order hyperbolic equation using forward differences.

Corrector step: In the corrector step, the predicted value  is corrected according to the equation

Note that the corrector step uses backward finite difference approximations for spatial derivative. Note also that the time-step used in the corrector step is in contrast to the  used in the predictor step.

Replacing the  term by the temporal average

to obtain the corrector step as

.

IN THIS PROJECT, I HAVE WORKED ON THE TWO FORMS OF NOZZLE FLOW EQUATION:-

1.) CONSERVATIVE FORM.

2.) NON-CONSERVATIVE FORM.

1.) NON-CONSERVATIVE FORM:- IN THIS FORM, WE HAVE CONSIDERED THE FLOW TO BE CONSEREVED i.e THE FLOW DOMAIN IS CONSTANT.

THE GOVERNING EQUATION IN THIS FORM IS:-

HENCE, THE GOVERNING EQUATION IN THE FORM OF DIFFERENCE EQUATTION IS:-

2.) CONSERVATION FORM:- IN THIS FORM, WE HAVE CONSIDERED THE FLOW TO BE UNCONSEREVED i.e THE FLOW DOMAIN IS MOVING.

THE GOVERNING EQUATION IN THIS FORM IS:-

, CONTINUITY

,MOMENTUM

,   ENERGY

NOW, IN ORDER TO RUN THE GOVERNING EQUATION, WE HAVE TO PROVIDE BOUNDARY CONDITION AND INTIAL CONDITION SO AS TO MAKE THE SIMULATION PROCESS MORE EASIER.

INTIAL CONDITION:-

1.) NON-CONSERVATIVE:-

FOR DENSITY,     rho=1-0.3146x

FOR TEMPERATURE,  T=1-0.2314X

FOR VELOCITY,   V=(0.1+1.09x)T^(1/2)

2.) CONSERVATIVE :-

NOZZLE SHAPE AND ITS CONDITION FOR BOTH :-

A=1+2.2(x-1.5)^2 , 0lexle3

BOUNDARY CONDITION:-

SUBSONIC INFLOW BOUNDARY CONDITION:-

V_1=2V_2-V_3

rho=1

T=1

SUBSONIC INFLOW BOUNDARY CONDITION:-

V_n=2*V_(n-1)-V_(n-2)

rho_n=2rho_(n-1)-rho_(n-2)

T_n=2T_(n-1)-T_(n-2)

CALCULATION OF TIME STEP:-

HERE TIME STEP ALSO PLAYS AN IMPORTANT ROLE IN BOTH CASE AS IT IMPACT IS HIGH ON THE CONVERGENCE OF VALUE . SO IN ORDER TO MAKE SIMULATION MORE CONVINIENT WE HAVE TO PROVIDE SUCH TIME STEP VALUE HERE WHICH WILL MAKE SIMULATION RUN EASIER AND WOULD MAKE US EASIER TO ATTAIN VALUE.

CRITERIA OF TIME STEP:-

Deltat=C*((Deltax)/(a +v))

where, 'a' is known as as the speed of sound and speed of fluid is known as 'v'.

'C' IS THE 'COURANT-NUMBER'. WHERE IT SHOULD BE LESS THAN 1 FOR CONVERGENCE OF VALUE.

AFTER CALCULATING THE VALUE OF TIME STEP AT EACH GRID POINT OF NOZZLE, WE FIND OUT MINIMUM VALUE OF TIME STEP FROM IT WHICH WILL DO CONVERGENCE EASILY AT ALL GRID POINTS.

i.e min(Deltat_1+Deltat_2+Deltat_3+....Deltat_n)

MAIN CODE:-

clear all
close all
clc

% Time steps:-
nt=[100 400 700 1400];
n=1;

while(n|=0)

n=input("Pick a number, any number for grid points for eg.31,61!\n say 0 for no output");
x=linspace(0,3,n);
dx=x(2)-x(1);

for j=1:4
[sim_tim_conv]=nozzle_flow_non_conservative(n,x,dx,nt(j));
end

for j=1:4
[sim_tim_non_conv]=nozzle_flow_conservative(n,x,dx,nt(j));
end

end

##### PLEASE INPUT VALUE IN COMMAND WINDOW  ####

FUNCTION DEFINED FOR CONSERVATIVE FORM:-

function [sim_tim]=nozzle_flow_conservative(n,x,dx,nt)

% calculation of Intial profile:-
a=1+(2.2*(x-1.5).^2);

for i=1:n
if (x(i)>=0&&x(i)<=0.5)
rho(i)=1;
T(i)=1;
end

if (x(i)>=0.5&&x(i)<=1.5)
rho(i)=1-(0.366*(x(i)-0.5));
T(i)=1.0-(0.167*(x(i)-0.5));
end

if (x(i)>=1.5&&x(i)<=3.5)
rho(i)=0.634-(0.3879*(x(i)-1.5));
T(i)=0.833-(0.3507*(x(i)-1.5));
end
end

% Intial values:-
v=0.59./(rho.*a);
gamaa=1.4;

% continuity equation:-

U_1=rho.*a;

% momentum equation:-
U_2=rho.*a.*v;

% energy equation:-
U_3=rho.*((T./(gamaa-1))+((gamaa/2)*(v.^2))).*a;

% Time steps:-

a_conv=sqrt(T);
c=0.5;
dt=(c*dx)./(a_conv+v);
dt=min(dt);
iter=0;

tic()
% Time loop:-
for t=1:nt

U_1_old=U_1;
U_2_old=U_2;
U_3_old=U_3;
iter=iter+1;

% predictor method:-

% flux terms:-
F_1=U_2;

F_2=((U_2.^2)./U_1)+((gamaa-1)/gamaa)*(U_3-((gamaa/2)*((U_2.^2)./U_1)));

F_3=gamaa*((U_2.*U_3)./U_1)-((gamaa*(gamaa-1))/2)*(U_2.^3./U_1.^2);

for i=2:n-1

% source term:-
J_2(i)=(1/gamaa)*rho(i)*T(i)*(a(i+1)-a(i))/dx;

% finding predicted values of solution terms using forward difference:-
U_1_dt_p(i)= -(F_1(i+1)-F_1(i))/dx;

U_2_dt_p(i)= -(F_2(i+1)-F_2(i))/dx+J_2(i);

U_3_dt_p(i)= -(F_3(i+1)-F_3(i))/dx;

% solution update:-
U_1(i)=U_1(i)+ U_1_dt_p(i)*dt;

U_2(i)=U_2(i)+ U_2_dt_p(i)*dt;

U_3(i)=U_3(i)+ U_3_dt_p(i)*dt;

end

% updated predicted values of density and temperature:-
rho=U_1./a;
T=(gamaa-1)*(U_3./U_1-(gamaa/2)*((U_2.^2/U_1)));
v=U_2./U_1;

% updated values of flux terms:-

F_1=U_2;

F_2=((U_2.^2)./U_1)+((gamaa-1)/gamaa)*(U_3-(gamaa/2)*(U_2.^2./U_1));

F_3=gamaa*((U_2.*U_3)./U_1)-((gamaa*(gamaa-1))/2)*(U_2.^3./U_1.^2);

% corrector step:-

for i=2:n-1

% updated value of source term:-

J_2(i)=(1/gamaa)*rho(i)*T(i)*(a(i)-a(i-1))/dx;

% finding solution terms using reaward differennce:-
U_1_dt_c(i)= -(F_1(i)-F_1(i-1))/dx;

U_2_dt_c(i)= -(F_2(i)-F_2(i-1))/dx +J_2(i);

U_3_dt_c(i)= -(F_3(i)-F_3(i-1))/dx;

end

% average time derivative is as folow:-

U_1_dt_av= 0.5*(U_1_dt_p+U_1_dt_c);

U_2_dt_av= 0.5*(U_2_dt_p+U_2_dt_c);

U_3_dt_av= 0.5*(U_3_dt_p+U_3_dt_c);

% final solution update:-
for j=2:n-1

U_1(j)= U_1_old(j)+U_1_dt_av(j)*dt;

U_2(j)= U_2_old(j)+U_2_dt_av(j)*dt;

U_3(j)= U_3_old(j)+U_3_dt_av(j)*dt;

end

## applying boundary conditions:-
% inflow boundary condition:-
U_2(1)=2*U_2(2)-U_2(3);

% outflow boundary condition:-
U_1(n)=2*U_1(n-1)-U_1(n-2);
U_2(n)=2*U_2(n-1)-U_2(n-2);
U_3(n)=2*U_3(n-1)-U_3(n-2);

end

sim_tim=toc();

% final value of density :-
rho=U_1./a;

% final value of temperature:-
T=(gamaa-1)*(U_3./U_1-(gamaa/2)*((U_2./U_1).^2));

% final value of velocity:-
v=U_2./U_1;

% final value of mass flow:-
mass_flow_rate=rho.*a.*v;

% final value of MACH number:-
M=v./a_conv;

% final value of pressure:-
P=rho.*T;

text=sprintf('CONSERVATIVE FLOW \n NUMBER OF GRID POINTS USED IN NOZZLE LENGTH=%d,\n NUMBER OF TIME STEPS=%d',n,t);

if(t==100)
figure(6)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==400)
figure(7)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==700)
figure(8)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==1400)
figure(9)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

FUNCTION DEFINED FOR NON-CONSERVATIVE FORM:-

function [sim_tim]=nozzle_flow_non_conservative(n,x,dx,nt)

% calculation of Intial profile:-
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^(1/2);
a=1+2.2*(x-1.5).^2;

gamaa=1.4;

% calculation of minimum time steps:-
a_conv=sqrt(T);
c=1/2;
dt=c*dx./(a_conv+v);
dt=min(dt);

% intial iteration counter:-
iter=0;

tic()
% Time loop:-
for t=1:nt
rho_old=rho;
v_old=v;
T_old=T;
iter=iter+1;

% predictor method:-

for i=2:n-1

dv_dx=(v(i+1)-v(i))/dx;
drho_dx=(rho(i+1)-rho(i))/dx;
da_dx=(log(a(i+1))-log(a(i)))/dx;
dT_dx=(T(i+1)-T(i))/dx;

% continuity equation
drho_dt_p(i)=-rho(i)*dv_dx-rho(i)*v(i)*da_dx-v(i)*drho_dx;

% momentum equation
dv_dt_p(i)=-v(i)*dv_dx-(1/gamaa)*(dT_dx+(T(i)/rho(i))*drho_dx);

% energy equation
dT_dt_p(i)=-v(i)*dT_dx-(gamaa-1)*T(i)*(dv_dx+v(i)*da_dx);

% Solution update
v(i)=v(i)+dv_dt_p(i)*dt;
rho(i)=rho(i)+drho_dt_p(i)*dt;
T(i)=T(i)+dT_dt_p(i)*dt;

end

% corrector method:-

for i=2:n-1
dv_dx=(v(i)-v(i-1))/dx;
drho_dx=(rho(i)-rho(i-1))/dx;
da_dx=(log(a(i))-log(a(i-1)))/dx;
dT_dx=(T(i)-T(i-1))/dx;

% continuity equation
drho_dt_c(i)=-rho(i)*dv_dx-rho(i)*v(i)*da_dx-v(i)*drho_dx;

% momentum equation
dv_dt_c(i)=-v(i)*dv_dx-1/gamaa*(dT_dx+(T(i)/rho(i))*drho_dx);

% energy equation
dT_dt_c(i)=-v(i)*dT_dx-(gamaa-1)*T(i)*(dv_dx+v(i)*da_dx);

end

% compute the average time derivative:-
drho_dt_av=0.5*(drho_dt_p+drho_dt_c);
dv_dt_av=0.5*(dv_dt_p+dv_dt_c);
dT_dt_av=0.5*(dT_dt_p+dT_dt_c);

% Final solution update:-

for j=2:n-1
rho(j)=rho_old(j)+drho_dt_av(j)*dt;
v(j)=v_old(j)+dv_dt_av(j)*dt;
T(j)=T_old(j)+dT_dt_av(j)*dt;
end

% apply boundary condition:-
% Inlet:-
v(1)=2*v(2)-v(3);

% Outlet:-
v(n)=2*v(n-1)-v(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
T(n)=2*T(n-1)-T(n-2);

% MACH value:-
M=v./a_conv;

% pressure force value:-
P=rho.*T;

% Throttle values for parameter:-
rho_throttle_value(t)=rho(((n-1)/2)+1);
T_throttle_value(t)=T(((n-1)/2)+1);
P_throttle_value(t)=P(((n-1)/2)+1);
M_throttle_value(t)=M(((n-1)/2)+1);

end
sim_tim=toc();

% mass flow rate:-
mass_flow_rate=rho.*v.*a;

text=sprintf('NON-CONSERVATIVE FLOW \n NUMBER OF GRID POINTS USED IN NOZZLE LENGTH=%d,\n NUMBER OF TIME STEPS=%d',n,t);

if(t==100)
figure(1)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==400)
figure(2)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==700)
figure(3)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==1400)
figure(4)
plot(x,mass_flow_rate,'b')
axis([0 3 0.4 0.7])
xlabel('length of nozzle (x/L)')
ylabel('mass flow rate(\rho*a*v)/(\rho_0*a_0*A)')
title(text)
end

if(t==1400)
figure(5)
subplot(4,1,1)
plot(rho_throttle_value,'b')
xlabel('Number of iteration')
ylabel('density (\rho)/(\rho_0)')
title(text)

subplot(4,1,2)
plot(T_throttle_value,'r')
xlabel('Number of iteration')
ylabel('temperature(T/T_0)')
title(text)

subplot(4,1,3)
plot(P_throttle_value,'g')
xlabel('Number of iteration')
ylabel('pressure (P/P_0)')
title(text)

subplot(4,1,4)
plot(M_throttle_value,'y')
xlabel('Number of iteration')
ylabel('mach number(M)')
title(text)

end



OUTPUT :-

1.) NON-CONSERVATIVE :-

2.) CONSERVATIVE :-

HERE I OBSERVED THAT THE VALUE OF MASS FLOW IS CONVERGING IN NON-CONSERVATIVE FORM BETTER THAN THE CONSERVATIVE FORM.

IN BOTH FORM, THE VALUE WAS TOO LARGE FOR LESS TIME STEP i.e UPTO 400 WHICH CAN BE EASILY SEEN IN THE PLOTTED GRAPH.

AT 700 TIME STEP, CONVERGENCE WAS MORE IN CASE OF NON-CONSERVATIVE STATE THAN THE CONSERVATIVE STATE.

AT 1400 TIME STEP, CONVERGENCE WAS HIGHLY ATTAINED THERE IN CONSERVATIVE WHILE IT WAS  STILL VARYING VERY SMALL MARGINALLY IN CASE OF NON-CONSERVATIVE FORM.

SO IT IS CONCLUDED THAT NON-CONSERVATIVE HAS ADVANTAGE OF ATTAINING CONVERGENCE IN LESS TIME STEP THAN THE CONSERVATIVE FORM.

AND NON-CONSERVATIVE TERM HAVE ONE MORE ADVANTAGE THAT IT IS MORE EXACT NEARER TO THE EXACT VALUE OF MASS FLOW.

WHEN BOTH CASE IS TESTED BY THE EFFECT OF VARYING GRID POINTS, IT IS NOTICED THAT THERE IS HUGE DIFFERENCE DIFFERENCE BETWEEN BOTH CASE WHEN COMPARED AT MINIMUM TIME STEP. WHERE AS, IN CASE OF HIGH TIME STEP OR AT AND AFTER THE TIME STEP VALUE AT WHICH CONVERGENCE OCCUR, IT HAVE VERY LESS DIFFERNCE i.e AROUND ORDER OF 10^-2.

SO IT IS CONCLUDED THAT THERE IS NO EFFECT OF GRID SIZE VARIATION ONCE THE CONVERGENCE IT IS ACHIEVED.

HERE WHEN BOTH IS COMPARED IN TERMS OF SIMULATION TIME, IT IS NOTICED THAT NON-CONSERVATIVE CONSUMES LESS TIME FOR SIMULATION WHEN COMPARED TO THE CONSERVATIVE FORM.

AT 31 GRID POINTS:-

SIMULATION ON PRANDTL MEYER SHOCK PROBLEM AND UNDERSTANDING THE CONCEPT BEHIND IT Shyam Rai · 2020-01-08 20:30:44

THEORY:- BOUNDARY CONDITION :- IN PHYSICS ,THESE ARE THOSE CONDITIONS WHICH ARE APPLIED AT THE BOUNDARY OF THE DOMAIN . IT SPECIFIES THE CONDITION AT BOUNDARY ON THE BASIS OF WHAT THE GOVERNING EQUATION DEFINES THE CONDITION OF DOMAIN ALONG WITH THE HELP OF INTIAL COND Read more

Transient simulation of flow over a throttle body Shyam Rai · 2019-11-09 08:36:53

OBJECTIVE:- STEADY STATE SIMULATION OF FLOW OVER A THROTTLE BODY . SOFTWARE USED :- 1. CONVERGE STUDIO:- TO SETUP THE MODEL.  2. CYGWIN:- FOR RUNNING THE SIMULATION. 3. OPENFOAM:- TO VISUALIZE THE DIFFERENT POST-PROCESSED RESULT.   THEORY:- HERE, WE HAV Read more

Steady state simulation of flow over a throttle body Shyam Rai · 2019-11-08 21:36:56

OBJECTIVE:- STEADY STATE SIMULATION OF FLOW OVER A THROTTLE BODY . SOFTWARE USED :- 1. CONVERGE STUDIO:- TO SETUP THE MODEL.  2. CYGWIN:- FOR RUNNING THE SIMULATION. 3. OPENFOAM:- TO VISUALIZE THE DIFFERENT POST-PROCESSED RESULT.   THEORY:- HERE, WE HAV Read more

COMPARISON OF WEDGE BOUNDARY CONDITION WITH SYMMETRIC BOUNDARY CONDITION AND HAIGEN-POISEUILLE EQUATION Shyam Rai · 2019-09-09 19:47:00

AIM:- IN THIS PROJECT SIMULATION IS DONE ON THE FLUID FLOW THROUGH PIPE USING OPENFOAM SOFTWARE AND MATLAB . THE TOPIC GIVEN BELOW IS COVERED IN THIS                      PROJECT. TO MAKE A PROGRAM IN MATLAB THAT Read more

Simulation of Flow through pipe in OpenFoam using Matlab software Shyam Rai · 2019-09-09 12:53:45

AIM:- IN THIS PROJECT SIMULATION IS DONE ON THE FLUID FLOW THROUGH PIPE USING OPENFOAM SOFTWARE AND MATLAB . THE TOPIC GIVEN BELOW IS COVERED IN THIS PROJECT. TO MAKE A PROGRAM IN MATLAB THAT CAN GENERATE MESH AUTOMATICALLY FOR ANY WEDGE ANGLE AND GRADING SCHEME. TO Read more

SIMULATION OF FLOW THROUGH THE BACKWARD FACING STEP BY USING BLOCKMESH IN OPENFOAM Shyam Rai · 2019-08-13 20:27:58

AIM :- IN THIS PROJECT, I HAVE WORKED ON SIMULATING THE FLOW OF FLUID THROUGH  THE DOMAIN (CONTROL VOLUME). DESCRIPTION:- THE DOMAIN IS NOTHING ELSE, IT IS JUST A CONTROL VOLUME THROUGH WHICH I HAVE SIMULATED THE RESULT OBTAINED FROM THE GOVERNING EQUATION. THE Read more

PROGRAM ON FLOW OF HEAT IN THE RECTANGULAR SLAB (2-D DOMAIN) USING 2-DIMENSIONAL HEAT CONDUCTION EQUATION IN STEADY AND TRANSITION STATE Shyam Rai · 2019-07-09 18:02:52

In this project, I have worked on the simulation of conduction of temperature heat in two dimensional space (2D domain).  The aim was to determine the rate of flow of  heat due to temperature difference over the space domain under steady condition and Read more

Program on deriving fourth order approximation for different difference schemes for second order derivative Shyam Rai · 2019-06-23 15:26:37

In this programming, I have derived the four order approximation of second order derivative using finite difference method. In Finite difference method,we generally convert ordinary differential equation into difference equation. In order to achieve so, we use simple t Read more

SIMULATION OF AN OSCILLATION OF A 2-D PENDULUM BY USING SECOND ORDER ORDINARY DIFFERENTIAL EQUATION Shyam Rai · 2019-06-21 19:15:09

In this report, I have simulated the oscillation of 2D pendulum and also generate a plot of \'angular_displacement\' & \'angular_velocity\' vs \"time\" in octave. Consider a pendulum which is having a string of length of \'1metre\' connected to a ball of mass,m=1kg Read more