UDK 519.61/.64:539.3
Original scientific article/Izvirni znanstveni članek
ISSN 1580-2949 MTAEC9, 47(6)789(2013)
NON-SINGULAR METHOD OF FUNDAMENTAL SOLUTIONS FOR THE DEFORMATION OF TWO-DIMENSIONAL ELASTIC BODIES IN CONTACT
NESINGULARNA METODA FUNDAMENTALNIH REŠITEV ZA DEFORMACIJO DVO-DIMENZIONALNIH STIKAJOČIH SE ELASTIČNIH TELES
Qingguo Liu1, Božidar Šarler1,2,3
!University of Nova Gorica, Nova Gorica, Slovenia 2IMT, Ljubljana, Slovenia 3Center of Excellence BIK, Solkan, Slovenia qingguo.liu@ung.si, bozidar.sarler@imt.si
Prejem rokopisa — received: 2013-03-11; sprejem za objavo - accepted for publication: 2013-10-10
The development of an effective new numerical method for the simulation of the micromechanics of multi-grain systems in contact is developed in the present paper. The method is based on the Method of Fundamental Solutions (MFS) for two-dimensional plane strain isotropic elasticity and employs the Kelvin Fundamental Solution (FS). The main drawback of MFS is the presence of an artificial boundary, outside the physical boundary, for positioning the source points of the FS, which is difficult or impossible in multi-body problems. In order to remove the singularities of the FS the point sources are replaced by the distributed sources over circular disks. The values of the distributed sources are calculated in a closed form in the case of the Dirichlet boundary conditions. In the case of the Neumann boundary conditions the respective values of the derivatives of the FS are calculated indirectly from the considerations of the solution of simple displacement fields. A problem of two, four and nine bodies in contact is tackled. The newly developed method is verified based on a comparison with the classic MFS. The numerical method will form a part of the microstructure-deformation model, coupled with the macroscopic thermo-mechanics simulation system for continuous casting, hot rolling and heat treatment.
Keywords: isotropic elasticity, plane strain, Navier's equation, displacement and traction boundary conditions, non-singular method of fundamental solutions, Kelvin's fundamental solution
V članku opisujemo učinkovito novo numerično metodo za simulacijo mikromehanike sistemov z več zrni v stiku. Ta temelji na metodi fundamentalnih rešitev (MFR) za dvo-dimenzionalne probleme izotropnih elastičnih ravninskih deformacij in uporablja Kelvinovo fundamentalno rešitev (FR). Bistvena slabost MFR je prisotnost fiktivnega roba zunaj fizikalnega roba za postavitev izvirnih točk FR, kar je težavno ali nemogoče pri problemih z več telesi. Z odstranitvijo nesingularnosti FR smo točkovne izvire nadomestili s porazdeljenimi izviri po krogih. Vrednosti porazdeljenih izvirov so izračunane v analitični obliki za Dirichletove robne pogoje. Pri Neumann-ovih robnih pogojih so vrednosti odvodov FR izračunane posredno pri upoštevanju rešitev preprostih polj premika. Obravnavani so sistemi z dvema, štirimi in devetimi telesi v stiku. Nova metoda je verificirana na podlagi primerjave s klasično MFR. Uporabljena bo v deformacijskem modelu mikrostukture, povezanim z makroskopskim termomehanskim sistemom za simulacijo kontinuirnega ulivanja, vročega valjanja in toplotne obdelave.
Ključne besede: izotropna elastičnost, ravninska deformacja, Navierove enačbe, robni pogoji premika in vlečenja, nesingularna metoda fundamentalnih rešitev, Kelvinova fundamentalna rešitev
1 INTRODUCTION
The physical modelling of metallurgical processes1
consists of modelling the relations between the process parameters and the macroscopic velocity, temperature,
concentration, and stress fields, and the relations bet-
ween these fields and the evolution of the microstructure. In such multi-scale modelling, the transport phenomena and solid mechanics on the level of microstructure play an important role and have to be properly computationally modelled.23 We have, in the recent years, developed a completely new generation of meshless methods, based on local collocation with radial basis functions, for
solving these models on different scales. The main advantage of the method is in their similar structure in 2 and 3D, no need for polygonisation, ease of coding, high accuracy, robustness and flexibility. The method has already been developed for modelling very complex
phenomena such as macro-segregation on the macroscopic level4 as well as dendritic growth on the micro-level.5 An extension of the meshless method, based on a collocation with a fundamental solution (Method of Fundamental Solutions (MFS)), for the simulation of the deformation of the multiple grains in an ideal mechanical contact is presented in the present paper. An extension of the represented method with anisotropic and plastic deformation capabilities will be used in our future ther-mo-mechanical calculations for the microstructure evolution in metallurgical processes.6
The main idea of MFS consists of approximating the solution of the partial differential equation by a linear combination of fundamental solutions, defined in source points. The expansion coefficients are calculated by collocation or least-squares fit of the boundary conditions. The fundamental solution is usually singular in the
source points and this is the reason why the source points are located outside the domain in the MFS. Then, the original problem is reduced to determining the unknown coefficients of the fundamental solutions and the coordinates of the source points by requiring the approximation to satisfy the boundary conditions and hence solving a non-linear problem. If the source points are a priori fixed, then the coefficients of the MFS approximation are determined by solving a linear problem. The MFS has become very popular in recent years because of its simplicity. Clearly, it is applicable when the fundamental solution of the partial differential operator of the governing equation (or of the system of governing equations) of the problem under consideration is known. The MFS has been widely used7 for the solution of problems in linear elasticity.
In the traditional MFS, the fictitious boundary, positioned outside the problem domain, is required to place the source points. This is very impractical or even impossible, particularly when solving muti-body problems. In recent years, various efforts have been made, aiming to remove this barrier in the MFS, so that the source points can be placed directly on the real boundary.8-12 In the present paper, we use a non-singular MFS based on8 to deal with the 2D multi-body isotropic elasticity problems. The application of a non-singular method of fundamental solutions (NMFS) in two-dimensional isotropic linear elasticity has been originally developed in.13 We extend these developments to multi-body problems in the present paper. We respectively use area-distributed sources covering the source points to replace the concentrated point sources. This NMFS approach also does not require information about the neighbouring points for each source point, thus it is a truly mesh-free boundary method. The derivatives of the fundamental solution in the distributed source points are calculated by adopting the methodology in9 from the Laplace to Kelvin fundamental solution.
The rest of the paper is structured as follows. The solution procedure is given for MFS and NMFS. Numeral results with 1, 4 and 9 bodies in contact are given, followed by conclusions and further research.
2 GOVERNING EQUATIONS OF ELASTICITY FOR THE MULTI-BODY PROBLEM
We consider a two-dimensional domain fi with the boundary T, divided into M sub-domains fi = fii U Qn U ... U fiM with boundaries T = (T U Tn U ... U Tm) -Ti-ii - ... - Ti-m - ... - T(M-I)-M as shown in Figure 1. Each of the domains is occupied by an isotropic, ideally elastic material with different material properties, in general. Let us introduce a two-dimensional Cartesian coordinate system with orthonormal base vectors ix and iy, and coordinates px and py of point P with the position vector p = Pxix + Pyiy. The solid is governed by Navier's equations for plane strain problems, which are the conditions for
equilibrium, expressed with the displacement u. The following governing equations are valid in the subdomain Qm,m = I, II, ..., M, p £ Qm.
2(1-vm) d2ux (p) d2ux (p) 1 d2Uy (p)
1-2v.
dp
dp
1-2v.
2(l-vm ) d 2 Uy (p) d 2 Uy (p)
■+ ^ , +
l
dpx dp y d2Ux(p)
1-2v.
dp
dp
1-2v.
= 0
= 0
(l)
dPx dPy
where Vm represents the Poisson ratio in the subdomain Qm. The boundary T is divided into two not necessarily connected parts T = TD + Tr. On the part TD the displacement (Dirichlet) boundary conditions are given, and on the part Tr the traction (Neumann) boundary conditions are given:
X i« a (p)+X h « ai (p)+. . m u çM (p) = «ç (p) Ç = x, y, p e T D
X1t Çi (p) + X il t Çii (p)+...+Z M t ÇM (p) = tç (p)
Ç = x, y, p e T T
where:
[1, p e T
Xm=
(2)
0, p Ž r „
(3)
On the interface between different regions, displacement continuity and traction equilibrium conditions have been assumed:
«Çm(p)-«Çk(p) = 0 Ç = X, y, p e T m ^ T t
tçm(p)+tçk(p) = 0 ç = X, y, p e T m n T t	(4)
m, k = i,ii,..., M
The strains Ç, ^ = x, y are related to the displacement gradients by:
dUç dU^ dp( dpç
(5)
		y		
rD				
—	£2t		OtI	
		0		X
n	P„		"M	rJ
				
Figure 1: A scheme of the multi-region problem. Each of the subdomains can have different elastic properties.
Slika 1: Shema problema z več območji. Vsako podobmocje ima lahko različne elastične lastnosti.
The stress components o^; £ = x, y are for the planestrain cases related to the strains through Hooke's law:
=X m <5 S (£ XX + £ yy ) + 2[ m £
m° K
(6)
where ^m = Em /2(1 + vm) is the shear modulus of elasticity, Em is a constant, known as the modulus of elasticity, or Young's modulus, Xm = 2 Vmßm / (1 - 2vm) is the Lamé constant, and ôçs is the Kronecker delta:
[1, £ = S
ô£s = k £ * S	(7)
3 SOLUTION PROCEDURE
The fields on each of the sub-domains are represented by collocation with fundamental solutions in the boundary points. The collocation needs to satisfy the boundary conditions between different regions and outer boundaries. In a numerical implementation of MFS and NMFS, we assume that one boundary collocation point belongs to only two regions at once. In order to keep the formulation simple we do not put the discretisation points on the corners that might belong to three or more regions at once. Explicit expressions for Kelvin's fundamental solution of elastostatics, used in the collocation, are given14 in a two-dimensional plane-strain situation by:
p, s) =
= x, y
{(3- *)„ri k+lää*
(8)
where the material properties depend on the position in a subdomain. Uçs (p,s) represents the displacement in the £ direction at point p due to a unit point force acting in the S direction at point s. r =
is the distance between the collocation point p and the source point s.
It can be shown that the following Ux and Uy satisfy the governing Eq. (1):
N	N
(p ) =X (p, p n )« n +X Vy ( p, p „ ) ß n
n= 1 N
n= 1 N
Uy (P ) =X Uyx (P, P B )ff B +X Uy (p, p „ )fi n	(9)
n=1	n=1
N
P i! A(p n, R)
n=1
where an and fin represent arbitrary constants and an = anm, fin = finm, N = Nm when p G Qm. Nm is the number of p G Tm, m = I, II, ..., M. A(pn,R) (Figure 2) represents a circle with radius R, centred around pn ■ pn = s represents points on the physical boundary. p and pn belong to the same sub-domain. The quantity Ug (p,pn) is singular when p = pn. We use the following de-sin-gularization technique, proposed by6, for evaluating the singular values:
Ug (P, P n )	P * P n
(P" P ' \R J U„ (P, Pn ^ P = Pn
(10)
A(s ,R)
The tractions can be expressed as:
N	N
tx ( P) = X Txx fo P n )« n + X Txy ( P, P n )ß n
n=1	n=1
NN
t y ( P) = X Tyx fo P n ) a n + X Tyy ( P, P n ) ß n
n=1	n=1
(11)
where:
Txx( P, Pn ) =
2ß(1 - V) dUxx(P, Pn ) + 2ßV_ dUyx(P, Pn )
1 - 2v dpx 1 - 2v dp
ß U(P, Pn ) + ß dUyx(P, Pn )
Txy( P, Pn ) =
dP y	dPx
2 fl(1 - V) dUxy(P, Pn ) 2[IV dUyy(P, Pn )
1 - 2v dpx 1 - 2v dp
dUxy(P, Pn ) dUyy(P, Pn ) [---+
dPy
Tyx ( P, Pn ) =
[
dPx
d U yx ( P, Pn ) jUx ( P, Pn )
[
dPx	dP y
2[(1 - V) dUyx(P, Pn ) + 2[V_ dUxx( P, Pn ) 1 - 2v dP	1 - 2v dPx
n„x +
Tyy ( ^ Pn ) =
dU yy ( P, Pn ) dUxy ( P, Pn )
dPx	dP y
2[(1 - V) dUyy(P, Pn ) 2[V dUxy(P, Pn )
n+
(12)
1 - 2v dpy 1 - 2v dpx tç = tçm, Tçç = Tgm, nnÇ = HnÇm, Ç, £ = x, y when p G Qm.
The coefficients an and fin are calculated from a system of (2Ni + 2Nii + ... + 2nm) x + (2ni + 2Nn + ... + 2nm) algebraic equations, obtained by collocating the boundary conditions:
Ax = b	(13)
Figure 2: Distributed source on a disk A(pn) Slika 2: Porazdeljeni izviri na disku A(pn)


+
nx
where A is composed of (p, pn) and Tg (p, pn) x is
composed of an and fin, and b is composed of uq, t- and 0. The explicit form of the elements of the algebraic equation system (13) can be found in.11
The diagonal terms Tg (pl, p l ) Q, £ = x, y, l =1, ...,
Ni + ... + Nm, in Eq. (13) are determined indirectly for the collocation points on Tr. For this purpose, the method proposed in7 for potential problems is applied to determine the diagonal coefficients of Eq. (13). In the approach, we first assume two simple solutions. The first simple solution is ux (p) = px + cx, uy (p) = 0 everywhere. The second simple solution is ux (p) = 0, uy (p) = py + c everywhere. We solve them for the corresponding a fiV and a(2), fi(2) using only the displacement boundary condition. From these two solutions, we can also know:
M\p) 5m<1)( p)

=1

( p)

( p)

(p)

dPy ( p)

dPx ( p)
dPy
dPy
dPy
= 0
dPy
day2"
= 0
(14)
( p)
dPy
= 1
By substituting a n°, fi ^ and a (n2\ fi (n2) and Eq. (9) into Eq. (14), we can obtain the diagonal terms T (pl,pi) Q, £ = x, y, l =1, ..., Ni + ... + Nm. The
constant c should be selected in such a way that all the points in the upper cases do not move the same distance. So that the denominators in the upper derivations are not zero.
By knowing all the elements Aij and bk of the system (13), we can determine all the values of an and fin. Subsequently, we can calculate the displacement for all the domain points using Eq. (9).
Figure 4: The deformation calculated with MFS and NMFS for a four-domain case EI = En = ... = E\V = 1N/m2, vi = vu = ... = viv = 0.3 and N = 240 (•: collocation points, o: source points, x: MFS solution, A: NMFS solution). The position of the source points in MFS is on squares around each of the square physical domains. Slika 4: Deformacja, izracunanana z MFR in NMFR, za primer s štirimi območji EI = EII = ... = EIV = 1N/m2, vI = vII = ... = vIV = 0,3 in N = 240 (•: kolokacijska točka, o: izvirna točka, x: MFR-rešitev, A: NMFR-rešitev). Pozicija izvirnih točk je na kvadratih okoli vsakega od štirih kvadratnih domen.
We consider a square with side a = 3m centred around px = 0m, py = 0m for testing the performance of the method. We distinguish three sub-examples. In the first one, the whole square is occupied by one material, with the material properties E = 1N/m2, v = 0.3, in the second one, the square is split into four parts with the same material properties as in the first example EI = EII = ... = Eiv = 1N/m2 Vi = Vii = ... = Viv = 0.3, and in the third one, the square is split into nine parts with the same
4 NUMERICAL EXAMPLES
Figure 3: The deformation calculated with MFS and NMFS for a one-domain case with E = 1N/m2, v = 0.3 and N =120 (•: collocation points, o: source points, x: MFS solution, A: NMFS solution) Slika 3: Deformacija, izračunana z MFR in NMFR, za primer z enim območjem z E = 1N/m2, v = 0,3 in N = 120 (•: kolokacijska točka, o: izvirna točka, x: MFR-rešitev, A: NMFR-rešitev)
Figure 5: The deformation calculated with MFS and NMFS for a nine-domain case EI = EII = ... = EIV = 1N/m2, vI = vII = ... = vIV = 0.3 and N = 360 (•: collocation points, o: source points, x: MFS solution, A: NMFS solution). The position of the source points in MFS is on squares around each of the square physical domains. Slika 5: Deformacija, izračunana z MFR in NMFR, za primer z devetimi območji EI = EII = ... = EIV = 1N/m2, vI = vII = ... = vIV = 0,3 in N = 360 (•: kolokacijske točke, o: izvirne točke, x: MFR-rešitev, A: NMFR-rešitev). Pozicija izvirnih točk je na kvadratih okoli vsakega od štirih kvadratnih domen.
material properties as in the first example Ei = En = ... = Eix = 1N/m2 vi = vii = ... = Vix = 0.3. We considered the solution of the Navier equations in this square subject to the boundary conditions ux = 0m, uy =-01mon the points of the north side with py = 1.5m; ux = 0m, uy = 01m on the south side with py = -1.5m and on the east px = 1.5m and west px = -1.5m sides tx = 0N/m2, ty = 0N/m2 is set. A plot of the deformation, calculated with the defined three sub-examples is shown in Figures 3, 4 and 5, respectively. The following parameters were used R = d/5, Ri = di/5,..., Rix = dix/5, where d, dm, m = I, II,..., IX are the smallest distances between two nodes on the boundary, Rm is the radius of the circle centred the point p„ G Tm, Cx = Cy = Cxi = Cyl = L = Cxix = Cyix = 4. The distance of the fictitious boundary from the true boundary for the MFS is set to Rm = 5d, Rmi = 5di, ..., Rmix = 5dix.
Figures 3, 4 and 5 show good agreement between the solution for a one-domain region and a solution recalculated with the four and nine regions in ideal mechanical contact and with the same material properties. The maximum absolute difference in displacements between the values in Figures 3 and 4 at the outer boundary are Aux = 0.0417m, Auy = 0.0886m, and between Figures 3 and 5 Aux = 0.0017m, Auy = 0.0012m, respectively.
5 CONCLUSION
A new, non-singular method of fundamental solutions13 is extended in the present paper to solve multi-dimensional linear elasticity problems. In this approach, the singular values of the fundamental solution are integrated over small circular discs, so that the coefficients in the system of equations can be evaluated analytically and consistently, leading to an extremely simple computer implementation of this method. The method essentially gives the same results as the classic MFS. It has the advantage that the artificial boundary is not present; however, at the expense of solving three times the systems of algebraic equations in comparison with only one solution in MFS. The main advantage of the method is that the discretisation is performed only on the boundary of the domain and no polygonisation is needed, like in the finite-element method. The NMFS method, presented in this paper, can be adapted or extended to handle many related problems, such as three-dimensional elasticity, anisotropic elasticity, and multi-body problems, which all represent directions for
our further investigations. The advantage of not having to generate the artificial boundary is particularly welcome in these types of problems. The method will be used in the future for the calculation of multigrain deformation problems in steel and aluminium alloys, with realistic grain shapes, obtained from the microscope images. The developed method is believed to represent a most simple state-of-the-art way to numerically cope with these types of problems.
Acknowledgement
This paper forms a part of the project L2-3651 Simulation and Optimisation of Casting, Rolling and Heat Treatment Processes for Competitive Production of Topmost Steels, supported by the Slovenian Research Agency (ARRS) and Store Steel company. The Centre of Excellence for Biosensors, Instrumentation and Process Control (COBIK) is an operation financed by the European Union, European Regional Development Fund and the Republic of Slovenia.
6 REFERENCES
'B. Sarler, R. Vertnik, S. Saletic, G. Manojlovic, J. Cesar, Berg-Huettenmaenn. Monatsh., 150 (2005), 300-306
2	T. Mura, Micromechanics of Defects in Solids, 2ed, Martinus Nijhoff Publishers, Netherlands 1987
3	M. Braccini, M. Dupeux, Mechanics of Solid Interfaces, Wiley, New York 2012
4G. Kosec, M. Zaloznik, B. Sarler, H. Ccombeau, CMC: Computers,
Materials & Continua, 22 (2011), 169-195 5 A. Z. Lorbiecka, B. Sarler, CMC: Computers, Materials & Continua,
18 (2010), 69-103 6U. Hanoglu, S. Islam, B. Sarler, Mater. Tehnol., 45 (2011) 6, 545-547
7Y. S. Smyrlis, Mathematics of Comptation, 78 (2009), 1399-1434
8	Y. J. Liu, Engineering Analysis with Boundary Elements, 34 (2010), 914-919
9	B. Sarler, Engineering Analysis with Boundary Elements, 33 (2009), 1374-1382
10	D. L. Young, K. H. Chen, J. T. Chen, J. H. Kao, CMES: Computer Modeling in Engineering & Sciences, 19 (2007), 197-222
11	D. L. Young, K. H. Chen, C. W. Lee, Journal of Computational Physics, 209 (2005), 290-321
12	W. Chen, F. Z. Wang, Engineering Analysis with Boundary Elements, 34 (2010), 530-532
13	Q. G. Liu, B. Sarler, CMES: Computer Modeling in Engineering and Sciences, 91 (2013), 235-266
14	D. E. Beskos, Boundary Element Methods in Mechanics, Elsevier, Amsterdam 1987, 33