- Email: [email protected]
Planar Slip Condition For Mesh Morphing Using Radial Basis Functions
Available online at www.sciencedirect.com
Available online at www.sciencedirect.com Available online at www.sciencedirect.com
ScienceDirect Procedia Engineering 20300 (2017) 349–361 Procedia Engineering (2017) 000–000 Procedia Engineering 00 (2017) 000–000
www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia
26th 26th International International Meshing Meshing Roundtable, Roundtable, IMR26, IMR26, 18-21 18-21 September September 2017, 2017, Barcelona, Barcelona, Spain Spain
Planar Planar Slip Slip Condition Condition For For Mesh Mesh Morphing Morphing Using Using Radial Radial Basis Basis Functions Functions a,∗ a,b a a St´ Rendu St´eephane phane Aubert Auberta,∗,, Franck Franck Mastrippolito Mastrippolitoa,b ,, Quentin Quentin Rendua ,, Martin Martin Buisson Buissona ,, Fr´ Fr´eed´ d´eeric ric b Ducros b Ducros a Univ. a Univ.
Lyon, Ecole Centrale de Lyon, Laboratoire de M´ecanique des Fluides et d’Acoustique UMR5509, Lyon, F-69134, France Lyon, Ecole Centrale de Lyon, Laboratoire de M´ecanique des´ Fluides et d’Acoustique UMR5509, Lyon, F-69134, France b Univ. Alpes, CEA, LITEN, DTBH, Laboratoire Echangeur et R´eacteurs, F-38000 Grenoble, France b Univ. Grenoble ´ Grenoble Alpes, CEA, LITEN, DTBH, Laboratoire Echangeur et R´eacteurs, F-38000 Grenoble, France
Abstract Abstract An original approach to morph meshes including slip planes is proposed. Extending the classical interpolation scheme based An original approach to morph meshes including slip planes is proposed. Extending the classical interpolation scheme based on Radial Basis Functions (RBF), it preserves the linear property, linking the deformation response to the known displacements, on Radial Basis Functions (RBF), it preserves the linear property, linking the deformation response to the known displacements, underlying the RBF formalism. No post-correction stage, nor extra not-moving control points, are required to enforce in plane underlying the RBF formalism. No post-correction stage, nor extra not-moving control points, are required to enforce in plane displacements. Accuracy and robustness are evaluated for a 2D configuration where a shock-wave impinges on a cross-flow displacements. Accuracy and robustness are evaluated for a 2D configuration where a shock-wave impinges on a cross-flow flexible panel. Analysis is based on deformation quality metrics derived from the displacement field gradient, keeping the proposed flexible panel. Analysis is based on deformation quality metrics derived from the displacement field gradient, keeping the proposed approach free from connectivity information. 3D applications are briefly illustrated by the tip of an aero engine fan subjected to approach free from connectivity information. 3D applications are briefly illustrated by the tip of an aero engine fan subjected to conjugated parameterized variations of tip clearance and stagger angle. conjugated parameterized variations of tip clearance and stagger angle. c 2017 The Authors. Published by Elsevier Ltd. c 2017 2017 The The Authors. Authors. Published Published by by Elsevier Elsevier Ltd. Ltd. © Peer-review under responsibility of the scientific committee of the 26th International Meshing Roundtable. Peer-review under responsibility of of the the scientific scientific committee committee of of the the 26th 26th International International Meshing Meshing Roundtable. Roundtable. Keywords: radial basis function mesh deformation, slip condition, CFD Keywords: radial basis function mesh deformation, slip condition, CFD
1. Introduction 1. Introduction In many CFD-applications, the mesh has to be moved after or during the simulation process. This feature is In many CFD-applications, the mesh has to be moved after or during the simulation process. This feature is typically implemented in optimization [1,2] and Fluid-Structure Interaction [3–5] solvers. A main issue is that the typically implemented in optimization [1,2] and Fluid-Structure Interaction [3–5] solvers. A main issue is that the mesh displacement is highly strained by grid-quality constrains (skewness, aspect ratio, volume positivity), boundary mesh displacement is highly strained by grid-quality constrains (skewness, aspect ratio, volume positivity), boundary conditions (wall, periodic, symmetric) or prescribed surface deformations. conditions (wall, periodic, symmetric) or prescribed surface deformations. Morphing a mesh using springs analogy [6,7], pseudo-solid approach [8] or free form deformation [9,10] is a well Morphing a mesh using springs analogy [6,7], pseudo-solid approach [8] or free form deformation [9,10] is a well studied subject. But to propagate displacements known only at clouds of control points to the whole volume mesh, studied subject. But to propagate displacements known only at clouds of control points to the whole volume mesh, the Radial Basis Function (RBF) interpolation method is more and more popular, after its introduction for external the Radial Basis Function (RBF) interpolation method is more and more popular, after its introduction for external aeroelasticity applications [1,4]. Based only on initial position and displacement of these control points, it defines aeroelasticity applications [1,4]. Based only on initial position and displacement of these control points, it defines by solving a linear system a continuous displacement field that requires no connectivity information. It is thus very by solving a linear system a continuous displacement field that requires no connectivity information. It is thus very flexible to deal with structured or unstructured grids, meshing fluid or solid domains. Besides, its memory requirement flexible to deal with structured or unstructured grids, meshing fluid or solid domains. Besides, its memory requirement ∗ ∗
Corresponding author. Tel.: +3-347-218-6741 Corresponding author. Tel.: +3-347-218-6741 E-mail address: [email protected], [email protected] E-mail address: [email protected], [email protected]
c 2017 The Authors. Published by Elsevier Ltd. 1877-7058 c 2017 The Authors. Published by Elsevier Ltd. 1877-7058 Peer-review©under responsibility of the scientific committee of the 26th International Meshing Roundtable. 1877-7058 2017responsibility The Authors. Published Elsevier of Ltd. Peer-review under of the scientificbycommittee the 26th International Meshing Roundtable. Peer-review under responsibility of the scientific committee of the 26th International Meshing Roundtable. 10.1016/j.proeng.2017.09.819
Aubert et al. /Engineering Procedia Engineering 203 (2017) 349–361 S. St´ephane Aubert et al. / Procedia 00 (2017) 000–000
350 2
Nomenclature b d ei n Nc Nm Ns r0 s smax smi S t x xci xmi x si γi Γ
Second unit tangent vector, 3i=1 bi ei (∈ R3 ). Distance measuring function (R3 → R+ ). ith basis vector (∈ R3 ). Unit normal vector, 3i=1 ni ei (∈ R3 ). Number of control points (∈ N). Number of moving points (∈ N). Number of sliding points (∈ N). Support radius (∈ R). Displacement function (R3 → R3 ). Largest displacement magnitude (∈ R). Displacement vector of ith moving point (∈ R3 ). Vector (∈ RNm ). First unit tangent vector, 3i=1 ti ei (∈ R3 ). Position vector (∈ R3 ). Position vector of ith control point (∈ R3 ). Position vector of ith moving point (∈ R3 ). Position vector of ith sliding point (∈ R3 ). Weight vector of ith control point (∈ R3 ). Vector (∈ RNc ).
δ ε φ φ Φm Φs Φt 0s 0 s×m 0m×c I s×s u · v u ∧ v u ⊗ v u ∇u ⇒
⇒
∇ u AT ∗
Tip gap height (∈ R). Strain tensor (∈ R3×3 ). Radial basis function (R+ → R). Derivative of φ (R+ → R). Matrix (∈ RNm ×Nc ). Matrix (∈ RNs ×Nc ). Matrix (∈ RNs ×Nc ). Null vector (∈ RNs ). Null matrix (∈ RNs ×Nm ). Null matrix (∈ RNm ×Nc ). Identity matrix (∈ RNs ×Ns ). Dot product, 3i=1 ui vi ((R3 , R3 ) → R). Cross product, [u j vk − uk v j ]i ((R3 , R3 ) → R3 ). Outer product, [ui v j ]i, j ((R3 , R3 ) → R3×3 ). √ Euclidean norm, u · u (R3 → R). Gradient operator, [∂u/∂xi ]i (R → R3 ).
Gradient operator, [∂ui/∂x j ]i, j (R3 → R3×3 ). Transpose operator, [a ji ]i, j (R3×3 → R3×3 ). After deformation state. Vector components written in ( e x , ey , ez ) space.
is mesh size independent [11]. However standard RBF formulation is isotropic. Mandatory geometrical properties, like planar shapes to define rotor/stator interfaces, are likely to be spoiled. Furthermore, using functions with compact support to localize mesh deformations, action radii should be at least three times larger than the control points displacement to smooth these over a long enough distance. It is thus tricky to setup RBF data for confined space configurations, like blade tip clearance and under-platform cavities. To circumvent such drawbacks, some classical approaches rely on adding enclosing boxes of not-moving control points [2] or implementing predictor/corrector strategies based on directional damping functions [1,12]. All try and preserve selected space regions by cancelling some undesired movements after the RBF setup. Unfortunately, there is no (retro)action of such corrections on the interpolation coefficients; beside these require many user defined extra data, which are case dependent. Thereafter, the proposed approach is firstly presented in section 2. It is based on the formulation of a slip condition, along one or more planes, for subsets of control points instead of specifying their whole displacement. Prescribed moving and sliding control points are taken into account together in the linear system defining the interpolation coefficients. So, no post-correction stage is required, nor extra not-moving control points. Deformation quality metrics inspired by continuum mechanics are also introduced. Based on the displacement field gradient, these keep the proposed approach free from connectivity information. Then, accuracy and robustness are analyzed in section 3.2 through a 2D configuration where a shock-wave impinges on a cross-flow flexible panel. 3D applications are briefly illustrated in section 3.3 by the tip of an aero engine fan subjected to conjugated parameterized variations of tip clearance and stagger angle.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
351 3
2. Formulation 2.1. Displacement field interpolation RBF-based mesh deformation is a node-by-node scheme free from mesh connectivity. It propagates by interpolation the known motion of a cloud of control points to the mesh nodes. The interpolated function is the displacement s of a node, defined as the weighted sum of radial basis functions by: s( x) =
Nc j=1
φ(d( x − xc j )) γ j
(1)
where x is the initial (i.e. before deformation) node position, Nc is the number of control points, φ is the chosen radial basis function, d is the chosen function measuring distance between two positions, xc j and γ j are respectively the initial position and weight of the jth control point. It is to notice that only initial positions appear explicitly in Eqn. (1). The prescribed displacement of control points and possibly other constraints are embedded in the weights { γi }1≤i≤Nc . When these are known, s is fully defined. Thereafter, the set of control points is split in two. Different conditions are associated to each subset to close the problem defining the weights. Let Nm ≤ Nc be the number of moving control points, i.e. with fully prescribed motion. Fixed, not moving, control points are included in this subset with zeroed displacement. Let N s = Nc − Nm be the number of sliding control points, i.e. with motion constrained to remain in a plane. The set of control points is assumed to be ordered, starting with the moving ones with index ranging from 1 to Nm , followed by the sliding ones with index ranging from Nm + 1 to Nm + N s = Nc . 2.2. Interpolation condition On one hand, the displacement s must satisfy the interpolation condition for each moving control point: ∀i ∈ [1, Nm ], ⇔ ⇔
∀k∈[1,3]
s( xmi ) = smi Nc φ(d( xmi − xc j )) γ j = smi
(2) (3)
j=1
Nc j=1
φ(d( xmi − xc j )) γ j · ek = smi · ek
(4)
where smi is the known displacement of the ith moving point initially at position xmi , and ( e1 , e2 , e3 ) is a basis of the 3D Euclidean space. Equation (4) is rewritten as: (5) ∀k ∈ [1, 3], Φm Γk = S k where matrix Φm and vectors Γk and S k are defined as: Φm = φ(d( xmi − xc j )) 1≤i≤Nm ,1≤ j≤Nc = φ(d( xci − xc j )) 1≤i≤Nm ,1≤ j≤Nc
∀k ∈ [1, 3],
Γk = γi · ek 1≤i≤Nc S k = smi · ek 1≤i≤Nm
(6) (7) (8) (9)
It is to notice that if Nm = Nc , Φm is a square matrix and the components of { γi }1≤i≤Nc are simply solution of the three linear systems defined by Eqn. (5). It is the classical RBF formulation where the Γk vectors are independent.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
352 4
2.3. Planar slip condition On the other hand, the displacement s must satisfy three scalar conditions, associated to the normal and tangential directions of the slip planes. For the sake of simplicity thereafter, with no loss of generality, only one such plane is considered. Let n be its unit normal vector and t a unit tangent vector, such that t · n = 0. Another unit tangent vector b is simply defined as b = n ∧ t. 2.3.1. Zeroed normal displacement The first condition on s states that its normal component s · n is zero for each sliding control point: ∀i ∈ [1, N s ], s( x si ) · n = 0 3 Nc e =0 n φ(d( x si − xc j )) γ j · ⇔ k k j=1 k=1 Nc 3 =0 nk φ(d( x − x )) γ · e ⇔ si cj j k j=1 k=1
(10) (11)
(12)
where x si is the initial position of the ith sliding point and (n1 , n2 , n3 ) are the components of n in ( e1 , e2 , e3 ) basis. Equation (12) is rewritten as: 3 nk Φ s Γk = 0 s (13) k=1
Ns
s
where 0 s is the null vector in R and matrix Φ is defined as: Φ s = φ(d( x si − xc j )) 1≤i≤N s ,1≤ j≤Nc = φ(d( xc(Nm +i) − xc j )) 1≤i≤N s ,1≤ j≤Nc
(14) (15)
2.3.2. Zeroed tangential contributions The second condition on s states that no sliding control point is contributing to the displacement function in the tangential direction t: ∀i ∈ [Nm + 1, Nm + N s ],
γi · t = 0 3 e =0 t ⇔ γi · k k
(16) (17)
k=1
⇔
3 k=1
tk γi · ek = 0
(18)
where (t1 , t2 , t3 ) are the components of t in ( e1 , e2 , e3 ) basis. Equation (18) is rewritten as: 3 tk Φt Γk = 0 s
(19)
k=1
where matrix Φt is defined as:
Φt = 0 s×m I s×s ∈ RNs ×Nc
(20)
with 0 s×m the null matrix in RNs ×Nm and I s×s the identity matrix in RNs ×Ns . The third and last condition on s states that no sliding control point is contributing to the displacement function in the other tangential direction b, i.e. substituting b to t in Eqn. (19): 3 k=1
bk Φt Γk = 0 s
(21)
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
353 5
where (b1 , b2 , b3 ) are the components of b in ( e1 , e2 , e3 ) basis. 2.4. Weights calculation Therefore, the components of { γi }1≤i≤Nc are solution of the (3Nc ) × (3Nc ) linear system assembled from Eqn. (5), Eqn. (13), Eqn. (19) and Eqn. (21): m S 1 Φ 0m×c 0m×c 0m×c Φm 0m×c S 2 Γ 0m×c 0m×c Φm 1 S 3 (22) Γ2 = n1 Φ s n2 Φ s n3 Φ s Γ 0 s 0 t Φt t Φt t Φt 3 2 3 s 1 t t t 0s b1 Φ b2 Φ b3 Φ with 0m×c the null matrix in RNm ×Nc . When multiple slip planes are considered, Eqn. (13), Eqn. (19) and Eqn. (21) must be replicated for each one of them with associated ( n, t, b) triplet.
It is to notice that the Γk vectors are coupled in Eqn. (22) by the matrix last three lines in the general case, compared to the classical RBF formulation. However if a single slip plane is considered, an uncoupled formulation is recovered in the basis ( e1 , e2 , e3 ) = ( n, t, b). The components of { γi }1≤i≤Nc are then again solution of three Nc × Nc linear systems:
∀k ∈ [2, 3],
Φm S Γ = 1 Φs 1 0s m Φ S Γ = k 0s 0 s×m I s×s k
(23) (24)
On one hand, Eqn. (23) can be interpreted as the classical RBF formulation projected on n, with the sliding points treated as fixed control points. On the other hand, Eqn. (24) can be interpreted as the classical RBF formulation projected on t or b and restricted to the first Nm control points, i.e. to the moving ones. 2.5. Algorithms Several operations are required to morph a mesh. Algorithm 1 summarizes the steps of the core one yielding the position { xi∗ }1≤i≤N of the deformed mesh nodes knowing the initial one { xi }1≤i≤N , the control points initial position { xci }1≤i≤Nc and the moving control points displacement { smi }1≤i≤Nm . After assembling and solving the linear system defined by Eqn. (22), the nodes displacement is interpolated with Eqn. (1) (line 5), then it is simply added to the nodes initial position (line 7). The number and position of the control points are key factors for the efficiency and accuracy of Alg. 1. On one hand, the interpolation/update step CPU time is proportional to N(Nc + 1) (lines 5 and 7), and computing the weights is of Nc3 complexity (line 3). On the other hand, too close control points can jeopardize the condition number of matrices Φm and Φ s , because of too small values of distance function d. Therefore a greedy procedure [13,14] is applied to build a quasi-optimal subset of the initially provided control points. It is presented in Alg. 2. The subset is by default initialized with the control point exhibiting the largest displacement smax (line 1). The idea is to apply the displacement field computed using only the control points of the current subset, to move all the provided control points, using Alg. 1 (line 3). The displacement of the subset control points is exactly interpolated. So, one control point not already chosen is exhibiting the largest error, defined from Eqn. (2) for moving points (line 5) and from Eqn. (10) for sliding points (line 6). If its error is greater than some tolerance, usually defined as a fraction of smax , it is added to the subset (line 9) and the whole process is repeated.
6 354
S. St´ephane Aubert et al. / Procedia 00 (2017) 000–000 Aubert et al. /Engineering Procedia Engineering 203 (2017) 349–361
Algorithm 1: Nodes displacement
1 2 3 4 5 6 7
Input: Nodes initial position { xi }1≤i≤N Input: Control points initial position { xci }1≤i≤Nc Input: Moving points displacement { smi }1≤i≤Nm Output: Nodes final position { xi∗ }1≤i≤N Compute the weights { γi }1≤i≤Nc : Assemble the matrix and RHS of Eqn. (22) Solve Eqn. (22) to get (Γ1 , Γ2 , Γ3 ) Interpolate the nodes displacement: c ∀i ∈ [1, N], s( xi ) = Nj=1 φ(d( xi − xc j )) γ j Update the nodes position: ∀i ∈ [1, N], xi∗ = xi + s( xi )
Algorithm 2: Control points optimization
1 2 3 4 5 6 7 8 9 10 11
Input: Control points initial position { xci }1≤i≤Nc Input: Moving points displacement { smi }1≤i≤Nm Output: Control points subset Initialize the control points subset repeat Using Alg. 1, move all the control points, input as nodes while the subset is input as control points Compute errors: ∗ ∀i ∈ [1, Nm ], eval ( xmi − xmi ) − smi ∗ ∀i ∈ [1, N s ], eval |( x si − x si ) · n| Find the point with the largest error if error > tolerance then Add this point to the control points subset end until error ≤ tolerance;
The equations involved in Alg. 1 and Alg. 2 are written in ( e1 , e2 , e3 ) vector space. In particular, normal vectors n are solely defined in this space. It can be convenient to choose it as different from the standard Cartesian coordinate system ( e x , ey , ez ). For example, a conical surface in ( e x , ey , ez ) space is planar in the associated cylindrical coordinate system. However, special care must be taken if switching between ( e x , ey , ez ) and ( e1 , e2 , e3 ) spaces is not a linear process. Algorithm 3 presents the steps to morph a mesh in such a case. Before using Alg. 2 then Alg. 1 (lines 9 and 10), all x → x) (line 3). Because the position vectors known in ( e x , ey , ez ) basis must be rewritten in ( e1 , e2 , e3 ) basis (e.g. displacement is the difference between two positions, it can not be simply converted if the relation between the two spaces is not linear. Therefore, the after deformation position of the moving points is firstly calculated (line 2). Then it is rewritten in the new basis (line 6). Finally, the displacement in ( e1 , e2 , e3 ) space is recovered by subtracting initial positions from final ones (line 8). The last step (line 12) brings back the moved nodes to the original ( e x , ey , ez ) space. Algorithm 3: Overall procedure
1 2 3 4
xi }1≤i≤N Input: Nodes initial position { xci }1≤i≤Nc Input: Control points initial position { Input: Moving points displacement { s mi }1≤i≤Nm xi∗ }1≤i≤N Output: Nodes final position { Compute the moving points final position in ( e x , ey , ez ): s ∀i ∈ [1, Nm ], x ∗ = x mi + mi mi
Switch from ( e x , ey , ez ) to ( e1 , e2 , e3 ) space: xi ∀i ∈ [1, N], eval xi from
5 6 7 8 9 10 11 12
xci eval xci from ∗ ∗ ∀i ∈ [1, Nm ], eval xmi from x mi Compute the moving points displacement in ( e1 , e2 , e3 ): ∗ ∀i ∈ [1, Nm ], smi = xmi − xmi Select an optimal subset of control points using Alg. 2 Move the nodes using Alg. 1 and this optimal subset Switch back from ( e1 , e2 , e3 ) to ( e x , ey , ez ) space: x ∗ from x ∗ ∀i ∈ [1, N], eval ∀i ∈ [1, Nc ],
i
i
2.6. Deformation quality metrics The quality of the displacement field s is classically evaluated through the quality of the morphed mesh x ∗ [3,11, 12]. For this, the size-skew quality metric [15] was widely used. This metric is based on the Jacobian determinants and measures the relative size change and distortion of an element respect to an original one. It also ensures that the computed mesh is valid (positive Jacobian determinants). However, such criteria are mesh-based. To keep the proposed formulation free from mesh connectivity, alternative metrics derived from continuum mechanics are introduced. The main purpose of these metrics is to evaluate the quality of the deformation field generated ⇒ by the RBF interpolation. These are based on the (infinitesimal) strain tensor ε defined as: ⇒ 1 ⇒ ⇒ ε( x) = ∇ s( x) + ∇ s( x)T (25) 2
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
with, differentiating Eqn. (1), ⇒
∇ s( x) =
Nc j=1
x − xc j ) φ (d( x − xc j )) γ j ⊗ ∇d(
355 7
(26)
Rotating the reference axes to the directions coincident with the principal strain directions, i.e. aligned with the ⇒ eigenvectors of ε, the strain tensor becomes: εI 0 0 ⇒ (27) ε( x) = 0 εII 0 0 0 εIII ⇒
where εI ≤ εII ≤ εIII are the ordered eigenvalues of ε. By analogy with Mohr’s strain circles [16,17], two metrics are introduced:
• The volumetric strain, defined as 13 (εI + εII + εIII ), quantifies volume variations around position x. Dilatation (resp. compression) is associated with positive (resp. negative) values. For small deformations, it approximates the relative variation in volume. • The maximum shear strain, defined as (εIII − εI ), quantifies orthogonality variations around position x. For small deformations, it approximates the largest angular variation between planes initially orthogonal, crossing at x. Displacement fields exhibiting low volumetric strains and low shear strains are expected to yield morphed meshes of quality close to the quality of the initial ones. Acceptable values of these metrics are a necessary but not sufficient condition to perform the update position step (Alg. 1 lines 6-7). Then if the deformed mesh is computed, the size-skew metric may be used to check the mesh quality. 3. Applications 3.1. Global setup Information pertaining to all the configurations follow. 3.1.1. Distance measuring function The unbiased Euclidean norm is used to measure the distance between two positions. It is furthermore scaled by a reference length r0 , related to the mesh size and to the imposed displacements amplitude. The influence of each control point is thus isotropic, iso-distance surfaces from a point being spherical in ( e1 , e2 , e3 ) space. Therefore: 1 d( u) = u (28) r0 1 u r02 d( u) , if d( u) > 0 ∇d( u) = (29) 0, if d( u) = 0
3.1.2. Radial basis function Several RBF expressions are available [18,19]. Based on previous studies [3,5,11], Wendland’s C 2 compactly supported function is used: (4d + 1)(1 − d)4 , if d < 1 φ(d) = (30) 0, if d ≥ 1 −20d(1 − d)3 , if d < 1 (31) φ (d) = 0, if d ≥ 1 φ is monotonically decaying from a control point as d increases, and it is null if the scaled distance is larger than 1. So r0 is referenced thereafter as the support radius, limiting in space the volume influenced by the control points.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
356 8
3.1.3. Linear system solver As only a few hundreds of control points are considered, standard LU factorization is used to solve the notsymmetrical linear system defined by Eqn. (22). 3.2. Deforming panel Algorithm 3 is firstly applied to a 2D configuration modeling a shock tube fluid-structure interaction case in which the shock wave impinges on a cross-flow flexible panel, initially at rest [20,21]. As the panel root is clamped, its tip moves back and forth after shock wave impact. Figure 1(b) presents the subdomain including the panel. Upper and lower boundaries are rigid walls. The shock wave enters this region through the left boundary. The panel initially straight and vertical (opened symbols) is firstly bending rightward (filled symbols). The challenge is to keep the upper wall flat while the panel tip is free to move in both directions and the tip clearance is varying in time. The influence of the initial distance δ between panel tip and upper wall on the predicted displacement field s is studied through two configurations. It is expected that the smaller δ is, the more constrained s should be. 3.2.1. Setup The 2D coordinate system used is such that ( e1 , e2 ) = ( e x , ey ). The panel is 1 mm thick and 50 mm high. Its mean camber line at rest is centered on x = 0.015 m (see Fig. 1); the modal displacement of a point of this line is defined as: y s(x, y) = 3y (32) −x The largest displacement at tip is thus smax = 7.83 mm. Panel thickness remains constant during the bending motion. Panel tip is rotating of about 17 deg. Support radius r0 is set to 40 mm, about 5 times smax . Greedy algorithm tolerance is set to 10−4 smax . The four corners, used to initialize the greedy algorithm, are prescribed as fixed control points. The bottom wall is sampled every 2 mm with fixed control points. The panel skin is sampled every 0.5 mm with moving control points. The left/top/right bounding lines are sampled every 1 mm with sliding control points.
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
0.02
100
Normalized normal displacement [-]
0.07
y [m]
y [m]
3.2.2. Larger tip gap configuration
0.02
0.01
0.01
0.00
0.00
−0.01
−0.01
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
(a) Without slip condition
(b) With slip condition
Left
Top
Right
10−1
10−2
10−3
10−4
10−5
10−6 0.00
0.05 0.10 0.15 Curvilinear abscissa [m]
0.20
(c) s · n/smax along boundaries
Fig. 1. Panel configuration with 15 mm tip gap - (a,b) Outline of the deformed subdomain (open symbols and dashed lines represent initial positions; filled symbols and solid lines represent deformed positions; circles ◦ mark moving points, diamonds♦ fixed points and squares sliding points) (c) Displacement normal to the bounding planes (solid line: without slip condition (a); dashed line: with slip condition (b)).
Figures 1 and 2 present the results achieved for an initial tip gap height δ = 15 mm, the panel motion increasing the tip gap by 15%. As φ(δ) 0.4, it is expected that the upper wall is moderately modified in this configuration.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361
357
S. Aubert et al. / Procedia Engineering 00 (2017) 000–000 0.070
0.065
9
0.07
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
y [m]
y [m]
0.055
0.050
y [m]
0.060
0.02
0.02 0.01
0.01 0.045
0.00
0.00
0.040
−0.01
−0.01
0.035
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
0.030 0.005
0.010
0.015 0.020 x [m]
0.025
0.030
(a) Deformed mesh in tip region
−0.3
−0.2
−0.1 0.1 Volumetric Strain [-]
0.2
(b) Variations in cell volume
0.3
0
9
18 27 Shear Strain [deg]
36
45
(c) Variations in cell orthogonality
Fig. 2. Panel configuration with 15 mm tip gap - Displacement field with slip condition - (a) Test mesh - (b) Volumetric strain (levels between -5% and 5% are not plotted) - (c) Maximum shear strain (levels below 5 deg are not plotted).
The outline of the deformed domain is drawn Fig. 1(a) when no sliding points are specified; are plotted also the 22 fixed points and the 48 moving points selected by the greedy algorithm. It could be seen that the panel tip moving rightward and downward is dragging down the unconstrained upper wall, creating a depression in the top bounding line. The left one is also slightly modified. To quantify these side effects, the displacement normal to the left/top/right bounding planes is plotted Fig. 1(c), normalized by smax , as a function of the curvilinear abscissa measured from the bottom left corner. The green vertical dashed lines mark the position of upper left and upper right corners. On left and top boundaries, the largest errors are thus of the order of 0.1smax . An insignificant displacement is recorded close to the bottom right corner. This is a side effect of the fixed points along the bottom wall. The outline of the deformed domain is drawn Fig. 1(b) when sliding points are added along the left/top/right bounding lines; are plotted also the 22 fixed points, the 48 moving points and the 28 sliding points selected by the greedy algorithm. No deformation, out of the imposed panel motion, could be seen. For this setup, Fig. 1(c) shows that the largest error on the s · n/smax indicator is below the greedy algorithm tolerance, marked by the red horizontal dashed line. The slip condition constraining the displacement field, the total number of selected control points is accordingly increased from setups (a) to (b). To visualize the achieved space deformation, the displacement field with slip condition is applied to an uniform cartesian test mesh. Grid spacing is set to 0.5 mm in both directions. Figure 2(a) shows the deformed mesh into the tip gap. The black vertical dashed line marks the panel initial position. Near the panel tip, a solid rotation coupled to the imposed translation is predicted, while along the upper wall, only a horizontal translation is predicted. This horizontal motion is smoothened and dampened moving away from panel tip up to upper wall. To quantify this displacement field quality, volumetric strain and shear strain fields are respectively plotted Fig. 2(b) and Fig. 2(c) over the same test mesh. Close to zero values are predicted in the vicinity of the panel skin. Such values are associated with near rigid body motions, with no shear, nor compression, nor dilatation. As a consequence for example, the quality of a boundary layer mesh attached to the panel would be preserved, in terms of orthogonality and aspect ratio. It is to notice that such a behavior is not predicted along the slip walls, where only the normal motion is constrained. The panel moving rightward, the mesh is subjected to an expansion on its left side, as indicated by volumetric strain positive values Fig. 2(b). The maximum value, slightly above 20%, is in agreement with the 15% increase of tip gap height; there is also a contribution to this extremum of the tip horizontal displacement. On the right side, the relative compression remains below 20% in absolute value. The shear strain is not as symmetrically distributed Fig. 2(c), the maximum values remaining however below 36 deg on both sides. On the left side, the extremum seems correlated with the slightly anisotropic dilatation, driven there by the horizontal motion. On the right side, the extremum is closer to the upper wall (compared to the highest compression region), where the volumetric strain is small.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
358 10
Figure 2(a) reveals that in this region, cells are elongated vertically and compressed horizontally. The overall volume is thus slightly modified there while the highest shear is significant. As a conclusion, to counter balance the 17 deg rotation of the panel tip and the 15% widening of the tip gap, the cells above the panel tip are subjected in this configuration to a mild anisotropic dilatation and to a mild shear.
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
0.02 0.01
100
Normalized normal displacement [-]
0.07
y [m]
y [m]
3.2.3. Smaller tip gap configuration
0.02 0.01
0.00
0.00
−0.01
−0.01
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
(a) Without slip condition
Left
Top
Right
10−1
10−2
10−3
10−4
10−5
10−6 0.00
0.05 0.10 0.15 Curvilinear abscissa [m]
(c) s · n/smax along boundaries
(b) With slip condition
Fig. 3. Panel configuration with 5 mm tip gap - (a,b) Outline of the deformed subdomain (open symbols and dashed lines represent initial positions; filled symbols and solid lines represent deformed positions; circles ◦ mark moving points, diamonds♦ fixed points and squares sliding points) (c) Displacement normal to the bounding planes (solid line: without slip condition (a); dashed line: with slip condition (b)).
0.070
0.065
0.07
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
y [m]
y [m]
0.055
0.050
y [m]
0.060
0.02
0.02 0.01
0.01 0.045
0.00
0.00
0.040
−0.01
−0.01
0.035
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
0.030 0.005
0.010
0.015 0.020 x [m]
0.025
0.030
(a) Deformed mesh in tip region
−0.3
−0.2
−0.1 0.1 Volumetric Strain [-]
0.2
(b) Variations in cell volume
0.3
0
9
18 27 Shear Strain [deg]
36
45
(c) Variations in cell orthogonality
Fig. 4. Panel configuration with 5 mm tip gap - Displacement field with slip condition - (a) Test mesh - (b) Volumetric strain (levels between -5% and 5% are not plotted) - (c) Maximum shear strain (levels below 5 deg are not plotted).
To analyze the influence of tip gap height δ on the predicted displacement field, the upper wall is lowered by 10 mm. Figures 3 and 4 present the results achieved then for δ = 5 mm, the panel motion increasing the tip gap by 45%. As φ(δ) 0.9, the upper wall should be significantly modified in this configuration. The outline of the deformed domain is drawn Fig. 3(a) when no sliding points are specified; are plotted also the 22 fixed points and the 48 moving points selected by the greedy algorithm. As expected, the depression in the top bounding line is more severe compared to Fig. 1(a). The s · n/smax indicator presented Fig. 3(c) shows that the largest error there is of the order of smax , while it is lowered to 0.01smax along the left bounding line. When sliding points are
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
100
Original Without slip condition With slip condition
75
Percentage of cells [%]
Percentage of cells [%]
50 25 0
0.90
0.92
0.94 0.96 Size-skew metric [-]
(a) Larger tip gap quality metric
0.98
1.00
100
359 11
Original Without slip condition With slip condition
75 50 25 0
0.90
0.92
0.94 0.96 Size-skew metric [-]
0.98
1.00
(b) Smaller tip gap quality metric
Fig. 5. Size-skew metric histograms - (a,b) The original mesh (black) have a maximal quality for all the elements. Morphing, without (white) and with (grey) slip condition, highlight a slightly reduction of the quality.
added along the left/top/right bounding lines, 22 fixed points, 47 moving points and 31 sliding points are selected by the greedy algorithm (see Fig. 3(b)). The numbers of selected control points are thus very marginally modified by the tip gap height, with −1 moving and +3 sliding points relatively to the δ = 15 mm configuration. Figure 4(a) presents the deformed mesh into the tip gap. Compared to Fig. 2(a), it shows that the horizontal translation applied to the cells above the panel tip is almost uniform. Looking at the volumetric strain and the shear strain fields Fig. 4(b) and Fig. 4(c), it appears that the near rigid body motion close to the panel skin is preserved, with levels close to zero. The most striking differences with Fig. 2 are, above the panel tip, a significant expansion conjugated with almost no shear. The volumetric strain maximum value, above 30%, is there in agreement with the relative increase of tip gap height. A mild shear strain at the same spot indicates near-isotropic deformation. On the contrary, the maximum shear strain on the right side is beyond 54 deg. Compared to the δ = 15 mm configuration, the vertical elongation there is more severe, increasing the deformation anisotropy. As a conclusion, to counter balance the 17 deg rotation of the panel tip and the 45% widening of the tip gap, the cells above the panel tip are subjected in this configuration to a strong dilatation, quite isotropic, coupled to a slight shear. 3.2.4. Quality metric analysis The size-skew metric histograms for the two configurations are shown on Fig. 5. In both cases, the initial quality is maximum (equal to 1) for all elements of the mesh. Indeed, the test mesh is an uniform cartesian grid with a grid spacing set to 0.5 mm in both directions. After the deformation without slip condition, 80% of the elements still have an excellent quality because the impacted region is confined in space. In the impacted area, the worst quality is about 0.90 for the larger tip gap (Fig. 5(a)) and 0.93 for the smaller one (Fig. 5(b)). The deformation using slip condition slightly increases the overall quality. 3.3. Parameterized variations of fan tip To illustrate 3D applications, the tip gap of a generic aero engine fan is considered (see Fig. 6(a,b)). To analyze in rotation varying conditions, it is subjected to conjugated parametrized variations of tip clearance (up to 1 mm) and stagger angle (in the range ±5 deg). Thereafter, a tip gap reduction of 0.5 mm and a rotation of 2 deg centered on the trailing edge are presented as preliminary results. The initial blade tip is drawn Fig. 6(a) with the control points locations. On Fig. 6(b), the after deformation surface is superimposed on the initial blade tip outline. The largest displacement is smax = 4.14 mm, located at the leading edge. Support radius r0 is set to 12 mm, about 3 times smax . Greedy algorithm tolerance is set to 10−3 smax . Moving control points are located only along the blade top profile. The 3D coordinate system used is such that ( e1 , e2 , e3 ) = ( eR , eRθ , ez ), with ez the rotation axis, eR and eRθ the associated radial and azimuthal directions. In this basis, the cylindrical shroud is planar, and Rθ is well adapted to distance measurement in azimuthal direction compared with R in radial direction. Figure 6(c) shows a meridional
12 360
S.St´ephane Aubert et al. / Procedia 00 (2017) 203 000–000 Aubert et al. / Engineering Procedia Engineering (2017) 349–361
0.2610 0.2605 0.2600
r[m]
0.2595 0.2590 0.2585 0.2580 0.2575 0.2570 −0.10 −0.08 −0.06 −0.04 −0.02 z[m]
(a) Front view with control points
(b) Morphed front view
0.00
0.02
0.04
(c) Meridional view
Fig. 6. Fan tip deformation - (a) Moving points (black) and sliding points (white) location on the blade (grey) and the shroud (alpha grey) before morphing - (b) Blade tip before (black line) and after (grey surface) morphing - (c) Outline of shroud and blade (black dashed lines represent initial positions; red dashed lines represent deformed positions without slip condition; solid black lines represent deformed positions with slip condition).
6 Percentage of cells [%]
0.260
R[m]
0.258 0.256 0.254 0.252 0.250 −0.010
−0.005
0.000 Rθ[m]
0.005
(a) Mesh before (grey) and after (black) morphing
0.010
Initial With slip condition
4
2
0 0.0
0.2
0.4 0.6 Skew metric [-]
0.8
1.0
(b) Skew metric histograms before (grey) and after (black) morphing
Fig. 7. Fan tip deformation - (a) Mid-chord axial view - (b) Quality metric histograms, evaluated for the mesh in the coordinate system ( e x , ey , ey )
view (z, R) of the blade tip, initially at R = 0.258 m, and the shroud at R = 0.260 m (black dashed line). The after deformation outline, without slip condition applied on shroud, is marked by a red dashed line. The trailing edge moves in radial direction of about 0.5 mm only, whereas the leading edge displacement is beyond 1 mm along eR . Indeed, the blade tip rotation centered on the trailing edge in ( e x , ey , ez ) space brings the leading edge closer to the shroud with increasing stagger angles. r0 being 6 times larger than tip clearance, the displacement without slip condition on shroud is significant. Figure 7(a) shows one every four nodes of the mesh cut at mid-chord (z = −0.03 m) with slip condition applied. The initial grid is drawn in grey, while the after deformation one is in black. 303 moving points and 266 sliding points are selected by the greedy algorithm. After morphing, the shroud radius remains constant; the block meshing the tip clearance appears smoothly compressed radially and shifted in azimuth. It is to notice that by using RBF with compact support, the modified region is confined in space. Values above the display threshold set to 5 deg are limited to a disc centered on blade tip, of about 9 mm in radius, i.e. 32 r0 . Figure 7(b) shows the skew metric [15] for the initial mesh (grey) and the deformed one (black) evaluated for the mesh in the coordinate system ( e x , ey , ey ). The size-skew metric was not employed due to the inadequacy of the relative size metric to evaluate the quality of this kind of meshes. Indeed, the mesh is formed by a hexaedral structured grid, and around the blade, cells are elongated and twisted to discretize the shape. However it is not prejudicial for the
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
361 13
CFD solution. For this same reason, only 64% of cells have a skew quality larger than 0.8 for the initial mesh. The mesh has an acceptable quality considering 95% of the cells have a quality larger than 0.50. The morphing degrades slightly the quality, such as 54% and 89% of the cells have a quality larger than 0.8 and 0.5, respectively. 4. Conclusion A new formulation extending the well-known interpolation scheme based on Radial Basis Functions to morph meshes has been detailed. It introduces, besides the classical moving control points, a subset of sliding control points for which conditions of zeroed normal displacement and zeroed tangential contributions are prescribed in a casespecific vector space. Multiple slip planes are handled together, while the linear property, linking the deformation response to the known displacements, underlying the RBF formalism is preserved. Test cases analysis was based on deformation quality metrics derived from the displacement field gradient, keeping the approach free from connectivity information. It shows that the proposed RBF formulation is well adapted to yield smooth displacements fields whatever the distance between moving and sliding control points is. Further works will focus on 3D configurations for which the radial basis function expression, coupled to the support radius value, are key components to find a compromise between efficiency, robustness and displacement field smoothness. References [1] S. Jakobsson, O. Amoignon, Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization, Computers & Fluids 36 (2007) 1119–1136. [2] M. Biancolini, I. Viola, M. Riotte, Sails trim optimisation using cfd and rbf mesh morphing, Computers & Fluids 93 (2014) 46–60. [3] A. de Boer, M. Van der Schoot, H. Bijl, Mesh deformation based on radial basis function interpolation, Computers & Structures 85 (2007) 784–795. [4] T. Rendall, C. Allen, Unified fluid–structure interpolation and mesh motion using radial basis functions, International Journal for Numerical Methods in Engineering 74 (2008) 1519–1559. [5] Y. Rozenberg, S. Aubert, G. Benefice, Fluid structure interaction problems in turbomachinery using rbf interpolation and greedy algorithm, in: ASME Turbo Expo 2014 Turbine Technical Conference and Exposition, number 26735 in GT2014, ASME, 2014, p. V02BT39A033. [6] J. T. Batina, Unsteady euler airfoil solutions using unstructured dynamic meshes, AIAA journal 28 (1990) 1381–1388. [7] C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional springs for two-dimensional dynamic unstructured fluid meshes, Computer methods in applied mechanics and engineering 163 (1998) 231–245. [8] E. J. Nielsen, W. K. Anderson, Recent improvements in aerodynamic design optimization on unstructured meshes, AIAA journal 40 (2002) 1155–1163. [9] T. W. Sederberg, S. R. Parry, Free-form deformation of solid geometric models, Computer Graphics 20 (1986) 151–160. [10] J. A. Samareh, Multidisciplinary Aerodynamic-Structural Shape Optimization Using Deformation (MASSOUD), Technical Report 4911, AIAA, 2000. [11] B. Mocanu, R. Tapu, E. Tapu, Mesh deformation with hard constraints, in: 2013 International Symposium on Signals, Circuits and Systems, ISSCS, IEEE, 2013, pp. 1–4. [12] A. K. Michler, Aircraft control surface deflection using rbf-based mesh deformation, International Journal for Numerical Methods in Engineering 88 (2011) 986–1007. [13] R. Schaback, H. Wendland, Adaptive greedy techniques for approximate solution of large rbf systems, Numerical Algorithms 24 (2000) 239–254. [14] T. C. Rendall, C. B. Allen, Efficient mesh motion using radial basis functions with data reduction algorithms, Journal of Computational Physics 228 (2009) 6231–6249. [15] P. M. Knupp, Algebraic mesh quality metrics for unstructured initial meshes, Finite Elements in Analysis and Design 39 (2003) 217–241. [16] R. H. Parry, Mohr circles, stress paths and geotechnics, Spon Press, 2004, pp. 34–37. [17] F. P. Beer, E. R. Johnston, J. T. DeWolf, Mechanics of materials, Tata McGraw-Hill, 2004, pp. 470–477. [18] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in computational Mathematics 4 (1995) 389–396. [19] M. D. Buhmann, Radial basis functions, Acta Numerica 2000 9 (2000) 1–38. [20] J. Giordano, G. Jourdan, Y. Burtschell, M. Medale, D. E. Zeitoun, L. Houas, Shock wave impacts on deforming panel, an application of fluid-structure interaction, Shock Waves 14 (2005) 103–110. [21] G. Benefice, Y. Rozenberg, S. Aubert, P. Ferrand, F. Thouverez, A partitioned strong coupling procedure for simulation of shock wave/flexible structure interaction, in: ASME 2014 4th Joint US-European Fluids Engineering Division Summer Meeting, number 21624 in FEDSM2014, ASME, 2014, p. V01BT12A006.
Available online at www.sciencedirect.com Available online at www.sciencedirect.com
ScienceDirect Procedia Engineering 20300 (2017) 349–361 Procedia Engineering (2017) 000–000 Procedia Engineering 00 (2017) 000–000
www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia
26th 26th International International Meshing Meshing Roundtable, Roundtable, IMR26, IMR26, 18-21 18-21 September September 2017, 2017, Barcelona, Barcelona, Spain Spain
Planar Planar Slip Slip Condition Condition For For Mesh Mesh Morphing Morphing Using Using Radial Radial Basis Basis Functions Functions a,∗ a,b a a St´ Rendu St´eephane phane Aubert Auberta,∗,, Franck Franck Mastrippolito Mastrippolitoa,b ,, Quentin Quentin Rendua ,, Martin Martin Buisson Buissona ,, Fr´ Fr´eed´ d´eeric ric b Ducros b Ducros a Univ. a Univ.
Lyon, Ecole Centrale de Lyon, Laboratoire de M´ecanique des Fluides et d’Acoustique UMR5509, Lyon, F-69134, France Lyon, Ecole Centrale de Lyon, Laboratoire de M´ecanique des´ Fluides et d’Acoustique UMR5509, Lyon, F-69134, France b Univ. Alpes, CEA, LITEN, DTBH, Laboratoire Echangeur et R´eacteurs, F-38000 Grenoble, France b Univ. Grenoble ´ Grenoble Alpes, CEA, LITEN, DTBH, Laboratoire Echangeur et R´eacteurs, F-38000 Grenoble, France
Abstract Abstract An original approach to morph meshes including slip planes is proposed. Extending the classical interpolation scheme based An original approach to morph meshes including slip planes is proposed. Extending the classical interpolation scheme based on Radial Basis Functions (RBF), it preserves the linear property, linking the deformation response to the known displacements, on Radial Basis Functions (RBF), it preserves the linear property, linking the deformation response to the known displacements, underlying the RBF formalism. No post-correction stage, nor extra not-moving control points, are required to enforce in plane underlying the RBF formalism. No post-correction stage, nor extra not-moving control points, are required to enforce in plane displacements. Accuracy and robustness are evaluated for a 2D configuration where a shock-wave impinges on a cross-flow displacements. Accuracy and robustness are evaluated for a 2D configuration where a shock-wave impinges on a cross-flow flexible panel. Analysis is based on deformation quality metrics derived from the displacement field gradient, keeping the proposed flexible panel. Analysis is based on deformation quality metrics derived from the displacement field gradient, keeping the proposed approach free from connectivity information. 3D applications are briefly illustrated by the tip of an aero engine fan subjected to approach free from connectivity information. 3D applications are briefly illustrated by the tip of an aero engine fan subjected to conjugated parameterized variations of tip clearance and stagger angle. conjugated parameterized variations of tip clearance and stagger angle. c 2017 The Authors. Published by Elsevier Ltd. c 2017 2017 The The Authors. Authors. Published Published by by Elsevier Elsevier Ltd. Ltd. © Peer-review under responsibility of the scientific committee of the 26th International Meshing Roundtable. Peer-review under responsibility of of the the scientific scientific committee committee of of the the 26th 26th International International Meshing Meshing Roundtable. Roundtable. Keywords: radial basis function mesh deformation, slip condition, CFD Keywords: radial basis function mesh deformation, slip condition, CFD
1. Introduction 1. Introduction In many CFD-applications, the mesh has to be moved after or during the simulation process. This feature is In many CFD-applications, the mesh has to be moved after or during the simulation process. This feature is typically implemented in optimization [1,2] and Fluid-Structure Interaction [3–5] solvers. A main issue is that the typically implemented in optimization [1,2] and Fluid-Structure Interaction [3–5] solvers. A main issue is that the mesh displacement is highly strained by grid-quality constrains (skewness, aspect ratio, volume positivity), boundary mesh displacement is highly strained by grid-quality constrains (skewness, aspect ratio, volume positivity), boundary conditions (wall, periodic, symmetric) or prescribed surface deformations. conditions (wall, periodic, symmetric) or prescribed surface deformations. Morphing a mesh using springs analogy [6,7], pseudo-solid approach [8] or free form deformation [9,10] is a well Morphing a mesh using springs analogy [6,7], pseudo-solid approach [8] or free form deformation [9,10] is a well studied subject. But to propagate displacements known only at clouds of control points to the whole volume mesh, studied subject. But to propagate displacements known only at clouds of control points to the whole volume mesh, the Radial Basis Function (RBF) interpolation method is more and more popular, after its introduction for external the Radial Basis Function (RBF) interpolation method is more and more popular, after its introduction for external aeroelasticity applications [1,4]. Based only on initial position and displacement of these control points, it defines aeroelasticity applications [1,4]. Based only on initial position and displacement of these control points, it defines by solving a linear system a continuous displacement field that requires no connectivity information. It is thus very by solving a linear system a continuous displacement field that requires no connectivity information. It is thus very flexible to deal with structured or unstructured grids, meshing fluid or solid domains. Besides, its memory requirement flexible to deal with structured or unstructured grids, meshing fluid or solid domains. Besides, its memory requirement ∗ ∗
Corresponding author. Tel.: +3-347-218-6741 Corresponding author. Tel.: +3-347-218-6741 E-mail address: [email protected], [email protected] E-mail address: [email protected], [email protected]
c 2017 The Authors. Published by Elsevier Ltd. 1877-7058 c 2017 The Authors. Published by Elsevier Ltd. 1877-7058 Peer-review©under responsibility of the scientific committee of the 26th International Meshing Roundtable. 1877-7058 2017responsibility The Authors. Published Elsevier of Ltd. Peer-review under of the scientificbycommittee the 26th International Meshing Roundtable. Peer-review under responsibility of the scientific committee of the 26th International Meshing Roundtable. 10.1016/j.proeng.2017.09.819
Aubert et al. /Engineering Procedia Engineering 203 (2017) 349–361 S. St´ephane Aubert et al. / Procedia 00 (2017) 000–000
350 2
Nomenclature b d ei n Nc Nm Ns r0 s smax smi S t x xci xmi x si γi Γ
Second unit tangent vector, 3i=1 bi ei (∈ R3 ). Distance measuring function (R3 → R+ ). ith basis vector (∈ R3 ). Unit normal vector, 3i=1 ni ei (∈ R3 ). Number of control points (∈ N). Number of moving points (∈ N). Number of sliding points (∈ N). Support radius (∈ R). Displacement function (R3 → R3 ). Largest displacement magnitude (∈ R). Displacement vector of ith moving point (∈ R3 ). Vector (∈ RNm ). First unit tangent vector, 3i=1 ti ei (∈ R3 ). Position vector (∈ R3 ). Position vector of ith control point (∈ R3 ). Position vector of ith moving point (∈ R3 ). Position vector of ith sliding point (∈ R3 ). Weight vector of ith control point (∈ R3 ). Vector (∈ RNc ).
δ ε φ φ Φm Φs Φt 0s 0 s×m 0m×c I s×s u · v u ∧ v u ⊗ v u ∇u ⇒
⇒
∇ u AT ∗
Tip gap height (∈ R). Strain tensor (∈ R3×3 ). Radial basis function (R+ → R). Derivative of φ (R+ → R). Matrix (∈ RNm ×Nc ). Matrix (∈ RNs ×Nc ). Matrix (∈ RNs ×Nc ). Null vector (∈ RNs ). Null matrix (∈ RNs ×Nm ). Null matrix (∈ RNm ×Nc ). Identity matrix (∈ RNs ×Ns ). Dot product, 3i=1 ui vi ((R3 , R3 ) → R). Cross product, [u j vk − uk v j ]i ((R3 , R3 ) → R3 ). Outer product, [ui v j ]i, j ((R3 , R3 ) → R3×3 ). √ Euclidean norm, u · u (R3 → R). Gradient operator, [∂u/∂xi ]i (R → R3 ).
Gradient operator, [∂ui/∂x j ]i, j (R3 → R3×3 ). Transpose operator, [a ji ]i, j (R3×3 → R3×3 ). After deformation state. Vector components written in ( e x , ey , ez ) space.
is mesh size independent [11]. However standard RBF formulation is isotropic. Mandatory geometrical properties, like planar shapes to define rotor/stator interfaces, are likely to be spoiled. Furthermore, using functions with compact support to localize mesh deformations, action radii should be at least three times larger than the control points displacement to smooth these over a long enough distance. It is thus tricky to setup RBF data for confined space configurations, like blade tip clearance and under-platform cavities. To circumvent such drawbacks, some classical approaches rely on adding enclosing boxes of not-moving control points [2] or implementing predictor/corrector strategies based on directional damping functions [1,12]. All try and preserve selected space regions by cancelling some undesired movements after the RBF setup. Unfortunately, there is no (retro)action of such corrections on the interpolation coefficients; beside these require many user defined extra data, which are case dependent. Thereafter, the proposed approach is firstly presented in section 2. It is based on the formulation of a slip condition, along one or more planes, for subsets of control points instead of specifying their whole displacement. Prescribed moving and sliding control points are taken into account together in the linear system defining the interpolation coefficients. So, no post-correction stage is required, nor extra not-moving control points. Deformation quality metrics inspired by continuum mechanics are also introduced. Based on the displacement field gradient, these keep the proposed approach free from connectivity information. Then, accuracy and robustness are analyzed in section 3.2 through a 2D configuration where a shock-wave impinges on a cross-flow flexible panel. 3D applications are briefly illustrated in section 3.3 by the tip of an aero engine fan subjected to conjugated parameterized variations of tip clearance and stagger angle.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
351 3
2. Formulation 2.1. Displacement field interpolation RBF-based mesh deformation is a node-by-node scheme free from mesh connectivity. It propagates by interpolation the known motion of a cloud of control points to the mesh nodes. The interpolated function is the displacement s of a node, defined as the weighted sum of radial basis functions by: s( x) =
Nc j=1
φ(d( x − xc j )) γ j
(1)
where x is the initial (i.e. before deformation) node position, Nc is the number of control points, φ is the chosen radial basis function, d is the chosen function measuring distance between two positions, xc j and γ j are respectively the initial position and weight of the jth control point. It is to notice that only initial positions appear explicitly in Eqn. (1). The prescribed displacement of control points and possibly other constraints are embedded in the weights { γi }1≤i≤Nc . When these are known, s is fully defined. Thereafter, the set of control points is split in two. Different conditions are associated to each subset to close the problem defining the weights. Let Nm ≤ Nc be the number of moving control points, i.e. with fully prescribed motion. Fixed, not moving, control points are included in this subset with zeroed displacement. Let N s = Nc − Nm be the number of sliding control points, i.e. with motion constrained to remain in a plane. The set of control points is assumed to be ordered, starting with the moving ones with index ranging from 1 to Nm , followed by the sliding ones with index ranging from Nm + 1 to Nm + N s = Nc . 2.2. Interpolation condition On one hand, the displacement s must satisfy the interpolation condition for each moving control point: ∀i ∈ [1, Nm ], ⇔ ⇔
∀k∈[1,3]
s( xmi ) = smi Nc φ(d( xmi − xc j )) γ j = smi
(2) (3)
j=1
Nc j=1
φ(d( xmi − xc j )) γ j · ek = smi · ek
(4)
where smi is the known displacement of the ith moving point initially at position xmi , and ( e1 , e2 , e3 ) is a basis of the 3D Euclidean space. Equation (4) is rewritten as: (5) ∀k ∈ [1, 3], Φm Γk = S k where matrix Φm and vectors Γk and S k are defined as: Φm = φ(d( xmi − xc j )) 1≤i≤Nm ,1≤ j≤Nc = φ(d( xci − xc j )) 1≤i≤Nm ,1≤ j≤Nc
∀k ∈ [1, 3],
Γk = γi · ek 1≤i≤Nc S k = smi · ek 1≤i≤Nm
(6) (7) (8) (9)
It is to notice that if Nm = Nc , Φm is a square matrix and the components of { γi }1≤i≤Nc are simply solution of the three linear systems defined by Eqn. (5). It is the classical RBF formulation where the Γk vectors are independent.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
352 4
2.3. Planar slip condition On the other hand, the displacement s must satisfy three scalar conditions, associated to the normal and tangential directions of the slip planes. For the sake of simplicity thereafter, with no loss of generality, only one such plane is considered. Let n be its unit normal vector and t a unit tangent vector, such that t · n = 0. Another unit tangent vector b is simply defined as b = n ∧ t. 2.3.1. Zeroed normal displacement The first condition on s states that its normal component s · n is zero for each sliding control point: ∀i ∈ [1, N s ], s( x si ) · n = 0 3 Nc e =0 n φ(d( x si − xc j )) γ j · ⇔ k k j=1 k=1 Nc 3 =0 nk φ(d( x − x )) γ · e ⇔ si cj j k j=1 k=1
(10) (11)
(12)
where x si is the initial position of the ith sliding point and (n1 , n2 , n3 ) are the components of n in ( e1 , e2 , e3 ) basis. Equation (12) is rewritten as: 3 nk Φ s Γk = 0 s (13) k=1
Ns
s
where 0 s is the null vector in R and matrix Φ is defined as: Φ s = φ(d( x si − xc j )) 1≤i≤N s ,1≤ j≤Nc = φ(d( xc(Nm +i) − xc j )) 1≤i≤N s ,1≤ j≤Nc
(14) (15)
2.3.2. Zeroed tangential contributions The second condition on s states that no sliding control point is contributing to the displacement function in the tangential direction t: ∀i ∈ [Nm + 1, Nm + N s ],
γi · t = 0 3 e =0 t ⇔ γi · k k
(16) (17)
k=1
⇔
3 k=1
tk γi · ek = 0
(18)
where (t1 , t2 , t3 ) are the components of t in ( e1 , e2 , e3 ) basis. Equation (18) is rewritten as: 3 tk Φt Γk = 0 s
(19)
k=1
where matrix Φt is defined as:
Φt = 0 s×m I s×s ∈ RNs ×Nc
(20)
with 0 s×m the null matrix in RNs ×Nm and I s×s the identity matrix in RNs ×Ns . The third and last condition on s states that no sliding control point is contributing to the displacement function in the other tangential direction b, i.e. substituting b to t in Eqn. (19): 3 k=1
bk Φt Γk = 0 s
(21)
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
353 5
where (b1 , b2 , b3 ) are the components of b in ( e1 , e2 , e3 ) basis. 2.4. Weights calculation Therefore, the components of { γi }1≤i≤Nc are solution of the (3Nc ) × (3Nc ) linear system assembled from Eqn. (5), Eqn. (13), Eqn. (19) and Eqn. (21): m S 1 Φ 0m×c 0m×c 0m×c Φm 0m×c S 2 Γ 0m×c 0m×c Φm 1 S 3 (22) Γ2 = n1 Φ s n2 Φ s n3 Φ s Γ 0 s 0 t Φt t Φt t Φt 3 2 3 s 1 t t t 0s b1 Φ b2 Φ b3 Φ with 0m×c the null matrix in RNm ×Nc . When multiple slip planes are considered, Eqn. (13), Eqn. (19) and Eqn. (21) must be replicated for each one of them with associated ( n, t, b) triplet.
It is to notice that the Γk vectors are coupled in Eqn. (22) by the matrix last three lines in the general case, compared to the classical RBF formulation. However if a single slip plane is considered, an uncoupled formulation is recovered in the basis ( e1 , e2 , e3 ) = ( n, t, b). The components of { γi }1≤i≤Nc are then again solution of three Nc × Nc linear systems:
∀k ∈ [2, 3],
Φm S Γ = 1 Φs 1 0s m Φ S Γ = k 0s 0 s×m I s×s k
(23) (24)
On one hand, Eqn. (23) can be interpreted as the classical RBF formulation projected on n, with the sliding points treated as fixed control points. On the other hand, Eqn. (24) can be interpreted as the classical RBF formulation projected on t or b and restricted to the first Nm control points, i.e. to the moving ones. 2.5. Algorithms Several operations are required to morph a mesh. Algorithm 1 summarizes the steps of the core one yielding the position { xi∗ }1≤i≤N of the deformed mesh nodes knowing the initial one { xi }1≤i≤N , the control points initial position { xci }1≤i≤Nc and the moving control points displacement { smi }1≤i≤Nm . After assembling and solving the linear system defined by Eqn. (22), the nodes displacement is interpolated with Eqn. (1) (line 5), then it is simply added to the nodes initial position (line 7). The number and position of the control points are key factors for the efficiency and accuracy of Alg. 1. On one hand, the interpolation/update step CPU time is proportional to N(Nc + 1) (lines 5 and 7), and computing the weights is of Nc3 complexity (line 3). On the other hand, too close control points can jeopardize the condition number of matrices Φm and Φ s , because of too small values of distance function d. Therefore a greedy procedure [13,14] is applied to build a quasi-optimal subset of the initially provided control points. It is presented in Alg. 2. The subset is by default initialized with the control point exhibiting the largest displacement smax (line 1). The idea is to apply the displacement field computed using only the control points of the current subset, to move all the provided control points, using Alg. 1 (line 3). The displacement of the subset control points is exactly interpolated. So, one control point not already chosen is exhibiting the largest error, defined from Eqn. (2) for moving points (line 5) and from Eqn. (10) for sliding points (line 6). If its error is greater than some tolerance, usually defined as a fraction of smax , it is added to the subset (line 9) and the whole process is repeated.
6 354
S. St´ephane Aubert et al. / Procedia 00 (2017) 000–000 Aubert et al. /Engineering Procedia Engineering 203 (2017) 349–361
Algorithm 1: Nodes displacement
1 2 3 4 5 6 7
Input: Nodes initial position { xi }1≤i≤N Input: Control points initial position { xci }1≤i≤Nc Input: Moving points displacement { smi }1≤i≤Nm Output: Nodes final position { xi∗ }1≤i≤N Compute the weights { γi }1≤i≤Nc : Assemble the matrix and RHS of Eqn. (22) Solve Eqn. (22) to get (Γ1 , Γ2 , Γ3 ) Interpolate the nodes displacement: c ∀i ∈ [1, N], s( xi ) = Nj=1 φ(d( xi − xc j )) γ j Update the nodes position: ∀i ∈ [1, N], xi∗ = xi + s( xi )
Algorithm 2: Control points optimization
1 2 3 4 5 6 7 8 9 10 11
Input: Control points initial position { xci }1≤i≤Nc Input: Moving points displacement { smi }1≤i≤Nm Output: Control points subset Initialize the control points subset repeat Using Alg. 1, move all the control points, input as nodes while the subset is input as control points Compute errors: ∗ ∀i ∈ [1, Nm ], eval ( xmi − xmi ) − smi ∗ ∀i ∈ [1, N s ], eval |( x si − x si ) · n| Find the point with the largest error if error > tolerance then Add this point to the control points subset end until error ≤ tolerance;
The equations involved in Alg. 1 and Alg. 2 are written in ( e1 , e2 , e3 ) vector space. In particular, normal vectors n are solely defined in this space. It can be convenient to choose it as different from the standard Cartesian coordinate system ( e x , ey , ez ). For example, a conical surface in ( e x , ey , ez ) space is planar in the associated cylindrical coordinate system. However, special care must be taken if switching between ( e x , ey , ez ) and ( e1 , e2 , e3 ) spaces is not a linear process. Algorithm 3 presents the steps to morph a mesh in such a case. Before using Alg. 2 then Alg. 1 (lines 9 and 10), all x → x) (line 3). Because the position vectors known in ( e x , ey , ez ) basis must be rewritten in ( e1 , e2 , e3 ) basis (e.g. displacement is the difference between two positions, it can not be simply converted if the relation between the two spaces is not linear. Therefore, the after deformation position of the moving points is firstly calculated (line 2). Then it is rewritten in the new basis (line 6). Finally, the displacement in ( e1 , e2 , e3 ) space is recovered by subtracting initial positions from final ones (line 8). The last step (line 12) brings back the moved nodes to the original ( e x , ey , ez ) space. Algorithm 3: Overall procedure
1 2 3 4
xi }1≤i≤N Input: Nodes initial position { xci }1≤i≤Nc Input: Control points initial position { Input: Moving points displacement { s mi }1≤i≤Nm xi∗ }1≤i≤N Output: Nodes final position { Compute the moving points final position in ( e x , ey , ez ): s ∀i ∈ [1, Nm ], x ∗ = x mi + mi mi
Switch from ( e x , ey , ez ) to ( e1 , e2 , e3 ) space: xi ∀i ∈ [1, N], eval xi from
5 6 7 8 9 10 11 12
xci eval xci from ∗ ∗ ∀i ∈ [1, Nm ], eval xmi from x mi Compute the moving points displacement in ( e1 , e2 , e3 ): ∗ ∀i ∈ [1, Nm ], smi = xmi − xmi Select an optimal subset of control points using Alg. 2 Move the nodes using Alg. 1 and this optimal subset Switch back from ( e1 , e2 , e3 ) to ( e x , ey , ez ) space: x ∗ from x ∗ ∀i ∈ [1, N], eval ∀i ∈ [1, Nc ],
i
i
2.6. Deformation quality metrics The quality of the displacement field s is classically evaluated through the quality of the morphed mesh x ∗ [3,11, 12]. For this, the size-skew quality metric [15] was widely used. This metric is based on the Jacobian determinants and measures the relative size change and distortion of an element respect to an original one. It also ensures that the computed mesh is valid (positive Jacobian determinants). However, such criteria are mesh-based. To keep the proposed formulation free from mesh connectivity, alternative metrics derived from continuum mechanics are introduced. The main purpose of these metrics is to evaluate the quality of the deformation field generated ⇒ by the RBF interpolation. These are based on the (infinitesimal) strain tensor ε defined as: ⇒ 1 ⇒ ⇒ ε( x) = ∇ s( x) + ∇ s( x)T (25) 2
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
with, differentiating Eqn. (1), ⇒
∇ s( x) =
Nc j=1
x − xc j ) φ (d( x − xc j )) γ j ⊗ ∇d(
355 7
(26)
Rotating the reference axes to the directions coincident with the principal strain directions, i.e. aligned with the ⇒ eigenvectors of ε, the strain tensor becomes: εI 0 0 ⇒ (27) ε( x) = 0 εII 0 0 0 εIII ⇒
where εI ≤ εII ≤ εIII are the ordered eigenvalues of ε. By analogy with Mohr’s strain circles [16,17], two metrics are introduced:
• The volumetric strain, defined as 13 (εI + εII + εIII ), quantifies volume variations around position x. Dilatation (resp. compression) is associated with positive (resp. negative) values. For small deformations, it approximates the relative variation in volume. • The maximum shear strain, defined as (εIII − εI ), quantifies orthogonality variations around position x. For small deformations, it approximates the largest angular variation between planes initially orthogonal, crossing at x. Displacement fields exhibiting low volumetric strains and low shear strains are expected to yield morphed meshes of quality close to the quality of the initial ones. Acceptable values of these metrics are a necessary but not sufficient condition to perform the update position step (Alg. 1 lines 6-7). Then if the deformed mesh is computed, the size-skew metric may be used to check the mesh quality. 3. Applications 3.1. Global setup Information pertaining to all the configurations follow. 3.1.1. Distance measuring function The unbiased Euclidean norm is used to measure the distance between two positions. It is furthermore scaled by a reference length r0 , related to the mesh size and to the imposed displacements amplitude. The influence of each control point is thus isotropic, iso-distance surfaces from a point being spherical in ( e1 , e2 , e3 ) space. Therefore: 1 d( u) = u (28) r0 1 u r02 d( u) , if d( u) > 0 ∇d( u) = (29) 0, if d( u) = 0
3.1.2. Radial basis function Several RBF expressions are available [18,19]. Based on previous studies [3,5,11], Wendland’s C 2 compactly supported function is used: (4d + 1)(1 − d)4 , if d < 1 φ(d) = (30) 0, if d ≥ 1 −20d(1 − d)3 , if d < 1 (31) φ (d) = 0, if d ≥ 1 φ is monotonically decaying from a control point as d increases, and it is null if the scaled distance is larger than 1. So r0 is referenced thereafter as the support radius, limiting in space the volume influenced by the control points.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
356 8
3.1.3. Linear system solver As only a few hundreds of control points are considered, standard LU factorization is used to solve the notsymmetrical linear system defined by Eqn. (22). 3.2. Deforming panel Algorithm 3 is firstly applied to a 2D configuration modeling a shock tube fluid-structure interaction case in which the shock wave impinges on a cross-flow flexible panel, initially at rest [20,21]. As the panel root is clamped, its tip moves back and forth after shock wave impact. Figure 1(b) presents the subdomain including the panel. Upper and lower boundaries are rigid walls. The shock wave enters this region through the left boundary. The panel initially straight and vertical (opened symbols) is firstly bending rightward (filled symbols). The challenge is to keep the upper wall flat while the panel tip is free to move in both directions and the tip clearance is varying in time. The influence of the initial distance δ between panel tip and upper wall on the predicted displacement field s is studied through two configurations. It is expected that the smaller δ is, the more constrained s should be. 3.2.1. Setup The 2D coordinate system used is such that ( e1 , e2 ) = ( e x , ey ). The panel is 1 mm thick and 50 mm high. Its mean camber line at rest is centered on x = 0.015 m (see Fig. 1); the modal displacement of a point of this line is defined as: y s(x, y) = 3y (32) −x The largest displacement at tip is thus smax = 7.83 mm. Panel thickness remains constant during the bending motion. Panel tip is rotating of about 17 deg. Support radius r0 is set to 40 mm, about 5 times smax . Greedy algorithm tolerance is set to 10−4 smax . The four corners, used to initialize the greedy algorithm, are prescribed as fixed control points. The bottom wall is sampled every 2 mm with fixed control points. The panel skin is sampled every 0.5 mm with moving control points. The left/top/right bounding lines are sampled every 1 mm with sliding control points.
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
0.02
100
Normalized normal displacement [-]
0.07
y [m]
y [m]
3.2.2. Larger tip gap configuration
0.02
0.01
0.01
0.00
0.00
−0.01
−0.01
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
(a) Without slip condition
(b) With slip condition
Left
Top
Right
10−1
10−2
10−3
10−4
10−5
10−6 0.00
0.05 0.10 0.15 Curvilinear abscissa [m]
0.20
(c) s · n/smax along boundaries
Fig. 1. Panel configuration with 15 mm tip gap - (a,b) Outline of the deformed subdomain (open symbols and dashed lines represent initial positions; filled symbols and solid lines represent deformed positions; circles ◦ mark moving points, diamonds♦ fixed points and squares sliding points) (c) Displacement normal to the bounding planes (solid line: without slip condition (a); dashed line: with slip condition (b)).
Figures 1 and 2 present the results achieved for an initial tip gap height δ = 15 mm, the panel motion increasing the tip gap by 15%. As φ(δ) 0.4, it is expected that the upper wall is moderately modified in this configuration.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361
357
S. Aubert et al. / Procedia Engineering 00 (2017) 000–000 0.070
0.065
9
0.07
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
y [m]
y [m]
0.055
0.050
y [m]
0.060
0.02
0.02 0.01
0.01 0.045
0.00
0.00
0.040
−0.01
−0.01
0.035
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
0.030 0.005
0.010
0.015 0.020 x [m]
0.025
0.030
(a) Deformed mesh in tip region
−0.3
−0.2
−0.1 0.1 Volumetric Strain [-]
0.2
(b) Variations in cell volume
0.3
0
9
18 27 Shear Strain [deg]
36
45
(c) Variations in cell orthogonality
Fig. 2. Panel configuration with 15 mm tip gap - Displacement field with slip condition - (a) Test mesh - (b) Volumetric strain (levels between -5% and 5% are not plotted) - (c) Maximum shear strain (levels below 5 deg are not plotted).
The outline of the deformed domain is drawn Fig. 1(a) when no sliding points are specified; are plotted also the 22 fixed points and the 48 moving points selected by the greedy algorithm. It could be seen that the panel tip moving rightward and downward is dragging down the unconstrained upper wall, creating a depression in the top bounding line. The left one is also slightly modified. To quantify these side effects, the displacement normal to the left/top/right bounding planes is plotted Fig. 1(c), normalized by smax , as a function of the curvilinear abscissa measured from the bottom left corner. The green vertical dashed lines mark the position of upper left and upper right corners. On left and top boundaries, the largest errors are thus of the order of 0.1smax . An insignificant displacement is recorded close to the bottom right corner. This is a side effect of the fixed points along the bottom wall. The outline of the deformed domain is drawn Fig. 1(b) when sliding points are added along the left/top/right bounding lines; are plotted also the 22 fixed points, the 48 moving points and the 28 sliding points selected by the greedy algorithm. No deformation, out of the imposed panel motion, could be seen. For this setup, Fig. 1(c) shows that the largest error on the s · n/smax indicator is below the greedy algorithm tolerance, marked by the red horizontal dashed line. The slip condition constraining the displacement field, the total number of selected control points is accordingly increased from setups (a) to (b). To visualize the achieved space deformation, the displacement field with slip condition is applied to an uniform cartesian test mesh. Grid spacing is set to 0.5 mm in both directions. Figure 2(a) shows the deformed mesh into the tip gap. The black vertical dashed line marks the panel initial position. Near the panel tip, a solid rotation coupled to the imposed translation is predicted, while along the upper wall, only a horizontal translation is predicted. This horizontal motion is smoothened and dampened moving away from panel tip up to upper wall. To quantify this displacement field quality, volumetric strain and shear strain fields are respectively plotted Fig. 2(b) and Fig. 2(c) over the same test mesh. Close to zero values are predicted in the vicinity of the panel skin. Such values are associated with near rigid body motions, with no shear, nor compression, nor dilatation. As a consequence for example, the quality of a boundary layer mesh attached to the panel would be preserved, in terms of orthogonality and aspect ratio. It is to notice that such a behavior is not predicted along the slip walls, where only the normal motion is constrained. The panel moving rightward, the mesh is subjected to an expansion on its left side, as indicated by volumetric strain positive values Fig. 2(b). The maximum value, slightly above 20%, is in agreement with the 15% increase of tip gap height; there is also a contribution to this extremum of the tip horizontal displacement. On the right side, the relative compression remains below 20% in absolute value. The shear strain is not as symmetrically distributed Fig. 2(c), the maximum values remaining however below 36 deg on both sides. On the left side, the extremum seems correlated with the slightly anisotropic dilatation, driven there by the horizontal motion. On the right side, the extremum is closer to the upper wall (compared to the highest compression region), where the volumetric strain is small.
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
358 10
Figure 2(a) reveals that in this region, cells are elongated vertically and compressed horizontally. The overall volume is thus slightly modified there while the highest shear is significant. As a conclusion, to counter balance the 17 deg rotation of the panel tip and the 15% widening of the tip gap, the cells above the panel tip are subjected in this configuration to a mild anisotropic dilatation and to a mild shear.
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
0.02 0.01
100
Normalized normal displacement [-]
0.07
y [m]
y [m]
3.2.3. Smaller tip gap configuration
0.02 0.01
0.00
0.00
−0.01
−0.01
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
(a) Without slip condition
Left
Top
Right
10−1
10−2
10−3
10−4
10−5
10−6 0.00
0.05 0.10 0.15 Curvilinear abscissa [m]
(c) s · n/smax along boundaries
(b) With slip condition
Fig. 3. Panel configuration with 5 mm tip gap - (a,b) Outline of the deformed subdomain (open symbols and dashed lines represent initial positions; filled symbols and solid lines represent deformed positions; circles ◦ mark moving points, diamonds♦ fixed points and squares sliding points) (c) Displacement normal to the bounding planes (solid line: without slip condition (a); dashed line: with slip condition (b)).
0.070
0.065
0.07
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
y [m]
y [m]
0.055
0.050
y [m]
0.060
0.02
0.02 0.01
0.01 0.045
0.00
0.00
0.040
−0.01
−0.01
0.035
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
−0.02 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x [m]
0.030 0.005
0.010
0.015 0.020 x [m]
0.025
0.030
(a) Deformed mesh in tip region
−0.3
−0.2
−0.1 0.1 Volumetric Strain [-]
0.2
(b) Variations in cell volume
0.3
0
9
18 27 Shear Strain [deg]
36
45
(c) Variations in cell orthogonality
Fig. 4. Panel configuration with 5 mm tip gap - Displacement field with slip condition - (a) Test mesh - (b) Volumetric strain (levels between -5% and 5% are not plotted) - (c) Maximum shear strain (levels below 5 deg are not plotted).
To analyze the influence of tip gap height δ on the predicted displacement field, the upper wall is lowered by 10 mm. Figures 3 and 4 present the results achieved then for δ = 5 mm, the panel motion increasing the tip gap by 45%. As φ(δ) 0.9, the upper wall should be significantly modified in this configuration. The outline of the deformed domain is drawn Fig. 3(a) when no sliding points are specified; are plotted also the 22 fixed points and the 48 moving points selected by the greedy algorithm. As expected, the depression in the top bounding line is more severe compared to Fig. 1(a). The s · n/smax indicator presented Fig. 3(c) shows that the largest error there is of the order of smax , while it is lowered to 0.01smax along the left bounding line. When sliding points are
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
100
Original Without slip condition With slip condition
75
Percentage of cells [%]
Percentage of cells [%]
50 25 0
0.90
0.92
0.94 0.96 Size-skew metric [-]
(a) Larger tip gap quality metric
0.98
1.00
100
359 11
Original Without slip condition With slip condition
75 50 25 0
0.90
0.92
0.94 0.96 Size-skew metric [-]
0.98
1.00
(b) Smaller tip gap quality metric
Fig. 5. Size-skew metric histograms - (a,b) The original mesh (black) have a maximal quality for all the elements. Morphing, without (white) and with (grey) slip condition, highlight a slightly reduction of the quality.
added along the left/top/right bounding lines, 22 fixed points, 47 moving points and 31 sliding points are selected by the greedy algorithm (see Fig. 3(b)). The numbers of selected control points are thus very marginally modified by the tip gap height, with −1 moving and +3 sliding points relatively to the δ = 15 mm configuration. Figure 4(a) presents the deformed mesh into the tip gap. Compared to Fig. 2(a), it shows that the horizontal translation applied to the cells above the panel tip is almost uniform. Looking at the volumetric strain and the shear strain fields Fig. 4(b) and Fig. 4(c), it appears that the near rigid body motion close to the panel skin is preserved, with levels close to zero. The most striking differences with Fig. 2 are, above the panel tip, a significant expansion conjugated with almost no shear. The volumetric strain maximum value, above 30%, is there in agreement with the relative increase of tip gap height. A mild shear strain at the same spot indicates near-isotropic deformation. On the contrary, the maximum shear strain on the right side is beyond 54 deg. Compared to the δ = 15 mm configuration, the vertical elongation there is more severe, increasing the deformation anisotropy. As a conclusion, to counter balance the 17 deg rotation of the panel tip and the 45% widening of the tip gap, the cells above the panel tip are subjected in this configuration to a strong dilatation, quite isotropic, coupled to a slight shear. 3.2.4. Quality metric analysis The size-skew metric histograms for the two configurations are shown on Fig. 5. In both cases, the initial quality is maximum (equal to 1) for all elements of the mesh. Indeed, the test mesh is an uniform cartesian grid with a grid spacing set to 0.5 mm in both directions. After the deformation without slip condition, 80% of the elements still have an excellent quality because the impacted region is confined in space. In the impacted area, the worst quality is about 0.90 for the larger tip gap (Fig. 5(a)) and 0.93 for the smaller one (Fig. 5(b)). The deformation using slip condition slightly increases the overall quality. 3.3. Parameterized variations of fan tip To illustrate 3D applications, the tip gap of a generic aero engine fan is considered (see Fig. 6(a,b)). To analyze in rotation varying conditions, it is subjected to conjugated parametrized variations of tip clearance (up to 1 mm) and stagger angle (in the range ±5 deg). Thereafter, a tip gap reduction of 0.5 mm and a rotation of 2 deg centered on the trailing edge are presented as preliminary results. The initial blade tip is drawn Fig. 6(a) with the control points locations. On Fig. 6(b), the after deformation surface is superimposed on the initial blade tip outline. The largest displacement is smax = 4.14 mm, located at the leading edge. Support radius r0 is set to 12 mm, about 3 times smax . Greedy algorithm tolerance is set to 10−3 smax . Moving control points are located only along the blade top profile. The 3D coordinate system used is such that ( e1 , e2 , e3 ) = ( eR , eRθ , ez ), with ez the rotation axis, eR and eRθ the associated radial and azimuthal directions. In this basis, the cylindrical shroud is planar, and Rθ is well adapted to distance measurement in azimuthal direction compared with R in radial direction. Figure 6(c) shows a meridional
12 360
S.St´ephane Aubert et al. / Procedia 00 (2017) 203 000–000 Aubert et al. / Engineering Procedia Engineering (2017) 349–361
0.2610 0.2605 0.2600
r[m]
0.2595 0.2590 0.2585 0.2580 0.2575 0.2570 −0.10 −0.08 −0.06 −0.04 −0.02 z[m]
(a) Front view with control points
(b) Morphed front view
0.00
0.02
0.04
(c) Meridional view
Fig. 6. Fan tip deformation - (a) Moving points (black) and sliding points (white) location on the blade (grey) and the shroud (alpha grey) before morphing - (b) Blade tip before (black line) and after (grey surface) morphing - (c) Outline of shroud and blade (black dashed lines represent initial positions; red dashed lines represent deformed positions without slip condition; solid black lines represent deformed positions with slip condition).
6 Percentage of cells [%]
0.260
R[m]
0.258 0.256 0.254 0.252 0.250 −0.010
−0.005
0.000 Rθ[m]
0.005
(a) Mesh before (grey) and after (black) morphing
0.010
Initial With slip condition
4
2
0 0.0
0.2
0.4 0.6 Skew metric [-]
0.8
1.0
(b) Skew metric histograms before (grey) and after (black) morphing
Fig. 7. Fan tip deformation - (a) Mid-chord axial view - (b) Quality metric histograms, evaluated for the mesh in the coordinate system ( e x , ey , ey )
view (z, R) of the blade tip, initially at R = 0.258 m, and the shroud at R = 0.260 m (black dashed line). The after deformation outline, without slip condition applied on shroud, is marked by a red dashed line. The trailing edge moves in radial direction of about 0.5 mm only, whereas the leading edge displacement is beyond 1 mm along eR . Indeed, the blade tip rotation centered on the trailing edge in ( e x , ey , ez ) space brings the leading edge closer to the shroud with increasing stagger angles. r0 being 6 times larger than tip clearance, the displacement without slip condition on shroud is significant. Figure 7(a) shows one every four nodes of the mesh cut at mid-chord (z = −0.03 m) with slip condition applied. The initial grid is drawn in grey, while the after deformation one is in black. 303 moving points and 266 sliding points are selected by the greedy algorithm. After morphing, the shroud radius remains constant; the block meshing the tip clearance appears smoothly compressed radially and shifted in azimuth. It is to notice that by using RBF with compact support, the modified region is confined in space. Values above the display threshold set to 5 deg are limited to a disc centered on blade tip, of about 9 mm in radius, i.e. 32 r0 . Figure 7(b) shows the skew metric [15] for the initial mesh (grey) and the deformed one (black) evaluated for the mesh in the coordinate system ( e x , ey , ey ). The size-skew metric was not employed due to the inadequacy of the relative size metric to evaluate the quality of this kind of meshes. Indeed, the mesh is formed by a hexaedral structured grid, and around the blade, cells are elongated and twisted to discretize the shape. However it is not prejudicial for the
St´ephane Aubert et al. / Procedia Engineering 203 (2017) 349–361 S. Aubert et al. / Procedia Engineering 00 (2017) 000–000
361 13
CFD solution. For this same reason, only 64% of cells have a skew quality larger than 0.8 for the initial mesh. The mesh has an acceptable quality considering 95% of the cells have a quality larger than 0.50. The morphing degrades slightly the quality, such as 54% and 89% of the cells have a quality larger than 0.8 and 0.5, respectively. 4. Conclusion A new formulation extending the well-known interpolation scheme based on Radial Basis Functions to morph meshes has been detailed. It introduces, besides the classical moving control points, a subset of sliding control points for which conditions of zeroed normal displacement and zeroed tangential contributions are prescribed in a casespecific vector space. Multiple slip planes are handled together, while the linear property, linking the deformation response to the known displacements, underlying the RBF formalism is preserved. Test cases analysis was based on deformation quality metrics derived from the displacement field gradient, keeping the approach free from connectivity information. It shows that the proposed RBF formulation is well adapted to yield smooth displacements fields whatever the distance between moving and sliding control points is. Further works will focus on 3D configurations for which the radial basis function expression, coupled to the support radius value, are key components to find a compromise between efficiency, robustness and displacement field smoothness. References [1] S. Jakobsson, O. Amoignon, Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization, Computers & Fluids 36 (2007) 1119–1136. [2] M. Biancolini, I. Viola, M. Riotte, Sails trim optimisation using cfd and rbf mesh morphing, Computers & Fluids 93 (2014) 46–60. [3] A. de Boer, M. Van der Schoot, H. Bijl, Mesh deformation based on radial basis function interpolation, Computers & Structures 85 (2007) 784–795. [4] T. Rendall, C. Allen, Unified fluid–structure interpolation and mesh motion using radial basis functions, International Journal for Numerical Methods in Engineering 74 (2008) 1519–1559. [5] Y. Rozenberg, S. Aubert, G. Benefice, Fluid structure interaction problems in turbomachinery using rbf interpolation and greedy algorithm, in: ASME Turbo Expo 2014 Turbine Technical Conference and Exposition, number 26735 in GT2014, ASME, 2014, p. V02BT39A033. [6] J. T. Batina, Unsteady euler airfoil solutions using unstructured dynamic meshes, AIAA journal 28 (1990) 1381–1388. [7] C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional springs for two-dimensional dynamic unstructured fluid meshes, Computer methods in applied mechanics and engineering 163 (1998) 231–245. [8] E. J. Nielsen, W. K. Anderson, Recent improvements in aerodynamic design optimization on unstructured meshes, AIAA journal 40 (2002) 1155–1163. [9] T. W. Sederberg, S. R. Parry, Free-form deformation of solid geometric models, Computer Graphics 20 (1986) 151–160. [10] J. A. Samareh, Multidisciplinary Aerodynamic-Structural Shape Optimization Using Deformation (MASSOUD), Technical Report 4911, AIAA, 2000. [11] B. Mocanu, R. Tapu, E. Tapu, Mesh deformation with hard constraints, in: 2013 International Symposium on Signals, Circuits and Systems, ISSCS, IEEE, 2013, pp. 1–4. [12] A. K. Michler, Aircraft control surface deflection using rbf-based mesh deformation, International Journal for Numerical Methods in Engineering 88 (2011) 986–1007. [13] R. Schaback, H. Wendland, Adaptive greedy techniques for approximate solution of large rbf systems, Numerical Algorithms 24 (2000) 239–254. [14] T. C. Rendall, C. B. Allen, Efficient mesh motion using radial basis functions with data reduction algorithms, Journal of Computational Physics 228 (2009) 6231–6249. [15] P. M. Knupp, Algebraic mesh quality metrics for unstructured initial meshes, Finite Elements in Analysis and Design 39 (2003) 217–241. [16] R. H. Parry, Mohr circles, stress paths and geotechnics, Spon Press, 2004, pp. 34–37. [17] F. P. Beer, E. R. Johnston, J. T. DeWolf, Mechanics of materials, Tata McGraw-Hill, 2004, pp. 470–477. [18] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in computational Mathematics 4 (1995) 389–396. [19] M. D. Buhmann, Radial basis functions, Acta Numerica 2000 9 (2000) 1–38. [20] J. Giordano, G. Jourdan, Y. Burtschell, M. Medale, D. E. Zeitoun, L. Houas, Shock wave impacts on deforming panel, an application of fluid-structure interaction, Shock Waves 14 (2005) 103–110. [21] G. Benefice, Y. Rozenberg, S. Aubert, P. Ferrand, F. Thouverez, A partitioned strong coupling procedure for simulation of shock wave/flexible structure interaction, in: ASME 2014 4th Joint US-European Fluids Engineering Division Summer Meeting, number 21624 in FEDSM2014, ASME, 2014, p. V01BT12A006.