Structured 3D Solid Mesh of Complex Thin Parts Using Dual Kriging Interpolation
H. Champliaud, F. Duchaine and V.N. Lê
Mechanical Engineering Department
École de Technologie Supérieure
Montreal (Quebec), H3C 1K3, Canada
To reduce computer time and required memories in finite element (F.E.) analysis of complex thin parts with 3D
solid elements, a structured brick type element meshing technique is proposed using dual kriging interpolation.
A program using this technique is written in MatLab® language to create data for F.E. models into files easily
readable by Ansys® F.E. code. Results show a reduction by factors of about 2 on number of nodes and up to 6
on number of elements compared to tetrahedral type models in meshing sheet metal shape.
Keywords: Forming, Sheet, Interpolation, Dual Kriging, Solid mesh, Brick type, Tetrahedral
One of major problems with finite element analysis in metal forming problem is the high number of nodes due
to meshing complex volumes into tetrahedral elements [1, 2]. This is mostly true when dealing with large
amount of computer time and disk space in iterative solution procedures for non-linear behaviors. The present
work deals with a structured way of meshing complex thin parts into brick type finite elements so as to reduce
model size, computer time and required disk space, accordingly. In the first part of the paper the dual kriging
methodology is recalled. Second, the steps of the 3D meshing MatLab® program are discussed. Finally,
examples of applications of the program are shown.
2 Dual kriging interpolation
Kriging is a generalization of the interpolation least squares method. The main interest of the method is that it
allows the fitting model to pass through the data points. The name of kriging is due to Matheron  who first
present in a rigorous mathematical form the work of Krige  who originally applied the method for evaluation
in mining exploitation.
In this paper the surface of complex thin parts is described by three coordinate functions x(Xi), y(Xi) and z(Xi)
which are based on dual kriging interpolation in such a way that they fit all the three dimensional sample points
Xi. Each function is decomposed into the sum of two terms, in the form a(Xi) + b(Xi), where a(Xi) represents
the average behavior, called the drift, and b(Xi) is an error term called the fluctuation. The Xi are the location
sample points recorded in a structured sequence on the surfaces of the part. The result is a parametric grid that
allows to mesh the part into shell or brick elements. The observations ui of a phenomena, at data points Xi along
a row, can be interpolated using the basic dual kriging model written in the form [5, 6]:
u(X) = ∑aipi(X) + ∑b jK(X−X j )
i =1 j=1
where X is the 3D point where the phenomena is evaluated, ai one coefficient of the M functions pi(X) forming
the drift, and bj one coefficient of the N functions K(|X – Xj|) forming the fluctuation (also called the
generalized covariance). |X – Xi| is the Euclidean distance between X and the data point Xi. M is the number of
functions needed to form the drift, and N is the number of the data points along a row (or a column) of data
points (see figure 1).
* 3D recorded Xi data points
* * P(s,t) *
* ¤ Kriging interpolation by row
* * Kriging interpolation by column
y * * *
Figure 1: Kriging interpolation of a surface
By using the normalized coordinates s and t for writing x, y and z along rows and columns, two parametric
equations for each coordinate are obtained by dual kriging. For example
for a column k: xk(s) = ∑asipi(s) + ∑bs jK(s−s j )
i =1 j=1
for a row m: xm(t) = ∑atipi(t) + ∑bt jK( t −t j )
i =1 j=1
where k = 1 to number of columns, and m = 1 to number of rows. The coefficients asi, bsj, ati and btj are
obtained by solving the equations xmk = xk(sm) = xm(tk) at each data point. Finally it can be shown that:
x(s,t) = [xk(s)]T[xm(t)]
and similar functions for y(s,t) and z(s,t). These functions are then used to construct a regular grid on the surface
with any desired refinement in s and t directions.
Three dimensional data of thin parts are recorded in files using a laser detection system or any coordinate
measuring device, see figure 2. The whole surface of the part is scanned by rows, with the same number of
samples points per row. This is a necessary criteria for structuring the mesh with only quadrilateral or brick
elements. If the numbers of samples points are different between rows, structured mesh can still be done such as
explained in step 2.
thin part controller
Figure 2: Example of recording data points with a coordinate measuring machine (CMM)
After copying the generated file in the database, the Matlab® program is run. A dialog box is opened allowing
the user to select data file to load. The kriging interpolation window is displayed next to provide some options
on fitting the surface, see figure 3. The recorded data are then interpolated row by row and column by column,
using the dual kriging method, in order to describe the complete surface. To obtain a regular grid over the
surface it is necessary that the number of data points be the same along each recorded sequence (row or
column). If it is not the case a redistribution of points is automatically done by kriging rows or columns
individually for filling data where sample points are missing.
Figure 3: User interface offering the kriging choices
If solid mesh option is selected in step 2, the normal vector at each point of the grid is computed, and a new
point is generated along this vector at a distance equal to the input thickness. The first surface can be assigned
inside, outside or middle (see figure 4). Polynomial functions are often used to construct the drift and the
fluctuation. These functions are useful to fit smooth surfaces for they are continuously differentiable and easily
handled to compute the normal unit vector at any point by (see figure 1):
∂OP ∧ ∂OP
n = ∂sr ∂t r
Figure 4: User interface offering the thickness options
When the fitting model is accepted, it can be meshed. The refinement of the mesh can be controlled in each
direction, depending of the desired precision of analysis results. The mesh options are presented in figure 5.
Figure 5: User interface for mesh options
Figure 6: Resulting mesh in Ansys®
Users can go back and forth among dialog boxes to change values and options until satisfactory mesh is
obtained. The final results are written into three files readable by Ansys® finite element code. These files are
node coordinate file, element connectivity file, and preprocessing set up log file. After reading these files into
Ansys®, the remaining steps to be completed are material properties, boundary conditions, loading conditions,
etc. An example of brick type mesh read by Ansys® from log file created by kriging interpolation program is
shown in figure 6. Actually, only log files for Ansys are generated, but it is not much work to create log files for
other codes such as Nastran®, Algor®, Catia®, Proengineer®, …
4 Meshing thin parts
Examples of dual kriging of a complex part is presented in figure 7. The part is mesh into brick elements with
desired patterns. The results show that the model size could be reduced by at least a factor of 2 on number of
nodes and up to 6 on number of elements compared to equivalent edge size tetrahedral finite element models.
The structured mesh program is written in Matlab® and then converted into a stand-alone application. The
program writes also necessary preprocessing commands in a log file for Ansys® finite element code, including
element attributes such as type and material.
Figure 7: Comparison between a tetrahedral mesh and a brick mesh
The methodology presented here is a meshing technique using the dual kriging interpolation. The methodology
is programmed with MatLab® which offers a powerful environment for matrix computations and graphics. The
program generates log files for preprocessing finite element analysis. The size of the problem are reduced by at
least a factor 2 in number of nodes, and up to a factor 6 in number of element.
Authors want to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for
financial support, and A. Croteau and M. Cheriet of the Automated Production Engineering Department, École
de Technologie Supérieure, for valuable helps.
1. Champliaud H., Analyse par éléments finis du sertissage de capsules d’étanchéité, Ph. D. dissertation,
Mechanical Engineering Department, École de Technologie Supérieure, (2000).
2. Champliaud H., Lê N.V., Finite Element Analysis of Crowning Sealing Caps, 2000 Ansys Conference,
3. Matheron G., The Intrinsic Random Functions and their Applications, Adv. Appl. Prob., 5, 439-468, (1973).
4. Krige, D.G., A Statistical Approach to some Basic Mine Valuation Problems on the Witwatersrand, J.
Chem. Metall. Min. Soc. S. Afr., 52, 119-139, (1951).
5. Trochu F., A Contouring Program Based on Dual Kriging Interpolation, Engineering with computers, vol.
9, pp.160-177, (1993).
6. Poirier, C., Tinawi,R. Finite Element Stress Tensor Field Interpolation and Manipulation using 3D Dual
Kriging. Computer & Structures, vol.40, No.2, pp. 211-225, (1991).