Showing posts with label Python. Show all posts
Showing posts with label Python. Show all posts

Thursday, 25 August 2011

I am looking for a job

News - I greatly appreciate your comments below and surprisingly plenty of information and valuable advices came to me during the last one year. Thanks to all your concern, I finally accepted a thermal expert position in ABB USA (www.abb.com) and am sitting in North Carolina.

I particularly thank the three companies in England, France and America that I personally consulted for and their encouraging recommendations. I cherish the experiences I have obtained from the part-time working with them.

Within my PhD study, I have published eight research papers. Four of them are in leading journals and their links are as follows:

I am graduating in September 2011 with a PhD degree from School of Electrical and Electronic Engineering, University of Manchester. I am looking for a numerical modelling job.

I spent the past more than three years on the CFD simulations which are coupling between electric, thermal and fluid dynamic fields, from which I have accumulated rich experiences on modelling packages such as SALOME, Code_Saturne, Elmer, ParaView and COMSOL etc. I share experiences in this blog, which has also gained worldwide popularity.

I have been programming in C/C++, C#, VB.NET and Matlab since 2002. I also have experiences in Python (the script widely used for automation of modelling and the software such as SALOME and ParaView etc) and Fortran (the user routine definition for Code_Saturne) etc. Meanwhile, probably as a bonus, I am a fan of open-source software, Linux and Mac OS X. I even wrote device drivers for Linux and ported the MySQL client tool running on desktops onto embedded Linux platforms.

I bear high interests in numerical modelling and would like to pursue a career in this field. If you have relevant opportunities, please drop me an email (salad00 at gmail dot com). I really appreciate your help very much.

Saturday, 6 February 2010

CFD tutorial: laminar flow along a 2D channel - Part I

When fluid flows along a 2D channel, which is between two infinite parallel plates, with a relatively low velocity (laminar flow), a hydraulic stable status can be achieved after a certain distance (entrance length). This stable flow status is named as hydraulically fully developed flow. Furthermore, if there is heat flux at the duct surfaces, the fluid near the wall will be warmed up, and a thermal boundary layer is then formed.

For a flow along a channel without heated, after it is fully developed, the velocity pattern at a cross area is parabolic. However, if the flow is also heated, the parabolic velocity pattern will be twisted because of the viscosity dependency of the fluid (this is true for many types of liquid). This phenomenon is pictorially shown as the figure below, in which (a) is the parabolic velocity profile, and (b) is the distorted profile in a heated duct due to the temperature dependency of the fluid viscosity.


Denote the duct height as H, duct length as L, fluid density as , and the average velocity in the duct region as U. Then the pressure drop between the two ends of the duct can be expressed as



in which is defined as friction coefficient. Actually this equation says the active pressure force has to balance the friction force due to shear stress, which is finally related to the fluid viscosity.

Anyway, I select this example to show a simple CFD calculation, because its simple geometry. Using SALOME we can very quickly build up a geometry and then mesh it. Probably an expected structured mesh of the channel can be like this. It is better to have dense mesh near the inlet and the walls of the duct.



Geometry and mesh using SALOME

According to the philosophy of executing commands on terminals, I'll use python scripts to re-perform the procedure of modelling in SALOME.

#######################################################################
# Geometry construction and meshing creation for a typical
# 2d channel flow between two infinite parallel plates.
#
# Written by: salad
# Manchester, UK
# 06/02/2010
#######################################################################

How to execute these scripts? New a blank study in SALOME and then a 'Python Console' can be seen at the bottom of the window. Input or copy scripts into the console and press Enter to execute them. Scripts can also be saved in a .py file, and with the menu item 'File >> Load Script...' a dialog 'Load python script' can help load script files and execute them as a batch. 'Ctrl+T' is the shortcut to call the dialog.

Ok, let's start. 1. Import the necessary libraries and define the related dimensions. We are going to build a duct with height 5 mm and length 100 mm.

import os
import sys
import salome
from math import *
from geompy import *
import copy

# basic unit
unit = 0.001
# extrusion length
extru = unit

duct_height = unit * 5
duct_length = unit * 100

'extru' is the extrusion length along the z axis. This length will be meshed with only one cell. Actually Code_Saturne only handles 3D meshes with FVM. Therefore, to simulate a 2D mesh, it is a good practice to extrude the 2D mesh with a 1 cell layer thickness.

2. Build the basic coordinate.

#
# basic coordinate
#
p0 = MakeVertex(0, 0, 0)
dx = MakeVectorDXDYDZ(unit, 0, 0)
dy = MakeVectorDXDYDZ(0, unit, 0)
dz = MakeVectorDXDYDZ(0, 0, unit)
addToStudy(p0, "p0")
addToStudy(dx, "dx")
addToStudy(dy, "dy")
addToStudy(dz, "dz")

print "basic coordinate built..."

3. Construct the channel geometry. The geometry modelling procedure can be briefly expressed as: points => edges => faces => extrusion. Points are constructed with coordinates, edges are defined by linking points, and faces are then defined by specifying its edges. Finally extrusion helps convert the 2D geometry to 3D.

#
# geometry construction
#

# points
p1 = MakeVertex(duct_length, 0, 0)
p2 = MakeVertex(0, duct_height, 0)
p3 = MakeVertex(duct_length, duct_height, 0)

# build edges
e0 = MakeEdge(p0, p1)
e1 = MakeEdge(p2, p3)
e2 = MakeEdge(p0, p2)
e3 = MakeEdge(p1, p3)

# build face
f_duct = MakeFaceWires([e0, e1, e2, e3], 1)

# extrude along the z axis
v_duct = MakePrismVecH(f_duct, dz, extru)
addToStudy(v_duct, "v_duct")

print "channel geometry constructed..."

4. Define the groups. Groups are necessary because they correspond to the boundaries used during CFD calculations. Here we define the geometry groups and these groups will then be encapsulated into mesh groups later. Mesh groups will be used to define boundary conditions when using Code_Saturne.

Firstly, define the face groups. Functions used in the python script can be looked up in the SALOME document.

# FACE groups definition:
# 1. inlet
# 2. outlet
# 3. bottom
# 4. top
# 5. sym

# sym is declared and filled
sub_faces = SubShapeAllSorted(v_duct, ShapeType["FACE"])
g_sym = CreateGroup(v_duct, ShapeType["FACE"])
UnionList(g_sym, sub_faces)

# inlet
sub_faces = GetShapesOnPlane(v_duct, ShapeType["FACE"], dx, GEOM.ST_ON)

g_inlet = CreateGroup(v_duct, ShapeType["FACE"])
UnionList(g_inlet, sub_faces)
DifferenceList(g_sym, sub_faces)

addToStudyInFather(v_duct, g_inlet, "inlet")

# outlet
sub_faces = GetShapesOnPlaneWithLocation(v_duct, ShapeType["FACE"], dx, p1, GEOM.ST_ON)

g_outlet = CreateGroup(v_duct, ShapeType["FACE"])
UnionList(g_outlet, sub_faces)
DifferenceList(g_sym, sub_faces)

addToStudyInFather(v_duct, g_outlet, "outlet")

# bottom
sub_faces = GetShapesOnPlane(v_duct, ShapeType["FACE"], dy, GEOM.ST_ON)

g_bottom = CreateGroup(v_duct, ShapeType["FACE"])
UnionList(g_bottom, sub_faces)
DifferenceList(g_sym, sub_faces)

addToStudyInFather(v_duct, g_bottom, "bottom")

# top
sub_faces = GetShapesOnPlaneWithLocation(v_duct, ShapeType["FACE"], dy, p2, GEOM.ST_ON)

g_top = CreateGroup(v_duct, ShapeType["FACE"])
UnionList(g_top, sub_faces)
DifferenceList(g_sym, sub_faces)

addToStudyInFather(v_duct, g_top, "top")

# sym is finally obtained
addToStudyInFather(v_duct, g_sym, "sym")

print "FACE groups defined..."

Then define the edge groups.

# EDGE groups definition:
# 1. let (inlet & outlet of the channel)
# 2. tb (top & bottom of the channel)
# 3. extru (the extrusion length)

# let
g_let = GetEdgesByLength(v_duct, duct_height, duct_height, 1, 1)
addToStudyInFather(v_duct, g_let, "let")

# tb
g_tb = GetEdgesByLength(v_duct, duct_length, duct_length, 1, 1)
addToStudyInFather(v_duct, g_tb, "tb")

# extru
g_extru = GetEdgesByLength(v_duct, extru, extru, 1, 1)
addToStudyInFather(v_duct, g_extru, "extru")

print "EDGE groups defined..."

5. Mesh. Till now the geometry is already produced. Meshing can be performed with the script below. Note that the expressions, which define the corresponding mesh density profiles, can be edited to any others if necessary.

#
# meshing
#
import smesh
import StdMeshers

mesh_d = smesh.Mesh(v_duct, "mesh_d")

print "prepare for meshing..."

# construction of mesh

# as default, a 1D edge is meshed with only 1 cell
# it is for the extrusion length
algo1d = mesh_d.Segment()
algo1d.NumberOfSegments(1)

# structured rectangular mesh is preferred for 2D faces
algo2d = mesh_d.Quadrangle()
algo2d.QuadranglePreference()

algo3d = mesh_d.Hexahedron()

# submesh

# for inlet & outlet, use parabolic mesh density profile
algo1d_let = mesh_d.Segment(g_let)
seg = algo1d_let.NumberOfSegments(50)
seg.SetDistrType(3)
seg.SetConversionMode(1)
seg.SetExpressionFunction('(t-0.5)^2+0.1')

# for top & bottom, use a decreasing profile from the inlet to the outlet
algo1d_tb = mesh_d.Segment(g_tb)
seg = algo1d_tb.NumberOfSegments(250)
seg.SetDistrType(3)
seg.SetConversionMode(0)
seg.SetExpressionFunction('(t-1)^4+0.1')

status = mesh_d.Compute()

print "mesh computed..."

The same with the geometry modelling, groups of mesh are also defined accordingly.

#
# mesh groups
#

mesh_d.GroupOnGeom(g_inlet, "inlet")
mesh_d.GroupOnGeom(g_outlet, "outlet")
mesh_d.GroupOnGeom(g_bottom, "bottom")
mesh_d.GroupOnGeom(g_top, "top")
mesh_d.GroupOnGeom(g_sym, "sym")

print "mesh groups defined..."

Arriving at this, you should be able to update the SALOME GUI to see the mesh generated, and also judge whether it is what you want.

# update
salome.sg.updateObjBrowser(1)

6. Output. After the mesh is produced, it has to be saved into a file before using it with Code_Saturne. In SALOME, right-click the object 'mesh_d' and from the context menu we can see that SALOME supports exporting to MED and UNV file formats. If your Code_Saturne was compiled with MED support, I recommend to use MED because of its smaller file size; or UNV can be used.


Finally a file 'mesh_d.med' (or 'mesh_d.unv') is obtained. This is the mesh file we will move to Code_Saturne's MESH directory.

CFD solving with Code_Saturne

Please read Part II of this article.

Post-processing with ParaView

Please read Part III of this article.