Saturday, 14 December 2024

Temperature Data Extraction, Visualization and Analysis using Scilab

Topic: Numerical Analysis of Temperature from Data Acquisition System (DAS);

Subject: Heat Transfer;

Tool: Scilab;

By: Gani Comia, Dec. 2024;

Data acquisition devices or systems record process data over time for monitoring and analysis purposes. This article presented how the Scilab software extracts, visualizes, and analyzes temperature numerical data from a DAS.

While DAS can actually show a graph of the data in real-time, other analysis can only be done after downloading the data soft copy in .CSV format. This blog post presented the two analyses made with the use of Scilab scripting. The first analysis made was to find the specific time for a particular temperature value, and the second one was the calculation of heat rate, or the temperature change per unit time during ramp-up.

The completed plot or graph of numerical data using Scilab is shown.

Fig. 1. Scilab Plot and Analysis of Temperature Data.

Below are the code snippets or scripts for the two sample analyses made. The completed script is shown after.

Analysis (1): Finding the Time for Specific Temperature.

// ----- analysis(1) : finding time @ temp = 232 degC -----
modOvenWallTemp = round(ovenWallTemp);
id = find(modOvenWallTemp >= 232,[-1]);  // max no. of indices
timeStd(1) = time(id(1));
timeStd(2) = time(id($));

Analysis (2): Temperature per Unit Time during Ramping-up.

// ----- analysis(2) : finding heat rate in degC / min -----
tempCond = [ovenWallTemp(1) temp(2)];
timeCond = [time(1) timeStd(1)];
heatRate = diff(tempCond)./diff(timeCond);


Scilab Script (Complete):

// Copyright (C) 2024 - Gani Comia
// Date of creation: 4 Dec 2024
// Heating Oven Temperature Data Analysis using Scilab
clear;clc;
// data extraction
clear importdata;

function [data]=importdata(filename)
    data = csvRead(filename, ",", ".", "double")
endfunction

[A] = importdata("temp_profile_before_repair_xform.csv");

// assigning data to the variables
time = A(2:$,13);
ovenWallTemp = A(2:$,2);
partTemp = A(2:$,9);

// data visualization
clf;
f=gcf();
f.figure_size=[900,600];
plot(time,ovenWallTemp,"-b","linewidth",1.75)
plot(time,partTemp,"-.r","linewidth",1)
legend(["Oven Wall","Workpiece"],"in_upper_right")
title("Analysis of Oven Temperature Profile from Data Logger")
xlabel("Time [ min ]")
ylabel("Temperature [ deg C ]")
xstring(245,100,"https://gani-mech-toolbox.blogspot.com")
a=gca(); 
a.data_bounds = [0 0;400 300]; 

// additional plot properties
t = [0:60:360, 400];
temp = [0 232 300];

// temp = 232 degC horizontal lines
plot([t(1),t($)],[temp(2),temp(2)],"--g","linewidth",1)
xstring(5,232,"Std = 232 C")

// hourly verical lines
for i = 2:length(t)-1
    plot([t(i),t(i)],[temp(1),temp(3)],"--y","linewidth",1)
end

// hourly label 
minute = 60;
hour = 1;
for i = 2:length(t)-1
    note = [string(hour),"hour"];
    xstring(minute,10,note,-90);
    minute = minute + 60;
    hour = hour + 1;
end

// ----- analysis(1) : finding time @ temp = 232 degC -----
modOvenWallTemp = round(ovenWallTemp);
id = find(modOvenWallTemp >= 232,[-1]);  // max no. of indices
timeStd(1) = time(id(1));
timeStd(2) = time(id($));

// plot of analysis(1)
timeStd = [timeStd(1) timeStd(2)];
stdTemp = [232 232];
plot(timeStd, stdTemp, "ro")
xstring(timeStd(1),stdTemp(1)-15,[string(timeStd(1)),"min"])
xstring(timeStd(2),stdTemp(2),[string(timeStd(2)),"min"])
ambient = ["Ini =",string(round(ovenWallTemp(1))),"C"];
xstring(0,ovenWallTemp(1)-15,ambient)

// ----- analysis(2) : finding heat rate in degC / min -----
tempCond = [ovenWallTemp(1) temp(2)];
timeCond = [time(1) timeStd(1)];
heatRate = diff(tempCond)./diff(timeCond);

// plot of analysis(2)
plot(timeCond,tempCond,"--","color","dimgray","linewidth",2)
noteHR = ["RATE:",string(round(heatRate)),"deg C / min"]
xstring(t(2)-20,110,noteHR,-58)

A readily available program script with some minor revision can facilitate visualization and analysis of numerical data for almost similar conditions. Program or code reusability leading to short lead times of engineering analysis is one of the benefits of using a programming tools in the engineering field.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.


Friday, 13 December 2024

Steady-State Heat Conduction Model using FDM Solution to 2D Laplace Equation

Topic: Finite Difference Method of Laplace Equation in 2D;

Subject: Heat Transfer and Numerical Methods;

Tool: Scilab;

By: Gani Comia, Dec. 2024;

The steady-state heat conduction model can be represented mathematically by Laplace equation. Shown is the Laplace equation in 2D:

$$\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\;=\;0$$

The solution, \(T(x,y)\), represents the temperature distribution. The Laplace equation can be discretized using the Finite Difference Method (FDM) to approximate the solution numerically. This article will present the Scilab script of the FDM solution to the Laplace equation.

The FDM solution for a uniform grid, \((\Delta x = \Delta y)\), with the boundary temperatures specified based on Dirichlet conditions, is shown:

$$T_{i,j}\;=\;\frac{T_{i+1,j}+T_{i-1,j}+T_{i,j+1}+T_{i,j-1}}{4}$$

The final value of  \(T(x,y)\) is the result of iteration based on the convergence criterion:

$$max\left\vert T_{i,j}^{new}\;-\;T_{i,j}^{old}\right\vert\;<\;\epsilon$$

Consider a sample 2D domain of dimensions 30x30 units with the following constant temperature conditions. 

Fig. 1. Sample 2D Domain with Constant Temperature.

Shown below is the temperature profile as represented by \(T(x,y)\) using the FDM solution to the Laplace equation.

Fig. 2. Temperature Profile \(T(x,y)\).

The Scilab script as a toolbox can be used to implement the FDM numerical solution for the steady-state heat conduction model represented mathematically using the Laplace equation.

Scilab Script

// Copyright (C) 2024 - Gani Comia
// Date of creation: 8 Dec 2024
// FDM Solution to 2D Laplace Eq, d2T/dx2 + d2T/dy2 = 0
clear;clc;
// domain and grid parameters
Lx = 30;                    // units, domain length at x-axis
Ly = 30;                    // units, domain length at y-axis
dx = 1;                     // grid spacing at x-axis
dy = 1;                     // grid spacing at y-axis
nx = Lx./dx;                // no. of grid points in x-direction
ny = Ly./dy;                // no. of grid points in y-direction
max_iter = 1000;            // max no. of iterations
tolerance = 1e-6;           // convergence tol.

// temperature field initial value
T = zeros(nx, ny);
// boundary value conditions
    T(:,1) = 0;             // T = 0 C @ y = 1
T(1,:) = 0;                 // T = 0 C @ x = 1
T(:,ny) = 0;                // T = 0 deg C @ y = 30
T(nx,:) = 100;              // T = 100 deg C @ x = 30

// finite difference method (Jacobi iteration)
for iter = 1:max_iter
    Tnew = T;
    for i = 2:nx-1
        for j = 2:ny-1
            Tnew(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1));
        end
    end
    // convergence check
    if norm(Tnew - T) < tolerance then
        break;
    end
    T = Tnew;
end
// calculation status 
mprintf("\n\tCalculation completed!\n")

// plot properties
clf;
f = gcf()
f.figure_size = [700,650]
ax = gca()
ax.tight_limits = ["on" "on" "off"];

// plotting results
surf(T);
title(["$Temperature\;Profile$","$T(x,y)$"],"fontsize",4);
xlabel("x-axis");
ylabel("y-axis");
zlabel("$T$");
colormap(jet);
colorbar(0,100);

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.

Saturday, 7 December 2024

Curve Fitting with Polynomials using Python

Topic: Method of Curve Fitting with Polynomial Functions;

Subject: Numerical Methods;

Tool: Python;

By: Gani Comia, Dec. 2024;

Curve fitting with polynomials finds a polynomial equation of a given degree that approximates the set of data points. Curve fitting is used in the following applications:

  • Data Smoothing
  • Predictive Modelling
  • Scientific Analysis

The general form of a polynomial function of degree n is:

$$P(x)\;=\;{a_n\,x^n}+{a_{n-1}\,x^{n-1}}+\ldots+{a_1\,x}+{a_0}$$

This article shows a Python script as toolbox to implement the method of curve fitting with polynomials. The script requires numerical data for both dependent and independent variables. For a given sample of data points in Fig.1., an approximate polynomial curve of 3rd degree shown in Fig.2. can be calculated and plotted using Python script.

$$xdata = [ 0 \quad 1 \quad 2 \quad 3 \quad 4 \quad 5 \quad 6 \quad 7 ]$$

$$ydata = [ 5 \quad 10 \quad 60 \quad 50 \quad 130 \quad 220 \quad 300 \quad 280 ]$$

Fig.1. Sample Data Points for Approximating Polynomial Curve.

Below is the approximate polynomial curve that closely fit the given data points.

Fig.2. Polynomial Curve of 3rd Degree from Data Points.

Python Script

# -*- coding: utf-8 -*-
"""
Created on Mon Nov  7 21:08:48 2024
@author: Gani Comia
"""
'''
Curve Fitting with Polynomials for a given xdata and ydata
'''
import numpy as np
import matplotlib.pyplot as plt

# given: xdata and ydata
xdata = np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7.])
ydata = np.array([5. , 10. , 60. , 50. , 130. , 220. , 300. , 280.])

# solution:
# solving for coefficients of polynomial fit for cubic (degree=3)
# z is an array of coefficients , highest first , i.e.
# z = np.array ([c_3, c_2 , c_1 , c_0])
z = np.polyfit(xdata, ydata, 3)
print("\n",z,"\n") 

# ‘poly1d‘ returns the polynomial function from coefficients    
p = np.poly1d(z)
print(p)

# polynomial's curve data points
xs = [0.1 * i for i in range (70)]  # from 0 to 7, step=0.1
ys = [p(x) for x in xs]   # evaluate p(x) for all x in list xs

# creating plot of data points and polynomial curve for visualization
plt.plot(xdata, ydata, 'o', label='data points') # data points
plt.plot(xs, ys, label='fitted curve')    # polynomial fit
plt.title('Curve Fitting with Polynomials')
plt.ylabel('y-value')
plt.xlabel('x-value')
plt.legend()
plt.grid()
plt.show()

On the above script, Python Numpy’s utilities, np.polyfit(), returns the polynomial coefficients of the given degree that approximately fits on the given data points.  Numpy’s np.poly1d() creates a polynomial functions for plotting and calculation.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the softcopy of the file.


C++ and Python for Numerical Solution to Spring-Mass-Damper Model

Topic: SDOF Spring-Mass-Damper System Subject: Numerical Methods & Mechanical Vibration Tools: C++, CodeBlocks, GitHub Copilot, LabP...