# To find slope of Force Residual plot¶

In Computational Fluid Dynamics simulations, to conclude whether the simulation is converged or not, we usually look at the residual values.
Looking at the Force plot residual, if the plot has reached a stabilized state (values flatlined), we can conclude simulation has converged, i.e. there won't be any change in averaged flow field with time.
For unsteady case, Since this is unstedy flow case, the residual plot of Total force will be oscillating around an average value.
But this entire oscillation will be settling down with the timesteps.
This is python code to check the residual plot if it has reached that stabilized state.
After ommiting initial transient stage and considered last time steps for this plot.
For this used Numpy polyfit(x,y,deg) function for finding best fit of polynomial of degree 'deg' to a set of data points given by the array of arguments x and y.
Also used utility poly1d to plot the above found polynomial.

In [4]:
import numpy as np
f=open('FrcMnt.out')

In [5]:
#Initiating lists
x=[]
y=[]

# Reading timestep and total force
for i in lines:
lines1=" ".join(i.split())
lines2=lines1.split(' ')
x.append(lines2[1])
y.append(lines2[11])

In [6]:
'Length of array {}'.format(len(x))

Out[6]:
'Length of array 9963'
In [7]:
# Converting from string to float

x=[float(i) for i in x]
y=[float(i) for i in y]

In [8]:
import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(x,y, 'ro')
fig.suptitle('Residual Plot', fontsize=20)
plt.xlabel('Timesteps', fontsize=18)
plt.ylabel('Total Force', fontsize=16)
plt.show()

<matplotlib.figure.Figure at 0x7f7b01fa57b8>

# Averaging after timestep, n¶

In [9]:
n=6000

In [10]:
x1=x[n:]
y1=y[n:]

In [11]:
# Using least squres function for curve fitting
fit = np.polyfit(x1,y1,1)
fit_fn = np.poly1d(fit)

# fit_fn is now a function which takes in x and returns an estimate for y

fig = plt.figure()
plt.plot(x1,y1, 'yo', x, fit_fn(x), '--k')
plt.xlim(n, 9500)
fig.suptitle('Residual Plot with Slope Line', fontsize=20)
plt.xlabel('Timesteps', fontsize=18)
plt.ylabel('Total Force', fontsize=16)
plt.show()

In [12]:
def slope(a,b):
return (fit_fn(x)[a]-fit_fn(x)[b])/(x[a]-x[b])


# Slope of the fitted line¶

In [13]:
print(round(slope(0,1),5))

-0.02118


Since this is unstedy flow case, the residual plot of Total force will be oscillating around an average value.
But this entire oscillation will be settling down with the timesteps.
Since the slope of line is small value ~0, we can say flow has stabilized.
Though the slope is small value, the line orientation looks steep as the y axis values has very high magnitude

Below cells are for styling, you can happily ignore them

In [2]:
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))