## 7 - 1D Supersonic Nozzle Flow Simulation - Macormack Method

This project main objective is the numerical solution of a 1D Supersonic Nozzle Flow using the Macormack Method.

Although analytically identical, the numerical solution of conservation and nonconservation forms of the governing equations is usually different. That is the reason why both ways are solved and compared in this project. For that, the number of iterations and computation time will be compared, as long as the final profile for the normalized mass flow rate. During this project,both a grid dependence test and a CFL-based adaptative time loop control was programmed to assure the fidelity of the numerical results.

1) EQUATIONS

The conservation and nonconservation equations are analytical identical. Both come from applying the integral conservation method to the control volume represented by the nozzle, and therefore the analytical solution is identical. However, in numerical analysis, the two sets of equations cant be used indistinguisible, and a deeper analysis must be carried out depending on the problem.

As a general rule when the control volume is fixed in space with the fluid moving through the governing equations will be in conservation form. However, when the control volume is moving with the fluid in a sense that same fluid particles are always remain inside the control volume, the calculation will generally performed in nonconservation form.

The derivation of these sets of equations for the current problem is not in the scope of this project. This derivation is pretty simple and can be found in any fluid mechanics book. In needed, the full derivation and the equations themselves can be found at the Chapter 7 of the book: John D. Anderson Jr. Computational Fluid Dynamics: The Basics with Applications.

2) SCRIPTS

Main Program

clear all
close all
clc

%Macormack Method for solving a 1D Supersonic Nozzle Flow
%Non-Conservative and Conservative Approach

%Inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;

%Data for the CFL Time Loop and Convergence Controler
C=0.5; %Courant Number
error=1e9;
error2=1e9;
tol=1e-5;

%Area Profile
A=1+2.2*(x-1.5).^2;
n_throat=find(A-1<0.001); %The value where A/A* is one (Or close)

%---------------NON-CONSERVATIVE FORM SOLVER-----------------------------
tic

[rho,v,T,p,M,massflow,rho_throat,T_throat,p_throat,M_throat,k]=NonConservative_1DSupersonicNozzleMacormack(n,x,dx,gamma,A,n_throat,C,error,tol);
time_nc = toc;
fprintf('Minimum Number of Iterations for Convergence (NonConservative Form): %f\n',k);
fprintf('Computation time (NonConservative Form): %f s\n',time_nc);

%TRANSIENT PROPERTIES AT THE THROAT

%Transient variation of the throat properties
Transient_Throat=figure(2);
subplot(4,1,1);
plot(rho_throat,'m');
title('Timewise variations of the properties at the nozzle throat (Non-Conservative Form)');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
axis([0 k 0.5 1]);
grid minor;

subplot(4,1,2);
plot(T_throat,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 k 0.6 1.2]);
grid minor;

subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 k 0.4 1]);
grid minor;

subplot(4,1,4);
plot(M_throat,'b');
legend('Mach Number');
ylabel('M');
axis([0 k 0.8 1.3]);
grid minor;
xlabel('Number of iterations');

%Steady profile of the properties at the nozzle
figure(3)
subplot(4,1,1);
plot(x,rho,'m');
title('Properties profiles in the 1D Supersonic Nozzle flow (Non-Conservative Form)');
ylabel('\rho/\rho_{o}');
legend('Density ratio','location','southwest');
grid minor;

subplot(4,1,2);
plot(x,T,'r');
ylabel('T/T_{o}');
legend('Temperature ratio','location','southwest');
grid minor;

subplot(4,1,3);
plot(x,p,'g');
ylabel('p/p_{o}');
legend('Pressure ratio','location','southwest');
grid minor;

subplot(4,1,4);
plot(x,M,'b');
ylabel('M');
legend('Mach Number','location','northwest');
grid minor;

%---------------CONSERVATIVE FORM SOLVER-----------------------------------

tic

[rho_con,v_con,T_con,p_con,M_con,massflow_con,rho_throat_con,T_throat_con,p_throat_con,M_throat_con,k2]=Conservative_1DSupersonicNozzleMacormack(n,x,dx,gamma,A,n_throat,C,error,tol);
time_c = toc;
fprintf('Minimum Number of Iterations for Convergence (Conservative Form): %f\n',k2);
fprintf('Computation time (Conservative Form): %f s\n',time_c);

%TRANSIENT PROPERTIES AT THE THROAT

%Transient variation of the throat properties
Transient_Throat_con=figure(5);
subplot(4,1,1);
plot(rho_throat_con,'m');
title('Timewise variations of the properties at the nozzle throat (Conservative Form)');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
axis([0 k2 0.5 1]);
grid minor;

subplot(4,1,2);
plot(T_throat_con,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 k2 0.6 1.2]);
grid minor;

subplot(4,1,3);
plot(p_throat_con,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 k2 0.4 1]);
grid minor;

subplot(4,1,4);
plot(M_throat_con,'b');
legend('Mach Number');
ylabel('M');
axis([0 k2 0.8 1.3]);
grid minor;
xlabel('Number of iterations');

%Steady profile of the properties at the nozzle
figure(6)
subplot(4,1,1);
plot(x,rho_con,'m');
title('Properties profiles in the 1D Supersonic Nozzle flow (Conservative Form)');
ylabel('\rho/\rho_{o}');
legend('Density ratio','location','southwest');
grid minor;

subplot(4,1,2);
plot(x,T_con,'r');
ylabel('T/T_{o}');
legend('Temperature ratio','location','southwest');
grid minor;

subplot(4,1,3);
plot(x,p_con,'g');
ylabel('p/p_{o}');
legend('Pressure ratio','location','southwest');
grid minor;

subplot(4,1,4);
plot(x,M_con,'b');
ylabel('M');
legend('Mach Number','location','northwest');
grid minor;

%Steady massflow rate profile (Comparisson Non-Conservative vs Conservative)
figure(7)
plot(x,massflow)
hold on
plot(x,massflow_con)
legend('Non-Conservative','Conservative')
grid minor

Nonconservation Solver

function [rho,v,T,p,M,massflow,rho_throat,T_throat,p_throat,M_throat,k]=NonConservative_1DSupersonicNozzleMacormack(n,x,dx,gamma,A,n_throat,C,error,tol)

%Macormack Method for solving a 1D Supersonic Nozzle Flow
%Non-Conservative Approach

%Calculate Initial Profiles
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;

%Initialization for the time loop
k=0;

%Outer time loop

while error>tol

k=k+1;

rho_old=rho;
v_old=v;
T_old=T;
a_old=sqrt(T_old);
massflow_old=rho_old.*v_old.*A;

%CFL Time Loop Control
TimeStep=C*dx./(a_old+v_old);
dt=min(TimeStep);

%Predictor Step
for j=2:n-1

dvdx_f=(v(j+1)-v(j))/dx;
drhodx_f=(rho(j+1)-rho(j))/dx;
dTdx_f=(T(j+1)-T(j))/dx;

%Continuity Equation

%Momentum Equation
dvdt_p(j)=-v(j)*dvdx_f -(1/gamma)*(dTdx_f + T(j)/rho(j)*drhodx_f);

%Energy Equation

%Solution Update
rho(j)=rho(j)+drhodt_p(j)*dt;
v(j)=v(j)+dvdt_p(j)*dt;
T(j)=T(j)+dTdt_p(j)*dt;
end

%Corrector Step
for j=2:n-1
%drho_dt_p= -rho(j)*(v(j+1)-v(j))/dx -rho(j)*v(j)*(log(A(j+1))-log(A(j)))/dx - v(j)*(rho(j+1)-rho(j))/dx
dvdx_b=(v(j)-v(j-1))/dx;
drhodx_b=(rho(j)-rho(j-1))/dx;
dTdx_b=(T(j)-T(j-1))/dx;

%Continuity Equation

%Momentum Equation
dvdt_c(j)=-v(j)*dvdx_b -(1/gamma)*(dTdx_b + T(j)/rho(j)*drhodx_b);

%Energy Equation

end

%Compute the Average Time derivative
drhodt_avg=0.5*(drhodt_p+drhodt_c);
dvdt_avg=0.5*(dvdt_p+dvdt_c);
dTdt_avg=0.5*(dTdt_p+dTdt_c);

for i=2:n-1

%Final Solution Update
rho(i)=rho_old(i)+drhodt_avg(i)*dt;
v(i)=v_old(i)+dvdt_avg(i)*dt;
T(i)=T_old(i)+dTdt_avg(i)*dt;

end

%Apply the boundary conditions

%Inlet
rho(1)=1;
v(1)=2*v(2)-v(3);
T(1)=1;

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

%Calculation of other variables of interest
M=v./sqrt(T); %Mach Number
massflow=rho.*v.*A; %Non-dimensional Mass Flow Rate
p=rho.*T; %Non-dimensional pressure

%Transient Throat Values
rho_throat(k)=rho(n_throat);
T_throat(k)=T(n_throat);
p_throat(k)=p(n_throat);
M_throat(k)=M(n_throat);

%Check for time convergence
error=max(abs(massflow-massflow_old));

%-----TRANSIENT MASS FLOW RATE PROFILE----------------
%We want to know how the massflow rate changes with time before the steady
%solution
TransientMassFlow=figure(1);
hold on;

if k== 1
plot(x,massflow,'--k');
else if k == 50
plot(x,massflow,'r');
else if k == 100
plot(x,massflow,'g');
else if k == 150
plot(x,massflow,'c');
else if k == 200
plot(x,massflow,'b');
else if k == 700
plot(x,massflow,'-.b','linewidth',1.25);
plot(x(end),massflow(end),'k*');
end
end
end
end
end
end

end

title('Transient Non-dimensional Mass Flow Rates(Non-Conservative Form)');
xlabel('x/L');
ylabel('\rhoVA/\rho_{o}V_{o}A_{o}');
legend('0\Deltat','50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat');
axis([0 3 0.2 1.8]);
grid minor;

end


Conservation Form Solver

function [rho_con,v_con,T_con,p_con,M_con,massflow_con,rho_throat_con,T_throat_con,p_throat_con,M_throat_con,k2]=Conservative_1DSupersonicNozzleMacormack(n,x,dx,gamma,A,n_throat,C,error2,tol)

%Macormack Method for solving a 1D Supersonic Nozzle Flow Conservative Approach

%Calculate Initial Profiles
for i = 1:length(x)
if x(i)>=0 && x(i)<0.5
rho_con(i)=1;
T_con(i)=1;
else if x(i)>=0.5 && x(i)< 1.5
rho_con(i)=1.0-0.366*(x(i)-0.5);
T_con(i)=1.0-0.167*(x(i)-0.5);
else if x(i)>=1.5 && x(i)<=3.0
rho_con(i)=0.634-0.3879*(x(i)-1.5);
T_con(i)=0.833-0.3507*(x(i)-1.5);
end
end
end
end
v_con =0.59./(rho_con.*A);
p_con=rho_con.*T_con;

%Calculate the initial solution vectors
U1=rho_con.*A;
U2=rho_con.*v_con.*A;
U3=rho_con.*A.*((T_con./(gamma - 1)) + (gamma/2)*v_con.^2);

%Initialization for the time loop
k2=0;

%Outer time loop with a convergence criteria
while error2>tol

%Convergence variables
k2=k2+1;
massflow_old_con=rho_con.*v_con.*A;

%We store the Old Solution Vectors
U1_old=U1;
U2_old=U2;
U3_old=U3;

%Defining the flux vectors:
F1= U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);

a_con=sqrt(T_con);
TimeStep=C*dx./(a_con+v_con);
dt=min(TimeStep);

%Predictor Step
for j=2:n-1

% Defining the source term:
J2(j)=(1/gamma).*rho_con(j).*T_con(j).*((A(j+1)-A(j))/dx);

%Continuity Equation
dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;

%Momentum Equation
dU2dt_p(j)=-(F2(j+1)-F2(j))/dx+J2(j);

%Energy Equation
dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;

%Solution Update
U1(j)=U1(j)+dU1dt_p(j)*dt;
U2(j)=U2(j)+dU2dt_p(j)*dt;
U3(j)=U3(j)+dU3dt_p(j)*dt;
end

% Predicted values of the primitive variables and flux vectors

rho_con=U1./A;
v_con=U2./U1;
T_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
p_con=rho_con.*T_con;

F1=U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);

%Corrector Step
for j=2:n-1

%Updating the source term:
J2(j)=(1/gamma)*rho_con(j).*T_con(j)*((A(j)-A(j-1))/dx);

%Continuity Equation
dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;

%Momentum Equation
dU2dt_c(j)=-(F2(j)-F2(j-1))/dx+J2(j);

%Energy Equation
dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
end

%Compute the Average Time derivative
dU1dt_avg=0.5*(dU1dt_p+dU1dt_c);
dU2dt_avg=0.5*(dU2dt_p+dU2dt_c);
dU3dt_avg=0.5*(dU3dt_p+dU3dt_c);

%Final Solution Update
for i=2:n-1

U1(i)=U1_old(i)+dU1dt_avg(i)*dt;
U2(i)=U2_old(i)+dU2dt_avg(i)*dt;
U3(i)=U3_old(i)+dU3dt_avg(i)*dt;

end

%Apply the boundary conditions

%Inlet
rho_con(1)=1;
T_con(1)=1;

U1(1)=rho_con(1)*A(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*(T_con(1)/(gamma-1)+(gamma/2)*v_con(1).^2);

%Outlet
U1(n)=2*U1(n-1)-U1(n-2);
U2(n)=2*U2(n-1)-U2(n-2);
U3(n)=2*U3(n-1)-U3(n-2);

%Calculation of the new primitive variables
rho_con=U1./A;
v_con=U2./U1;
T_con=(gamma -1)*(U3./U1-(gamma/2)*v_con.^2);
p_con=rho_con.*T_con;

%Calculation of other variables of interest
M_con=v_con./sqrt(T_con); %Mach Number
massflow_con=rho_con.*v_con.*A; %Non-dimensional Mass Flow Rate
p_con=rho_con.*T_con; %Non-dimensional pressure

%Transient Throat Values
rho_throat_con(k2)=rho_con(n_throat);
T_throat_con(k2)=T_con(n_throat);
p_throat_con(k2)=p_con(n_throat);
M_throat_con(k2)=M_con(n_throat);

%Check for time convergence
error2=max(abs(massflow_con-massflow_old_con));

%-----TRANSIENT MASS FLOW RATE PROFILE----------------
%We want to know how the massflow rate changes with time before the steady
%solution
TransientMassFlow_Con=figure(4);
hold on;

if k2 == 100
plot(x,massflow_con,'g');
else if k2 == 200
plot(x,massflow_con,'r');
else if k2 == 700
plot(x,massflow_con,'-.b','linewidth',1.25);
end
end
end
end

title('Transient Non-dimensional Mass Flow Rates(Non-Conservative Form)');
xlabel('x/L');
ylabel('\rhoVA/\rho_{o}V_{o}A_{o}');
legend('100\Deltat','200\Deltat','700\Deltat');
axis([0 3 0.5 0.8]);
grid minor;

end


3) NONCONSERVATION

3.1 - TIMEWISE VARIATION

3.1.1 - Throat Variables Change

This section shows the variation of the main properties at the throat during the iterations. The purpose of this section is to see the stability of the scheme, and to notice any strangely high or low results in the steady values. The non-dimensional values of the properties at the throat are well-known and that is why this section was chosen to check these calculations.

Here, the stability of the nonconservation form solver has been demosntrated. The slight disturbances during the first few iterations are soon dumped and the solution becomes stable. Although it will be compared lately, here we want to hihglight the coherence of the results, as the Mach=1 is achieved at the throat, as it was supposed to be.

3.1.2 - Mass Flow Profile Evolution Throuhg Time

Another useful parameter to check the solution convergence is the normalized mass flow rate. Under steady conditions, the mass flow rate of the nozzle should be constant at any of the nodes. This is actually the error parameter that was used to check for convergence in the time loop.

Here, note that the first dashed curve is the variation of the profile at the initial conditions. The strange looking is simply the product of the initial values given to the density, velocity and the parabolic variation of the nozzle area ratio. Note that after 50 times setps, the mass flow has considerably changed, and it keeps changing radically until around 200 steps. In here, we can see that the solution is approximating to the final solution. Around 700 iterations were enought to achieve the almost-constant mass flow rate profile. Also note that this convergention was to values really close to the analytical result of 0.579.

3.1.3 - Time and Iterations needed

Minimum Number of Iterations for Convergence (NonConservative Form): 767
Computation time (NonConservative Form): 3.275671 s


3.2.1 - Important variables Steadt Profile

Once the solver was tested to be stable, the final steady profiles were computed and are shown here.

What we can see in these graphs, is that they follow the tendency expected at a nozzle working under these conditions. The density, temperature and pressure are dropping through the exit of the nozzle, at the expend of an acceleration of the fluid (increase in the Mach Number). Although no compared here, these results are almost identical to the analytical solution of the problem, listed on the Reference book.

4) CONSERVATION

4.1 - TIMEWISE VARIATION

4.1.1 - Throat Variables Change

As in the nonconservation solver, this section was used to check for convergence in the numerical scheme. Note also that the change until the steady solutions does not go through the same path. That is because of the different numerical schemes (iteration solvers).

We can see that the properties are stable real soon in the simulation. The number of iterations is however higher, because the convergence condition was stated in the mass flow rate instead of these parameters. Again, the solution is stable and it seems to achieve promising results.

4.1.2 - Mass Flow Profile Evolution Throuhg Time

Once again, the mass flow rate profile provides us an extra insight to the mechanics of the timewise variation of the flow and its approach to the steady.

In here we can appreaciate that the solution converges in a softer way than with the nonconservation solver. Even from the beggining, the initial conditions lead to a more similar profile to the solution.

4.1.3 - TIME AND ITERATIONS TO CONVERGE

Minimum Number of Iterations for Convergence (Conservative Form): 951
Computation time (Conservative Form): 3.393511 s

In here, it can be seen that the computation time needed for the conservation solver is a bit  higher than to the nonconservation scheme. However, it does not seem like an important difference that would make one decide for one of them, things like accuracy and stability should be checked first.

Also note that the number of iterations is based on the mass flow rate criteria stablished, and not the primitive variables. Probably, changing to these last ones would make the conservation scheme to need less iterations (for what it can be seen in the previous section graphs)

4.2.1 - Important Variables Steady Profile

Once again, the solution achieved is not only stable, but also pretty accurate. This profile is almost identical to the once achieved in the nonconservation analysis

4.2.2 - Steady Mass Flow Profiles Comparisson

The stability of both schemes is undenniable. Although both seem to be reliable, it is now time to analyze accuracy and time, in order to decide if there is a more fitting scheme for this application.

In here, the steady mass flow profiles are compared. This figure illustrates a general advantage of the conservation form of the equations. This forms always does a better job of preserving math throughout the flow field. This is mainly because the mass flow itself is one of the variables at the equations. This does not mean that the conservation form is a better fit for the numerical scheme, as it will be seen in the following section.

5) GRID DEPENDENCE TEST

The previous discussion does not neccesarily mean the superiority of the conservation form, and further analysis must be done. In this path, the grid dependence test is usually critical to decide whether a problem can be solved faster with one numerical scheme or another (or if just can be solved with them).

This test consists on changing the number of grids on the mesh, to check if the results are dependent of the number of points. Any numerical scheme should has as an objective manage to find grid independence. In this case, we change the numer of grids from 31 to 61, and performed the same calculations. The result, showed in the table below, stated that the problem is already independent of the number of computation points taken in the control volume, at least to the precision needed here. Also, this results showed that the nonconservation solver is actually more accurate at computing the primitive variables.

 Density Ratio Temperature Ratio Pressure Ratio Mach Number Analytical Solution 0.634 0.833 0.528 1.0 Nonconservation (31 points) 0.639 0.836 0.534 0.999 Conservation (31 points) 0.648 0.839 0.544 0.987 Nonconservation (61 points) 0.637 0.835 0.532 0.999 Conservation (61 points) 0.644 0.838 0.540 0.989

### Autoignition delay - Cantera Jorge Martinez · 2018-12-11 03:00:54

The purpose of this project is to show the effect of initial temperature and pressure in the ignition delay in a combustion process. For this purpose, a Constant Volume Reactor and the corresponding network objects are created in Cantera, and the combustion process is s Read more

### Shock Tube Simulation using Converge Jorge Martinez · 2018-11-27 02:11:31

The purpose of this project, is to setup and analyze a transient shock tube simulation using Converge Studio. The project will analyze how the interaction of the original shock waves and the reflected waves affect the fluid in the tube, and how the computational mesh ha Read more

### Conjugate Heat Transfer Simulation Jorge Martinez · 2018-10-23 03:33:38

The purpose of this project is to simulate the conjugate heat transfer in a pipe using Converge Studio. For this, an air flow is simulated through an aluminum pipe, which is receiveing a constant heat flux of 10,000 W/m2 at the outside wall. The objective is to understa Read more

### Prandtl-Meyer Shock Problem - Converge Studio Jorge Martinez · 2018-10-20 03:51:53

This project objective is to simulate the Prandtl Meyer shock wave phenomena using Converge Studio. First, a quick literature review about shock waves and their boundary conditions is provided. Then , the problem is set up and solved, focusing in how the different param Read more

### Multivariate Newton-Raphson method Jorge Martinez · 2018-10-18 00:46:05

The purpose of this assignment is to create a Python program including a Multivariate Newton Rhapson Solver, to solve a non-linear coupled differential system. This method is really useful for stiff systems, where the explicit solver are unstable. The system of equatio Read more

### Simulation of Flow Over a Throttle Body Converge Studio - Transient Jorge Martinez · 2018-10-13 22:13:18

The objective of this part of the project, is to simulate the transient flow through an elbow pipe with a throttle valve in the middle. This valve will be rotating from zero to 25 degrees, and then stay steady.  Case Setup. In the steady state part of the project Read more

### Simulation of Flow Over a Throttle Body - Steady State Jorge Martinez · 2018-10-12 03:30:10

This project purpose is to simulate the flow over an elbow body that contains a throttle using CONVERGE Studio. In this first part of the project, an steady state analysis will be set up. This means that the throttle will not move. This case helps set up the real analys Read more

### Heat Recuperator - Combustion Efficiency Jorge Martinez · 2018-10-10 05:34:46

In this challenge, the efficiency of using a heat recuperator to pre-heat the air is analyzed. For that, we simulate a furnace burning methane and air mixture, and we analyze first how the pre-heating affects at the AFL, and then the total heat/work that could be extrac Read more

### Adiabatic Flame Temperature using Python and Cantera - Reactor with Heat Loss Jorge Martinez · 2018-10-09 02:16:08

This project main objective is to calculate the Adiabatic Flame Temperature using Python and Cantera. In the first part of the project, the effect of the equivalence ratio for the methane combustion in a constant volume chamber is analyzed, and the python and cantera re Read more

### Backward Facing Step Flow - CONVERGE Studio Jorge Martinez · 2018-10-03 02:52:16

This projects aims to simulate a backward facing step flow using CONVERGE Studio. The flow is creating when a not reacting specie, air (23% O2 and 77% N2 in mass fraction), suffers a pressure gradient between the ends of the channel. The velocity, pressure and recirulat Read more