Tahsin Khajah

Xavier Antoine

P A Bordas

Isogeometric finite element analysis of time-harmonic exterior acoustic scattering problems

Keywords: Isogeometric Analysis, time-harmonic acoustics, Scattering, Helmholtz problem, high frequency, pollution error

We present an isogeometric analysis of time-harmonic exterior acoustic problems. The infinite space is truncated by a fictitious boundary and (simple) absorbing boundary conditions are applied. The truncation error is included in the exact solution so that the reported error is an indicator of the performance of the isogeometric analysis, in particular of the related pollution error. Numerical results performed with high-order basis functions (third or fourth orders) showed no visible pollution error even for very high frequencies. This property combined with exact geometrical representation makes isogeometric analysis a very promising platform to solve high-frequency acoustic problems.

Introduction

The purpose of this paper is to investigate the performance of Non-Uniform Rational B-Splines (NURBS)-based Isogeometric Analysis (IGA) finite element methods (IGAFEM) for solving time-harmonic acoustic wave scattering and acoustic wave propagation problems, most particularly for small wavelengths. To fix the notations, let us consider Ω -as a d-dimensional (d = 1, 2, 3) scattering bounded domain of R d , with shape Γ := ∂Ω -. We introduce the associated exterior (i.e. unbounded) domain of propagation Ω + := R d /Ω -. Then, solving the scattering problem leads to computing the wave field u, solution to the Boundary-Value Problem (BVP): given an incident (plane wave) field u inc , find u such that ∆u + k 2 u = 0, in Ω + ,

∂ nΓ u = g := -∂ nΓ u inc , on Γ, lim |x|→+∞ |x| (d-1) ∇u • x |x| -iku = 0, (1)
where ∆ is the Laplacian operator, ∇ the gradient operator and n Γ is the outward-directed unit normal vector to Ω -. The spatial variable is x = (x, y) in 2D and x = (x, y, z) in 3D. The wavenumber k is related to the wavelength λ by the relation: λ := 2π/k. The boundary condition is a Neumann (sound-hard) boundary condition (but other boundary conditions could also be considered).

Denoting by a • b the hermitian inner-product of two complex-valued vector fields a and b, then the last equation of system [START_REF]Error estimation and adaptivity for the finite element method in acoustics: 2d and 3d applications[END_REF] is known as the Sommerfeld's radiation condition [START_REF] Colton | Integral equation methods in scattering theory[END_REF][START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF] at infinity, which presents the outgoing wave to the domain. Various numerical methods have been developed to accurately solve the exterior scattering problem, e.g. finite elements [START_REF] Bouillard | Error estimation and adaptivity for the finite element method in acoustics: 2d and 3d applications[END_REF][START_REF] Briot | Efficient implementation of highorder finite elements for helmholtz problems[END_REF][START_REF] Gabard | International journal for numerical methods in engineering[END_REF][START_REF] Giorgiani | Highorder continuous and discontinuous galerkin methods for wave problems[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF][START_REF] Ihlenburg | Finite element solution of the Helmholtz equation with high wave number .2. The h-p version of the FEM[END_REF][START_REF] Ortiz | An improved partition of unity finite element model for diffraction problems[END_REF][START_REF] Thompson | A review of finite-element methods for time-harmonic acoustics[END_REF]Thompson and Pinsky], boundary elements [START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF][START_REF] Antoine | Generalized Combined Field Integral Equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF][START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF][START_REF] Chandler-Wilde | Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering[END_REF][START_REF] Colton | Integral equation methods in scattering theory[END_REF][START_REF] Darbas | Combining analytic preconditioner and fast multipole method for the 3-D Helmholtz equation[END_REF][START_REF] Liu | Fast Multipole Boundary Element Method[END_REF][START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF], enriched and various modified (wave-based, hybrid, asymptotic) finite elements [START_REF] Antoine | Phase reduction models for improving the accuracy of the finite element solution of time-harmonic scattering problems i: General approach and low-order models[END_REF][START_REF] Astley | Finite element formulations for acoustical radiation[END_REF][START_REF] Cessenat | Application of an ultra weak variational formulation of elliptic pdes to the two-dimensional helmholtz problem[END_REF][START_REF] Farhat | The discontinuous enrichment method[END_REF][START_REF] Gabard | A comparison of wave-based discontinuous Galerkin, ultra-weak and least-square methods for wave problems[END_REF][START_REF] Geuzaine | An amplitude formulation to reduce the pollution error in the finite element solution of time-harmonic scattering problems[END_REF][START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF][START_REF] Harari | Numerical investigations of stabilized finite element computations for acoustics[END_REF][START_REF] Huttunen | Comparison of two wave element methods for the helmholtz problem[END_REF][START_REF] Kechroud | A performance study of plane wave finite element methods with a pad-type artificial boundary condition in acoustic scattering[END_REF][START_REF] Laghrouche | Short wave modelling using special finite elements[END_REF]. A challenging and still outstanding question for numerics and applications is related to the so-called high-frequency regime, where the wavelength λ is very small compared with the characteristic length of the scatterer Ω - [START_REF] Babuska | Is the pollution effect of the fem avoidable for the helmholtz equation considering high wave numbers[END_REF][START_REF] Thompson | A review of finite-element methods for time-harmonic acoustics[END_REF].

A first "exact" method consists in writing an integral equation on the surface Γ to represent the exterior field, and then solving this by means of a boundary element method, combined with fast evaluation algorithms (e.g. Fast Multilevel Multipole Methods (MFMM) [START_REF] Darbas | Combining analytic preconditioner and fast multipole method for the 3-D Helmholtz equation[END_REF][START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Liu | Fast Multipole Boundary Element Method[END_REF], Adaptive Cross Approximation (ACA) [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF][START_REF] Zhao | The adaptive cross approximation algorithm for accelerated method of moments computations of emc problems[END_REF]) and preconditioned (matrix-free) Krylov subspace iterative solvers (GMRES) [START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF][START_REF] Antoine | Generalized Combined Field Integral Equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF][START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. This approach has the very interesting property to reduce the initial problem to a finite (d -1)-dimensional problem (set on Γ) and to be relatively stable thanks to the frequency regime (typically between n λ = 5 to 10 points per wavelength are used, with n λ = λ/h, and h the meshsize). Nevertheless, the method is nontrivial to adapt to complex geometries, when the boundary conditions are modified or when a high-order of accuracy is required for the solution. In particular, developing efficient algorithms when the geometry is described with high accuracy remains an open question, most particularly for large wave numbers k since the number of degrees of freedom is potentially very large.

Another standard and very attractive engineering-type approach consists in truncating the exterior domain Ω + through the introduction of a fictitious boundary/layer Σ that surrounds the scatterer Ω -, resulting in a bounded computational domain Ω b with inner boundary Γ and outer boundary/layer Σ. When Σ is a boundary, an Absorbing Boundary Condition (ABC) [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF][START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF][START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF][START_REF] Hagstrom | New results on absorbing layers and radiation boundary conditions[END_REF][START_REF] Kechroud | Numerical accuracy of a Padé-type non-reflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency[END_REF][START_REF] Medvinsky | On surface radiation conditions for an ellipse[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF][START_REF] Turkel | Boundary conditions and iterative schemes for the Helmholtz equation in unbounded region[END_REF] is set on the fictitious boundary which introduces a truncation error even at the continuous level. This can be reduced by considering high-order (local) ABCs. A very popular alternative is, rather, to consider a surrounding absorbing layer and modifying the nonphysical media to obtain a Perfectly Matched Layer [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bermudez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Bermúdez | Perfectly matched layers for time-harmonic second order elliptic problems[END_REF][START_REF] Chew | A 3D Perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF].

In this paper, since our goal is to understand the quality of the NURBSapproximations used in Isogeometric Analysis (IGA) to approximate the wave field, we will consider simple boundary conditions where we can separate the continuous truncation error from the numerical approximation. High-order ABCs as well as PML will be treated in future work. When the problem is truncated, any numerical method can be used such as e.g. finite-element methods, the resulting linear system being then solved through well-adapted domain decomposition techniques [START_REF] Boubendir | A quasi-optimal nonoverlapping domain decomposition algorithm for the Helmholtz equation[END_REF][START_REF] Després | Décomposition de domaine et problème de Helmholtz[END_REF][START_REF] Després | A domain decomposition method for the harmonic Maxwell equations[END_REF][START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF]. In the high frequency regime, it is well-known that the numerical solution suffers from a phase shift due to numerical dispersion which is called the pollution error [START_REF]Error estimation and adaptivity for the finite element method in acoustics: 2d and 3d applications[END_REF][START_REF] Babuska | Is the pollution effect of the fem avoidable for the helmholtz equation considering high wave numbers[END_REF][START_REF] Ihlenburg | Finite element solution of the helmholtz equation with high wave number part i: The h-version of the fem[END_REF][START_REF] Ihlenburg | Finite element solution of the helmholtz equation with high wave number part ii: The h-p version of the fem[END_REF] of the finite element method. To maintain a prescribed pollution error, the mesh density n λ should be increased faster than the wave number leading to high computational cost for high frequency problems. Therefore, the FEM is limited to a upper frequency bound for which the computational cost becomes prohibitive and increasing the order of the polynomial basis functions used in conventional FEM is required to reduce the pollution error, but it does not fully eliminate it.

Original techniques were recently proposed to reduce the pollution error. One approach is to replace the basis functions, usually a polynomial, with plane wave functions [START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF][START_REF] Laghrouche | Short wave modelling using special finite elements[END_REF]. Since pollution error can be viewed as a phase shift, another hybrid approach is to asymptotically approximate the phase of the solution and to reformulate the problem into an equivalent problem with a slowly varying unknown [5-7, 9, 10, 19, 38, 73]. Nevertheless, until now, the only viable solution for complex problems is to consider a sufficiently high-order polynomial basis into the FEM together with high-order meshes to represent accurately the geometry Γ.

The aim of the present paper is to analyze the quality of alternative approx-imation methods based on B-splines and NURBS-based Isogeometic Analysis (IGA) and to compare the properties of these approximation schemes to more traditional, Lagrange-based high-order FEM. The goal of the paper is, specifically, to investigate the properties of the methods for large wave numbers k. Isogeometric analysis (IGA) was introduced over 10 years ago to streamline the transition from Computer-Aided Design to Analysis [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF][START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]. The central idea of the approach is to use the same shape functions to approximate the field variables as those used to describe the geometry of the domain. As the geometry of computational domains is typically provided by Computer-Aided Design (CAD) software, it is sensible to use NURBS shape functions, which are the most commonly used functions in CAD. Since 2005, IGA has been developed within the finite element, boundary element [START_REF] Lian | Stress analysis without meshing: Isogeometric boundary-element method[END_REF][START_REF] Simpson | An isogeometric boundary element method for elastostatic analysis: 2d implementation aspects[END_REF][START_REF] Simpson | A twodimensional isogeometric boundary element method for elastostatic analysis[END_REF] and was enriched through partition of unity [START_REF] Nguyen | Isogeometric analysis: An overview and computer implementation aspects[END_REF] IGA benefits from exact (and smooth) geometry representation and no loss of boundary details. This is fundamental in problems involving wave propagation or turbulence, as spurious wave scattering and vortices typically appear close to rough boundaries. IGA therefore alleviates the meshing burden as if the CAD model is modified, any modifications are immediately inherited by the shape functions used for analysis. When coupled to boundary element approaches, this suppresses completely the generation and regeneration of the mesh without any regeneration of the mesh [START_REF] Lian | Implementation of regularized isogeometric boundary element methods for gradient-based shape optimization in two-dimensional linear elasticity[END_REF][START_REF] Lian | Stress analysis without meshing: Isogeometric boundary-element method[END_REF]. When used within a finite element approach, volume parameterization is still required, but the shape optimization process is simplified [START_REF] Fusseder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF][START_REF] Wall | Isogeometric structural shape optimization[END_REF]. Since CAD models are generated from the boundary data, IGA Boundary Element Method (IGABEM) is promising in bridging the gap between CAD and analysis offering a unique shape optimization package [START_REF] Coox | An isogeometric indirect boundary element method for solving acoustic problems in open-boundary domains[END_REF][START_REF] Kostas | Ship-hull shape optimization with a t-spline based bemisogeometric solver[END_REF][START_REF] Lian | Sensitivity analysis and shape optimization through a t-spline isogeometric boundary element method[END_REF][START_REF] Scott | Isogeometric boundary element analysis using unstructured tsplines[END_REF].

IGA additionally offers higher inter-patch continuity, large support sizes and superior refinement capabilities compared to Lagrange-based approximation schemes. A recent review of IGA and its computer implementation aspects is available in [START_REF] Nguyen | Isogeometric analysis: An overview and computer implementation aspects[END_REF]. GEOPDE is another open source IGA research tool [START_REF] De Falco | Geopdes: A research tool for isogeometric analysis of {PDEs}[END_REF] which provides a flexible coding structure. The application of IGA and its performance in solving one-and two-dimensional eigenvalue and interior Helmholtz problems have been studied [START_REF] Coox | A performance study of nurbs-based isogeometric analysis for interior two-dimensional timeharmonic acoustics[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF][START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method {NURBS}[END_REF][START_REF] Hughes | Finite element and {NURBS} approximations of eigenvalue, boundary-value, and initial-value problems[END_REF] in which IGA showed superiority to conventional FEM. However, the performance of IGA in solving exterior scattering problems and related pollution error is still not thoroughly investigated.

In this paper, we study the performance of IGA for solving time-harmonic exterior scattering problems in one, two and three dimensions. First, in Section 2, we briefly present the Isogeometric Analysis and the related NURBS and B-Spline basis functions. Then in Section 3, we consider a simple but meaningful one-dimensional scattering problem and provide a pseudo-code to solve it in IGA context. In Section 4, we consider two-dimensional problems and study the performance of IGA in estimating the scattered field of a duct and a circular cylinder. The in-house IGA and Lagrange-based FEM codes developed in this study were written in Matlab R following the GEOPDE platform for two-and three-dimensional examples. Finally, section 5 conclude the paper.

Isogeometric Analysis (IGA)

In IGA, both the physical model and the solution space are constructed by B-spline/NURBS functions. The model is constructed using a set of control points and knot vectors in each spatial direction. A knot vector is defined as a set of non-decreasing coordinates and represented by ξ = {ξ 1 , ξ 2 , ..., ξ (n+p+1) }, where ξ i is the i th knot and i is the knot index, i = 1, 2, • • • , n + p + 1, where p is the polynomial order and n is the number of basis functions. Usually an open knot vector is used in IGA where the first and the last knot values appear p + 1 times. The elements in IGA are defined as non-zero knot intervals. If the knots are spaced equally within the knot vector, the knot vector is called uniform otherwise it is non-uniform.

B-Splines

A B-Spline basis function is defined recursively by the Cox-de Boor recursion formula starting with the zeroth order (p = 0) basis function:

for p = 0, N i,p (ξ) = 1 ξ i ≤ ξ ≤ ξ i+1 , 0 otherwise, (2)
for p = 1, 2, 3, • • • N i,p (ξ) = ξ -ξ i ξ i+p -ξ i N i,p-1 (ξ)+ ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 N i+1,p-1 (ξ), (3)
where 0/0 is defined to be zero. To compute the tangent and the normal functions, it is required to calculate the derivative of B-Splines. Efficient algorithms to calculate B-Splines and their derivatives can be found in the NURBS Book by Piegl and Tiller [START_REF] Piegl | Introduction to finite element analysis : formulation, verification and validation[END_REF]. The first-order B-Spline functions are identical to their Lagrangian (FEM) counterparts. Also, the B-Splines number of required shape functions for a specific order is similar to that of the FEM. For example, three shape functions are required to construct second-order Lagrangian or B-Spline shape functions. In contrast to FEM, the shape functions in IGA are nonnegative. The first derivative of the B-spline basis function is calculated using the following recursive formula:

d dξ N i,p (ξ) = p ξ i+p -ξ i N i,p-1 (ξ) - p ξ i+p+1 -ξ i+1 N i+1,p-1 (ξ). (4)
These functions exhibit C p-mi continuity across knot ξ i where m i is the multiplicity of the knot ξ i . The B-Spline basis function corresponding to p = 1, • • • , p = 5 are shown in Fig. 1 where the knot vectors are defined to represent one element (knot span)

ξ = { 0, • • • , 0 p+1 times , • • • , 1, • • • , 1 p+1 times }. (5)
The quadratic B-Spline basis function for a 1D domain divided into four ele- ments is shown in Fig. 2 where the knot vector is defined as

ξ = {0, 0, 0, 1/4, 1/2, 3/4, 1, 1, 1}.
The four non-zero knot spans correspond to the four elements of this onedimensional domain.

B-Spline curves

The piecewise-polynomial B-Spline curves in R d are constructed as a linear combination of B-Spline basis functions:

C(ξ) = n i=1 N i,p (ξ)B i . (6)
The vector-valued coefficients of the basis functions in Eq.6 are control points

(B i ∈ R d , i = 1, 2 • • • , n). Since C(ξ)
is representing a curve, the corresponding control points, B i , are analogous to nodal coordinates in FEM. The resulting curve is C p-1 -continuous everywhere except at the location of the repeated knots with m i multiplicity, where it is C p-mi -continuous. Therefore, when utilizing open knot vectors, the resulting curve is C 0 -continuous, tangent to the control polygon, and interpolatory at the start and the end of the curve and where the multiplicity is equal to the polynomial order. The support of a B-Spline basis function is local. In other words, moving a single control point can affect the geometry of no more than p + 1 elements of the curve.

B-Spline surfaces and volumes

A tensor product B-Spline surface is defined as:

S(ξ, η) = n i=1 m j=1 N i,p (ξ)M j,q (η)B i,j , (7)
where B i,j , i = 1, 2, • • • , n, j = 1, 2, • • • , m are control points, and N i,p (ξ) and M j,q (η) are univariate B-Spline functions of orders p and q corresponding to the knot vectors

Ξ 1 = {ξ 1 , ξ 2 , • • • , ξ n+p+1 } and Ξ 2 = {η 1 , η 2 , • • • , η m+q+1 },
respectively. Similarly, a B-Spline volume is constructed as a tensor product of three B-Spline basis functions:

V(ξ, η, ζ) = n i=1 m j=1 l k=1 N i,p (ξ)M j,q (η)L k,r (ζ)B i,j,k . (8)

NURBS

There are certain geometries, such as circles and ellipsoids, which B-Spline functions and Lagrange polynomials cannot represent exactly. In the CAD community, B-Spline functions are replaced with a more general form to overcome this shortcoming. This new form is called Non-Uniform Rational B-Splines (NURBS). NURBS shape functions are defined as:

R p i (ξ) = N i,p (ξ)w i W (ξ) = N i,p (ξ)w i n i=1 N i,p (ξ)w i , (9)
where {N i,p } is a set of B-Spline basis functions and {w i } is a set of positive NURBS weights. By appropriate weight selection both polynomials and circular arcs can be described. NURBS curves are defined as follows:

C(ξ) = n i=1 R p i (ξ)B i . (10)
Similarly, rational surfaces and volumes are defined in terms of rational basis functions as:

R p,q i,j (ξ, η) = N i,p (ξ)M j,q (η)w i,j n i=1 m j=1 N i,p (ξ)M j,q (η)w i,j , (11)
R p,q i,j (ξ, η, ζ) = N i,p (ξ)M j,q (η)L k,r (ζ)w i,j,k n i=1 m j=1 l k=1 N i,p (ξ)M j,q (η)L k,r (ζ)w i,j,k . ( 12
)
If the weights are all equal, NURBS basis functions will reduce to their B-Spline counterparts and the corresponding curve becomes a non-rational polynomial again. The geometry of an element in IGA can be expressed as:

x e ( ξ) = nen a=1 P e a R e a ( ξ), (13)
where n en = (p+1) dp and d p is the spatial dimension. The field u(x) is similarly expressed as:

u( ξ) = nen a=1 d e a R e a ( ξ), (14)
where d e a R e a is the set of control variables. As a consequence of the parametric definition of basis functions in IGA, it is necessary to consider two mapping for integration

Ω f (x)dΩ = n el e=1 Ω e f (x)dΩ = n el e=1 e Ω f (x(ξ))|J ξ |d Ω = n el e=1 Ω f (x( φe ( ξ))|J ξ ||J ξ|d Ω, (15)
where Ω and Ω are representing the parent and parameter spaces respectively, and Ω is the physical space.

A one-dimensional toy scattering problem

Now, let us come back to our physical problem by considering the simple one-dimensional Neumann scattering problem solved in [START_REF] Antoine | Phase reduction models for improving the accuracy of the finite element solution of time-harmonic scattering problems i: General approach and low-order models[END_REF]:

∂ 2 x u + k 2 u = 0, in Ω b = (0, 1), ∂ x u = ik, at Γ = {0}, ∂ x u -iku = 0, at Σ = {1}. (16)
The variational formulation is given as:

Ω b {∂ x u∂ x v -k 2 uv}dx -ik(uv)(1) = -(gv)(0), ∀v ∈ H 1 (Ω b ). ( 17
)
The solution of the one-dimensional BVP ( 16) is: u ex (x) = e ikx , which represents the scattering of an incident plane wave by the left half space. An exact transparent boundary condition is considered for the fictitious boundary at Σ = {1} using the Dirichlet-to-Neumann operator Λ = ik on Σ, resulting in the transparent boundary condition

∂ x u = Λu.
It is shown that the Finite Element approximated solution of the BVP described in system ( 16) suffers from the pollution error [START_REF] Antoine | Phase reduction models for improving the accuracy of the finite element solution of time-harmonic scattering problems i: General approach and low-order models[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF][START_REF] Ihlenburg | Finite element solution of the helmholtz equation with high wave number part ii: The h-p version of the fem[END_REF][START_REF] Thompson | A review of finite-element methods for time-harmonic acoustics[END_REF]. The pollution error has a direct relation with the wavenumber k and increases with frequency. It is the polynomial basis used in conventional FEM which is inadequate to represent the wavefield. Since the main difference of the IGA with the conventional FEM is in the selection of the basis function, it is interesting to study the performance of IGA specially in high frequencies and its pollution error. We divide the computational bounded domain Ω b into N h uniform nonzero knot spans (elements) so that each segment length is 1/h. We define the density of discretization as the number of elements per wavelength and denote it with n λ = λ/h. A Matlab code was prepared to solve the above example in IGA context. The stiffness matrix for a single element is calculated and assembled into the global matrix following these steps:

• Set K = 0
• loop over number of elements • Loop over Gaussian points ξi , wi 1. find the parametric coordinate ξ = φe ( ξj ) 2. compute the basis function R e a (a = 1, ...p + 1) at point ξ 3. define vector R = [R e 1 , R e 2 , ..., R e p+1 ] 4. compute the basis function derivatives R e a,ξ (a = 1, ...p + 1) at point ξ 5. define vector

R ξ = [R e 1,ξ , R e 2,ξ , ..., R e p+1,ξ ] 6. compute |J ξ | = R ξ P 7. compute |J ξ | = 0.5(ξ i+1 -ξ i ) 8. compute shape function derivatives R x = J -1 ξ R T ξ 9. Set k e = k e + R x R T x |J ξ ||J ξ | wj -k 2 RR T |J ξ ||J ξ | wj 10.
End of both loops then, the Neumann boundary condition is applied and the resulting linear system is solved to find the estimated solution in IGA. we plot the approximate solution u h for k = 40 as well as the exact solution in Fig. 3. The density of discretization for IGA approximations is n λ = 10. As expected, IGA (p = 1) also suffers from pollution error. B-Spline functions generated with p = 0 and Figure 3: The pollution error is the numerical solutions (u h ) of IGA (p = 1), is visible when compared with the exact solution (u ex ), for k = 40 and n λ = 10. p = 1 orders will result in the same piecewise constant and linear functions as standard FEM shape functions. The differences are found for cases using higher order shape functions. The p th order B-Spline function has p -1 continuous derivatives across the element boundaries. Each B-Spline basis is point-wise non-negative over the entire domain as a result, all of the entries of the mass matrix will be positive. The support of a B-Spline function of order p is always p + 1 knot spans. Therefore, higher order B-Spline function has support over much larger portion of the domain when compared to classical FEM. The total number of functions that any given function shares support with (including itself) is 2p + 1 regardless of whether a FEM basis or a B-Splines is used. Hence, using higher order B-Spline basis functions will provide improved support compared to standard FEM without increasing the number of required shape functions. Increasing the order of Lagrangian polynomials will increase the amplitude of oscillations in FEM; this problem is eliminated in IGA as a result of non-negativity and non-interpolatory nature of B-Splines shape functions.

Next we investigate the effect of increasing the basis function order on the accuracy of the IGA approximation again for discretization density n λ = 10. As expected increasing the order p reduces the error drastically as shown in Fig. 4. There is no visible pollution error in the IGA approximation if we keep all the parameters unchanged and increase the order of the basis function. One can reduce the error further by increasing the discretization density n λ . For the numerical solution f calc (x), x ∈ Ω b , we define the L 2 -error (in dimension d) as

2 = Ω b |f calc (x) -f ex (x)| 2 dx 1/2 . ( 18
)
The evolution of the L 2 -error vs. the discretization density n λ is reported in Fig. 5 for the wavenumber k = 40. For completeness, we also plot the curves relative to the function x -p , for x = n λ . We can see that the error curves fits very well the polynomial curves. This lets to think that the pollution error is extremely weak, most particularly when p is large enough (p ≥ 3) even for a small density n λ . Next we present the evolution of the L 2 -error with the wavenumber k in Fig. 6 for a fixed density n λ = 20, where for p = 3 and higher the pollution error is not visible and the error remains constant even for extremely high frequencies.

Figure 6: Evolution of the error 2 vs. the wavenumber k, for n λ = 20 and various approximation orders p = 1, ..., 5.

Two-dimensional examples

In this section, we evaluate the performance of IGA in solving two-dimensional acoustic problems. First, we consider a two-dimensional duct problem with rigid walls [START_REF] Huttunen | Comparison of two wave element methods for the helmholtz problem[END_REF]. Next, we propose two examples of cylindrical disk scattering problems where we obtain, first the scattered field of a particular mode, and then the scattering of the disk subject to an incident plane wave. By separating the truncation error from the numerical basis approximation, we fairly investigate the performance of IGA in solving the two-dimensional sound-hard exterior scattering problem and numerically analyze the related pollution error.

The duct problem

The duct model is shown in Fig. 7 where Ω b = [0, 2] × [0, 1]. We denote the outward boundary unit normal with v and assume that the lower and upper walls are rigid. We solve the Helmholtz problem stated in Eq. 19 for the acoustic pressure u ∆u + k 2 u = 0, in Ω b , ∂u ∂v = cos(mπy), on x = 0,

∂u ∂v + iku = 0, on x = 2, ∂u ∂v = 0, on y = 0, 1, (19)
where m ∈ N is the mode number. An inhomogeneous boundary condition is applied on the inlet boundary (x = 0) and an absorbing (and transparent for m = 0) boundary condition is set on the outlet boundary (x = 2). Since the boundaries at y = 0, 1 are assumed to be perfectly rigid, the normal derivative of the acoustic pressure vanishes on these boundaries. The exact solution of problem [START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF] with ABC is as follows

u ex (x, y) = cos(mπy)(A 1 e -ikxx + A 1 e ikxx ), (20)
where k x = k 2 -(mπ) 2 and the coefficients A 1 and A 2 are obtained from

k x -k x (k -k x )e -2ikx (k + k x )e 2ikx A 1 A 2 = 1 0 .
Solving this 2 × 2 linear system, we are able to get the expression of the solution to the duct problem with ABC. This is a very interesting point since we therefore can only focus the analysis on the finite-dimensional approximation without including the (continuous) truncation error. The cut-off frequency is m cut-off = k/π. If the mode is such that m ≤ m cut-off , the solution is representing propagating modes and if m > m cut-off the solution corresponds to evanescent modes. The real parts of the exact and estimated IGA solutions are presented in Fig. 8a and8b, respectively, for k = 40, m = 2 (propagative mode), p = 3, and n λ = 10. The corresponding absolute error |u ex -u h | is plotted in Fig. 9, the maximal value being of the order of 10 -6 . For the propagative mode m = 2, the evolution of the L 2 -error with respect to the discretization density n λ is shown in Fig. 10 for k = 40, and according to the wave number k in Fig. 11, where n λ = 10. We can see that the error is very low, decreasing strongly with p. In addition, the error does not seem to depend on k as observed in Fig. 11, meaning that the pollution error is negligible.

(a) (u ex ) (b) (u h )

The sound-hard circular cylinder problem

To further study the performance of IGA, we analyze the scattering by a 2D circular cylinder Ω -= D 0 of radius R 0 centered at the origin, with boundary Γ = C 0 (circle of radius R 0 ). The scatterer is surrounded by an outer fictitious circular boundary Σ = C 1 , again centered at the origin and with radius R 1 > R 0 . Hence, the computational domain Ω b is the annulus bounded between the inner C 0 and outer C 1 boundaries. On the exterior boundary, we set the symmetrical second-order Bayliss-Turkel ABC [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF] given by

∂ nΣ u + (-ik + κ 2 - κ 2 8(κ -ik) )u - 1 2(κ -ik) ∂ 2 s u = 0, on Σ, (21)
setting n Σ at the outwardly directed unit normal to Σ,

∂ nΣ := ∂ r the normal derivative, κ = 1/R 1 is the curvature of Σ := C 1 and ∂ 2 s := R -2 1 ∂ 2 φ
is the second-order curvilinear derivative on C 1 . In the above expressions, the polar coordinate system is denoted by (r, φ).

Scattering by a sound-hard circular cylinder : mode-by-mode analysis

We consider an incident wave with a fixed mode m which is described as

u inc m (x) = J m (kr)e imφ , m ∈ Z, (22)
where J m is the m-th order Bessel's function. The exact exterior modal solution u ex m to the truncated scattering problem of the inner sound-hard circular cylinder C 0 in polar coordinate (r, φ) is given as

u ex m (x) = a m H (1) m (kr) + b m H (2) m (kr) e imφ , r ≥ R 0 , m ∈ Z. (23)
The functions H

m ) is the first-kind (respectively secondkind) Hankel function of order m. The Neumann boundary condition is applied on C 0 and the ABC on C 1 , resulting in the following linear system of equations to obtain the two unknown coefficients a m and b m

a m = - A m 22 J m (kR 0 ) D m , b m = - A m 21 J m (kR 0 ) D m , (24)
where

D m = A m 11 A m 22 -A m 21 A m 12 and                A m 11 = H (1)
m (kR 0 ), A m 21 = kH (1)
m (kR 1 ) -B m H (1)
m (kR 1 ),

A m 12 = H (2)
m (kR 0 ),

A m 22 = kH (2)
m (kR 1 ) -B m H (2)
m (kR 1 ),

B m = -α m m 2 R 2 1 + β m , α m = - 1 2ik 1 + i kR 1 -1 , β m = -ik + 1 2R 1 + 1 8iR 1 (1 + kR 1 )
.

The notation f := ∂ r f designates the radial derivative of a given function f (r). Quite similarly to the duct problem, the modes m ∈ N such that |m| < kR 0 are propagating. For |m| > kR 0 , the modes are evanescent (and therefore are not visible in the far-field since they do not propagate). A special case corresponds to |m| ≈ kR 0 , which is related to grazing modes that are tangent to the scatterer (generating glancing rays). The modal analysis is meaningful since we can understand the accuracy of the IGA approximations thanks to the spatial frequencies m and the dimensionless wavenumber kR 0 . In addition, this also helps in clarifying what can be expected for the case of the scattering of a plane wave by the disk since the exact solution is built as a modal series expansion of the elementary mode solutions. We fix R 0 = 1 and R 1 = 2. The IGA model was generated using four identical patches as shown in Fig. 12. We report the real parts of u ex m and u h on Fig. 13 where the wavenumber is k = 40 and m = 2 (which is a propagative mode since m = 2 ≤ kR 0 = 40). For IGA, the order of the method is p = 3 for a density of 16, for m = 2, with various approximation orders p of IGA and n λ = 10. We can observe that the accuracy does not really depend on k for a fixed density when p ≥ 2. In fact, one question is to know how this error depends on m for a We can remark that the error almost does not depend on m, the value m = 2 corresponding more or less to the worst situation. When |m| ≥ kR 0 , the error starts decreasing strongly. Indeed, the solution is then evanescent and rather corresponds to the solution of a positive definite problem (similar to a shifted laplacian). Such solutions are easier to compute with high accurary than for a propagative solution and do not participate to the far-field. From these remarks, and since the solution of the scattering of a plane wave by a disk is computed as the superposition of these modes solutions, we can expect that similar graphs to Fig. 16 could be obtained for the full plane wave problem (see subsection 4.2.2).

Scattering by a sound-hard circular cylinder: plane wave scattering

In the previous section, the scattering problem of a sound-hard circular cylinder was found for a given harmonics m. The solution for the full plane wave incidence can be constructed by superposition as a Fourier series expansion of these harmonics. In this section, we consider an incident plane wave u inc (x) = e ikd•x , where d is the incidence direction d = (cos(θ inc ), sin(θ inc )) T and θ inc is the scattering angle. Because of the symmetry of the problem, we fix the incidence direction to d = (1, 0) T and the scatterer as the unit disk with R 0 = 1. Similar to the previous example, the second-order Bayliss-Turkel ABC ( 21) is placed on the circle with radius R 1 = 2. We consider the following exact solution [START_REF] Kechroud | A performance study of plane wave finite element methods with a pad-type artificial boundary condition in acoustic scattering[END_REF] to analyze the pollution and approximation errors and to avoid the domain truncation error

u ex = m∈Z (i) m u ex m , (25)
where u ex m is given by Eq. ( 23). To get an accurate reference solution, we truncate the above series expansion by summing up on m from -m max to m max := [kR 0 ] + 30 (where [r] denotes the integer part of a real-valued positive number r). We report on Fig. 18 the real parts of u ex (see Fig. 18a) and u h (see Fig. 18b), where the wavenumber is k = 40. For IGA, the order of the method is p = 3 for a density of discretization points per wavelength n λ = 10. We clearly see that the wavefield is accurately computed since the absolute error |u ex m -u h | is below 10 -4 as it can be deduced from Fig. 19a. The evolution of the error with respect to the discretization density n λ for p = 1, • • • , 5, is shown in Fig. 20 where the wavenumber k is equal to 40. This shows that if p is larger than 3, then the 2 -error starts being small (e.g. less than 10 -2 ) even for small densities n λ (typically n λ ≥ 6). In particular, the error gap between p = 2 and p = 3 is important. In addition, for n λ = 5, we can see on Fig. 19b that the error remains acceptable but of the order of 2 -3%, where the maximal error is located where the amplitude of the wavefield The evolution of the error 2 according to the wavenumber k is depicted in Fig. 21. For this test, we fix n λ = 10 and p = 3 for IGA. As we can see, the error almost does not depend on k, which means that the pollution error is very small. Again, a very good accuracy is obtained for p ≥ 3, the p-th order of the method being then visible for p larger than 3. The dependence with respect to k seems more present for the lowest orders IGA approximations, i.e. p = 1 and p = 2. We consider now a much higher wavenumber, i.e. k = 200. We report on Fig. 22b the real part of the numerical IGA solution for p = 3 and the low density of discretization points per wavelength n λ = 5. The solution can be compared with the exact solution available on Fig. 22a. We immediately see that the two solutions are the same. This is confirmed on Fig. 22c where we report the absolute error |u ex -u h | between the two solutions. This error is of the order of 10 -3 , which is very small. Finally, we report now the results for a very high wavenumber k = 500 on Fig. 23. We provide the real parts of the exact solution u ex (Fig. 23a) and IGA numerical solution u h (Fig. 23b). For IGA, the order of the method is p = 3 for a density of discretization points per wavelength n λ = 5. Even if there is no visible difference, plotting the absolute error on Fig. 23c shows an order of about 6%. This can be made much smaller by increasing the order of approximation to p = 4 by keeping n λ = 5 as seen on Fig. 23d, leading to an error of about 10 -3 .

Conclusions

The exterior Helmholtz problem was solved for one-, and two-dimensional problems using IGA. An absorbing boundary condition is applied on a fictitious boundary to truncate the infinite space. This truncation introduces an error in the numerical solution. To better study the performance of IGA, the truncation error was included in the analytical solution so that the pollution error is distinguishable. It is shown that IGA is less prone to pollution error even for large wavenumbers. Also, the pollution error was considerably reduced by increasing the order of the basis functions which can be conveniently done in IGA, even for small densities of discretization points per wavelenght. The smooth representation of the boundaries in IGA and the low pollution error makes IGA a suitable platform for scattering analysis, in particular for prospecting high frequency problems. For high wave numbers it will be interesting to investigate the use of different spaces for geometry and field approximations. Low order geometry description can be used with high order field approximations to investigate super and sub-geometric approaches, also known as GIFT [START_REF] Xu | Geometry-independent field approximation for spline-based finite element methods[END_REF].

(a) p = 2 (b) p = 3 ( 5 Figure 1 :

 2351 Figure 1: The B-Spline basis functions, p = 1 • • • , p = 5.

Figure 2 :

 2 Figure 2: The quadratic (p = 2) basis function of a 1D domain defined over ξ = {0, 0, 0, 1/4, 1/2, 3/4, 1, 1, 1} and representing four elements.

Figure 4 :

 4 Figure 4: Absolute error |u ex -u h | for k = 40, n λ = 10 and various approximation orders p = 1, ..., 5.

Figure 5 :

 5 Figure 5: Evolution of the error 2 vs. the discretization density n λ for k = 40 and p = 1, ..., 5 (the slopes of the functions x -p are plotted for comparison).

Figure 7 :

 7 Figure 7: The duct domain.

Figure 8 :

 8 Figure 8: Comparing (a) the real parts of the exact solution u ex and (b) the numerical solution u h for k = 40, m = 2, p = 3, and n λ = 10.

Figure 9 :

 9 Figure 9: Absolute error |u ex -u h | for k = 40, m = 2, p = 3 and n λ = 10.

Figure 10 :

 10 Figure 10: Evolution of the error 2 vs. the discretization density n λ , for k = 40, m = 2 and p = 1, ..., 5.

Figure 11 :

 11 Figure 11: Evolution of the error 2 vs. the wavenumber k, for m = 2, n λ = 10 and p = 1, ..., 5.

Figure 12 :

 12 Figure 12: Outline of the patches.

Figure 13 :

 13 Figure 13: Comparing (a) the real parts of the exact solution u ex m and (b) the numerical IGA solution u h for k = 40, m = 2, p = 3, and n λ = 10.

Figure 14 :

 14 Figure 14: Absolute error |u ex m -u h | for k = 40, m = 2, p = 3 and n λ = 10.the effect of increasing the discretization density n λ for various orders p of the underlying IGA basis function on the L 2 -error 2 for k = 40 and m = 2. We see that the error curves follow the slopes of the curves x -p (for x = n λ ). The evolution of the error 2 with respect to the wavenumber k is plotted on Fig.

Figure 15 :

 15 Figure 15: Evolution of the error 2 vs. the discretization density n λ for k = 40, m = 2 and p = 1, ..., 5.

Figure 16 :

 16 Figure 16: Evolution of the error 2 vs. the wavenumber k for m = 2, n λ = 10 and p = 1, ..., 5.

Figure 17 :

 17 Figure 17: Evolution of the error 2 vs. the mode m, for k = 40, n λ = 10 and p = 1, ..., 5.

Figure 18 :

 18 Figure 18: Comparing the real parts of (a) the exact solution u ex and (b) the numerical IGA solution u h for k = 40, p = 3, and n λ = 10.

  (a) n λ = 10 (b) n λ = 5

Figure 19 :

 19 Figure 19: Absolute error |u ex -u h | for k = 40 and p = 3. The discretization density is (a) n λ = 10 and (b) n λ = 5.

Figure 20 :

 20 Figure 20: Evolution of the error 2 vs. the discretization density n λ for k = 40 and p = 1, ..., 5.

Figure 21 :

 21 Figure 21: Evolution of the error 2 vs. the wavenumber k, for n λ = 10 and p = 1, ..., 5.

Figure 22 :

 22 Figure 22: Comparing the real parts of (a) the exact solution u ex and (b) the numerical IGA solution u h , and (c) absolute error |u ex -u h |, for k = 200, p = 3, and n λ = 5.

Figure 23 :

 23 Figure 23: Comparing the real parts of (a) the exact solution u ex and (b) the numerical IGA solution u h for k = 500, p = 3, and n λ = 5. Absolute error |u ex -u h | for (c) p = 3; n λ = 5, (d) p = 4; n λ = 5.


Shuning Li

Hazim El-Mounayri

Weijie Zhang

Bill Schindel
email: schindel@ictt.com

Jason Sherey

Implementation of Systems Engineering Model into Product Lifecycle Management Platform

Keywords: Product Life Cycle Management (PLM), Systems Engineering, Model-Based Systems Engineering, General Manufacturing Model, Systematica Metamodel

Manufacturing companies are facing the challenge of increasing product complexity while at the same time reducing cost and time in a highly competitive global market. Product Lifecycle Management (PLM) and Systems Engineering have the potential to provide solutions for these challenges. The two concepts not only share many common characteristics, but also complement each other. Even though systems engineering and PLM have become closely related in the past few years, implementation of systems engineering models into a PLM platform has rarely been conducted. In this study, the key portion of a model-based system engineering model was implemented into a PLM platform, and the implementation was validated. The results shows that the current implementation can help capture and reflect stakeholders' requirements and changes in product design process promptly and accurately; and reduce the time and potential mistakes in creating a new systems engineering model.

Introduction

In today's world, companies are facing the challenges of increasing product complexity while at the same time reducing cost and time in a highly competitive global market. A direct consequence of the increasing product complexity is a more complex manufacturing system. Product Lifecycle Management (PLM) and Systems Engineering have the potential to help companies manage the more and more complex manufacturing system, and avoid costly product development and launching, as well as reducing product failure. These two concepts not only share many common characteristics [START_REF] Messaadia | PLM as linkage process in a systems engineering framework[END_REF], but also complement each other. On one hand, PLM requires management of the entire product process; it must meet the challenge of synchronizing disciplines involved in complex product systems during the production process. Systems engineering methodologies provide ways to synchronize disciplines during design, simulation, testing, verification, and validation based on multidisciplinary functions for an industrial company [START_REF] Messaadia | PLM as linkage process in a systems engineering framework[END_REF]. On the other hand, PLM can provide a reliable data source to system engineers. This is especially helpful for those fast changing data like requirements. PLM also connected all the information within the processes together, which provides system engineers a tool to trace changes and preform analysis on the affections of the changes. Researchers started trying to integrate PLM and systems engineering to overcome the weaknesses of traditional systems engineering tools and achieve the benefits of using PLM platform [START_REF] Messaadia | PLM as linkage process in a systems engineering framework[END_REF] [START_REF] Bajaj | Slim: collaborative model-based systems engineering workspace for nextgeneration complex systems[END_REF]. But the current integrations either only used PLM as a data source or implemented a traditional systems engineering model. Using PLM only as a data source cannot get the full benefits of a PLM platform, and systems engineers and their works were still isolated from the rest of lifecycle processes and information. Traditional systems engineering models were document centric models which were not easy to implement into a PLM system.

The objectives of this study were to implement the key portion of a model-based system engineering model into a PLM platform, and validate the implementation by generating a product specified model from the general model.

Methodology

The systems engineering model used in this study is a general manufacturing model. This model was built by ICTT System Sciences using Systematica Metamodel (S* Metamodel) methodology [START_REF] Schindel | What is the Smallest Model of a System?[END_REF]. The Fig. 1 depict key element of the S* Metamodel. In this study, the trace of Feature, Functional Interaction, Functional Role, and Design Component was being focused. Features were packages of behavior or performance of a system that have stakeholder value [START_REF] Schindel | What is the Smallest Model of a System?[END_REF]. S* Models were aimed at covering all the stakeholders not only just users or customers. Feature attributes were features' parameters that expressed stakeholder valuations in stakeholder language [START_REF] Schindel | What is the Smallest Model of a System?[END_REF]. For example, the Cruise Control Feature for a car had feature attributes in fuel economy and speed variation. Because Features and Feature Attributes covered all stakeholders' value and interest, they affected all the design decisions, trade-offs, and optimization should be made in accordance. In a system, there always existed interactions between components; an interaction means exchange of energy, force, mass or information [START_REF] Schindel | System interactions: Making the heart of systems more visible[END_REF] which lead to change of state. The Systematica called that interaction a functional interaction, and each component played a functional role in that interaction [START_REF] Schindel | System interactions: Making the heart of systems more visible[END_REF]. The functional role was also called logical system in S* Metamodel. Functional Roles were described by their behavior, and the role attributes were parameters of functional role which had technical valuations [START_REF] Schindel | Pattern-based systems engineering: An extension of modelbased SE[END_REF].

Fig. 1. S* Metamodel [3]

The PLM platform used in this study is Teamcenter ® 10.1. Teamcenter is the most widely used PLM software system in the world [START_REF] Gecevska | Product lifecycle management tools[END_REF]. It helped companies deliver complex products to the market by connecting people with products and process in order to enhance productivity and integrate global operations. Teamcenter ® was the first PLM solution to integrate systems engineering within an entire product lifecycle. It provided a close loop systems engineering environment. The systems engineering environment employed systems engineering methodology to allow an engineer to establish systems requirement, then defined and validated all components and subsystems in the contact of the entire system's lifecycle. The benefit was that products meeted customers' value satisfaction and understood the entire impact of design decisions in the early stages of the lifecycle.

The research approach in this study consisted of three steps: mapping; implementation; and specialization. The first step was to map the S* Metamodel elements into Teamcenter ® classes including relationships. The second step was to build the model into Teamcenter ® . The last step was to generate oil filter end seal compression manufacturing model by specializing the general model.

Mapping

The mapping step transferred systems engineering model elements into PLM classes. The S* model elements, like feature, interface, functional interaction, were mapped Teamcenter ® classes. All these new classes were children to Logical Block class. The relationship between two S* model elements were mapped to either structural relationships or trace link relationships. Trace link relationship was used to represent relationships in hierarchies. For example, features in systems engineering model were created based on stakeholders needs. So the relationship between a need (requirement) and a feature was a trace link type of relationship (A "Need" object is the source of a "Feature" object). Trace links provide traceability between structure elements, and traceability defined as one object was precedent than another object. Most of the relationships were structural relationships which can be presented in a tree view, like the bill of material (BOM) view of a product structure. All of these classes were built by using Teamcenter ® Business Modeler Integrated Devolvement Environment (BMIDE)

Implementation

The implementation step built the general systems engineering model into the PLM platform. A set of features, interactions, and other S* metamodel elements in the general manufacturing model was created in Teamcenter ® . The relationships between feature and feature attribute, feature and interaction, interaction and logical system were created, and presented in a tree view. One trace link relationship between a super class and its children was also created.

Specialization

The oil filter end seal production process was selected to demonstrate the specialization step in this study. Systems engineers usually started with identifying interactions when they built the model. The interaction in this specific model should be specialized from "Transform Material" functional interaction in the general model. A new interaction called "Perform Compression Bonding" was created for the specialized model based on the "Transform Material" functional interaction in the general model.

Once the interaction was decided, the logical systems were modeled. Filter media, bonding compound, and end cap are necessary parts for production, so local airspace, manufacturing system and other manufacturing logical systems were inherited from general production pattern directly. Functional role attributes were created and attached to logical systems, those role attributes were considered as technical valuation of production logical systems.

Based on stakeholders needs, production and manufacturing features were created. Material transformation capability, production capability features and other general production features for manufacturing were able to be populated from the general model directly. And a specialized oil filter production feature that contains engine lubricant filtration feature and reliability (production) features were created, and added to the model. Feature attribute based on stakeholders valuations under each feature were also created.

The Physical Systems were not created in the general manufacturing model because it is difficult to summarize physical system into a general model. But physical systems were necessary for specialized model. In this oil filter end seal case, the physical system included: end seal adhesive, accordion filtration component, and filter cap component. Physical system attributes were created based on the understanding of each physical part and the related physical laws. The CAD models of these physical parts can be linked to physical system elements in the model as references and for change effects analysis.

Results

General Manufacturing Model

The general manufacturing model included groups of feature (feature framework), interactions (interaction framework), logical systems (system environment), and other S* model elements. Fig. 2 showed an example of the functional interaction framework in tree view. The structural relationships were also created, and presented in a tree view. Fig. 3 was an example of the relationships between features and interactions. Microsoft ® Office tools were commonly used by systems engineers, and the reports or diagrams generated in Office tools were still the major deliveries for systems engineers. These reports and diagrams can be generated automatically after proper configuration in Teamcenter ® . Fig. 4 gave an example of an automatically generated diagram based on the implemented in Visio.

Oil Filter End Seal Bonding Compression Manufacturing Model

Corresponding to the implemented general model, the groups of S* metamodel elements for this specific manufacturing system were created during the specialization step. Fig. 5 showed an example of the specialized functional interaction framework. Fig. 6 was an example of a specialized relationship (relationships between interaction and logical Systems) in the tree view. Fig. 7 gave an example of the automatically generated model in Visio. The current implementation only covered the key portion in the general manufacturing model. This portion of the general model had the core elements and their relationships. The implementation of the portion provided a foundation for evaluation the implementation methodology and the implemented model. It also can help present the value of implementing systems engineering models into PLM platform. The rest of the model will be implemented in the future.

One of the most difficult parts of this study is to model structural relationships. There were so much different structural relationships in the model, and these relationships were normally n: n relationships. The whole model was a basically complex network. The current solution to provide a clear presentation of these relationships was to use the tree view. This solution was working well for the general model, but potential problems existed when the model was built for a higher complexity manufacturing process.

The specialization in this study was done manually. However, an automatic specialization should be very beneficial. There were two possible ways to automate the specialization step, the first one was to develop a new tool to generate product specific models from the general model based on systems engineering's inputs; and the second one was to consider the model as a product structure and use product configuration functions in the PLM platform to generate new models. Further studies will be needed to identify a better solution.

Conclusions

A model-based systems engineering model can be implemented into a PLM platform.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Functional Interaction Framework in Tree View

Fig. 4 .

 4 Fig. 4. Interaction Overview Diagram in Visio

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Specialized Functional Interaction Framework


Pradorn Sureephong

Apitchaka Singjai

SaaS for Education: A Case Study of Google Apps in Software Engineering Class

Keywords: Cloud Computing, Google apps, Learning Cloud Services, SaaS, Software Selection

Cloud computing technology in higher education brings cost efficiency and flexibility into the organizations: software as a service (SaaS) solutions require low infrastructural investments and migrate IT resources to the Internet. Adapting to the rapidly changing environment, students require new technology, new methods, new instruments and even new learning techniques. Google apps is the choice for the students. This research aims to assess the Google apps that were selected from the previous works. The researchers applied Google apps in software engineering class of College of Art Media and Technology (CAMT).

Introduction

Cloud computing is one of the applicable technologies which can improve end-user productivity and reduce driven-overhead by offering services. Cloud computing technology has been widely defined in many different ways. Cloud computing not provides only services through internet but also the provision of several commodities "as-a-service" [START_REF] Accorsi | Business Process as a Service: Chances for Remote Auditing[END_REF][START_REF] Van Den Heuvel | Leveraging Business Process as a Service with Blueprinting[END_REF]. On another word, cloud computing can provide the provisioning services such as Infrastructure as a Service (IaaS), Platform as a Service (PaaS), Software as a Service (SaaS), and Business Process as a Service (BPaaS).

Cloud computing is also information technology communication for education environment [START_REF] Xin | Cloud Computing for Education: A Case Study of Using Google Docs in MBA Group Project[END_REF][START_REF] Pocatilu | Cloud Computing Benefits for E-learning Solutions[END_REF]. The students can get the benefits from this technology, especially, Software as a Service (SaaS) which can improve flexibility and accessibility [START_REF] Xin | Cloud Computing for Education: A Case Study of Using Google Docs in MBA Group Project[END_REF][START_REF] Godse | An Approach for Selecting Software-as-a-Service (SaaS) Product[END_REF]. Cost cutting is a key benefit because the users don't have to install the software on their own computer. SaaS change the way how to deliver the software to customers because the customers can access the services through internet [START_REF] Nitu | Configurability in SaaS (software as a service) applications[END_REF]. The organizations which are interested in the SaaS model would like to subscribe SaaS more than develop a new one [START_REF] Hirzalla | Realizing Business Agility Requirements through SOA and Cloud Computing[END_REF]. If the organizations would like to apply SaaS, they need to select appropriate services. Different SaaS providers propose different services to their educational solution. They offer not only different services; but also, different solutions. On the other hand, these solutions are not completely different, they still have common services. These reasons lead competition among SaaS providers occur. SaaS providers are growing up rapidly due to continuous growing of cloud computing technology. Consequence, SaaS vendors can offer several alternative solutions that able to meet SaaS user needs. Normally, SaaS providers provide standardized services to subscribers but the subscribers always have unique requirements for their own businesses. If the standardized services could not fulfill their requirements, customization and configuration are applicable [START_REF] Nitu | Configurability in SaaS (software as a service) applications[END_REF][START_REF] Wei | Software as a Service: Configuration and Customization Perspectives[END_REF]. Selection issue comes up with decision making. The decision making applies the prioritization to clearly define the selection [START_REF] Godse | An Approach for Selecting Software-as-a-Service (SaaS) Product[END_REF][START_REF] Aasem | Analysis and optimization of software requirements prioritization techniques[END_REF].

Higher education and cloud computing bring cost effectiveness to the educational institution [START_REF] Vitkar | Cloud Based Model for E-learning in Higher Education[END_REF]. For the higher education, the cost of being in this level is increasing compare to the lower level but the university's budget is decreasing [START_REF] Mircea | SOA, BPM and Cloud Computing: Connected for Innovation in Higher Education[END_REF][START_REF] Guin | A smart Architectural Concept for making of University Educatin System using Cloud Computing Paradigm[END_REF]. The students in the higher education also ask for improving the teaching methods, new technologies, new instrument and new learning technique [START_REF] Mircea | SOA, BPM and Cloud Computing: Connected for Innovation in Higher Education[END_REF][START_REF] Vitkar | Cloud Based Model for E-learning in Higher Education[END_REF]. Some education institutions do not have resources or infrastructure to support their students [START_REF] Pocatilu | Cloud Computing Benefits for E-learning Solutions[END_REF][START_REF] Peng | Strategies for e-education[END_REF]. In order to the needs of student and the problems of education institutions, adapting software services can help with this issue. Reusing software service leads us more efficient and more adaptable to the needs [START_REF] Mircea | SOA, BPM and Cloud Computing: Connected for Innovation in Higher Education[END_REF][START_REF] Guin | A smart Architectural Concept for making of University Educatin System using Cloud Computing Paradigm[END_REF]. Moreover, this research also adopt learning cloud services framework to build up the service checklists of students because this framework is under the perspective of learners [START_REF] Cheng-Sian | The Implications of Learning Cloud for Education: From the Perspectives of Learners[END_REF].

This research aims to assess the Google apps that were selected from the previous work. There are two previous works that have been related to this research paper: provider selection and service prioritization. Recommendation matrix for selecting education cloud service providers: configuration and customization perspectives, is applied to select the cloud provider [START_REF] Apitchaka | Prioritization Technique for Learning Cloud Services[END_REF]. The imperativeness of services is prioritized by using combined techniques between 100 dollars test and grouping [16].

In the following, we begin by showing the literature review of software selection process that we apply for the research methodology (Section 2). Next, we explain how this research were conducted (Section 3), this is followed by a results and discussions (Section 4) and conclusion (Section 5).

Software Selection Process

In 2009 Syed Ahsan Fahmi and Choi summarized the existing methods which are OTSO, CSSP, CRE and CBR [START_REF] Fahmi | A study on software component selection methods[END_REF]. The first one, Off-The Shelf-Option (OTSO) is the general method to choose the component. The aim of this method is evaluate the existing component-of-the-shelf (COTs) based on the requirement of the organization. This method consists of 6 tasks which are searching, screening, evaluation, analysis, deployment and assessment [START_REF] Constantine | SysML Modeling of Off-the-Shelf-Option Acquisition for Risk Mitigation in MilitaryPrograms[END_REF]. The second one, COTS Software Selection Process (CSSP), the aim of this method is to choose the right COTs for the organization. The important thing is this method has defined the criteria of input and output of each task clearly. The method is composed of 6 tasks which are evaluation team, goal identify, plan, filtering from vendor information, filtering from vendor demonstration and documentation [START_REF] Lin | COTS Software Selection Process[END_REF]. The third one, Cots-Based Requirement Engineering (CRE) is the systematic method and reusable follow the requirements. The highlight of this method is supporting non-functional requirement to choose software. The main point of this method is the requirement which meets the customer will be kept and others requirements are eliminated [START_REF] Alves | CRE: A Systematic Method for COTS Components Selection[END_REF]. The last one, Cots-Based Requirement Engineering (CBR), contains 4 steps which are retrieve, reuse, revise and retain. This method suitable with the organization where have experience about choosing the software before, therefore; the organization can select the software base on the experience [START_REF] Aamodt | Case-Based Reasoning: Foundational Issues, Methodological Variations, and System Approaches[END_REF]. Each software selection process has own appropriated context. OTSO is proper with the organization where decision making depends on authorized people. CSSP I suitable for the organization where explicit results of each selection phase is required. CRE is for the organization where they concern about users; satisfaction. CBR fits with the organization where they have experiences about choosing software. On the other hand, these software selection processes have the phases in common. As a result researcher classifies the common task of each method into 3 tasks which are selection, analysis and assesment as shown in Fig. 1. The Eliminate, Combine Rearrange, and Simplify (ECRS) method is adapted to conclude the research methodology in this research [START_REF] Haitao | Submerged Arc Welding Procedure Improvement based on Human-Machine Operation Analysis[END_REF]. Selection is about searching or identifying what kind of software that we need. Analysis is the designing process that we could filter or manage the software from the previous process. Assessment is the process that validate the selected software.

Research Methodology

The research methodology consists of three main processes: provider selection process, service analysis process and software as a services assessment process. In the process of selection, the principles of the learning cloud service framework are utilized. In the service analysis process, we use two techniques of prioritization: the 100 dollars test and grouping. This process classifies the services into groups by using K-means algorithm. In the course of software-as-s-service assessment, a requirement document is created to build up the prototype of the learning supported system for students. In this research we focus on the third process (SaaS Assessment) as shown in Fig. 2., there are 3 phases which are matching, designing and applying. Firstly, matching phase, the prioritized services and cloud service provider are matched. Secondly, Designing phase is the phase that rank the selected services into the web portal. Lastly, these selected services should be applied in the organization.

Fig. 2. The Overview of Proposed Framework

Fig. 2. Describes the proposed framework. The first step is about selection provider through recommendation matrix. The result is an appropriate provider for the organization. The second step is analyzing the services of the provider by adopting the learning cloud service requirement to classify these services into groups. After that let the student prioritize the services follow 100 dollars test technique. The K-means algorithm applied to identify the proper number of group of the services. The prioritized services are formed into the document and the SaaS is assessed by the students in the last step.

Matching Prioritized Services with the Selected Cloud Provider

The selected cloud provider in this research is Google. In addition, Google apps are the services that we will match with the prioritized services. The prioritized services are from learning cloud service checklist. The requirements of service checklist adopt the requirements of learning cloud service framework. 5 categories of services are personal learning, community learning, collaborative learning, communicative learning, and multi-rich material learning.

Designing Web Portal

The document is written to communicate with developer who can create the prototype and let the student try the collection of services that has been chosen. Regarding the matching prioritized services with the selected cloud provider, we know which Google apps will be used. Moreover, each service and each category has own priority.

With the priority, we can design the web portal.

Applying Learning Cloud Services in the Class

First of all, train the students about how to use the Google apps that has been chosen. And then, let the students try these Google apps for one semester .After that, the students could return their feedback. There are 5 scales that the student can assign to each services about their usability. The score is 5 if that services is strongly useful for their learning. The score is 4 if that services is useful for their learning. The score is 3 if that services is not affected for their learning. The score is 2 if that services is useless for their learning. The score is 1 if that services is strongly useless for their learning.

Result and Discussion

The researchers applied this approach in College of Art Media and Technology (CAMT) with 42 students. For the provider selection process, Google were chosen as a provider because Google got the highest recommendation score from the recommendation matrix for selecting education cloud service provider [16]. This matrix is built up from the configuration and customization and take the provider who provide the services for educational institution to be candidate. Moreover, the criteria of this matrix are from learning cloud service framework. Next the service analysis process apply two prioritization techniques. Numerical assignment technique allows students to answer 'need' and 'no need' and 'no opinion'. cumulative voting technique or 100 dollars technique is adapt in 'need' answer by allowing the student spend money on service that they need and they can spend money 100 dollars in total. Then the prioritized services are from the service analysis process. Services are classified into 3 groups as follows: critical, standard and optional. The optional services are not applied in this research.

4.1

Matching Prioritized Services with the Selected Cloud Provider Fig. 3. The Diagram for the Learning Support Services Fig. 3. is showing the relation between the services and software as a service of the provider. There are 14 learning cloud services as follows: translation, search, synchronization, file storage and sharing, wiki, discussion board, social network, blogging, internet telephony, e-mail, sms, chat, video on demand, and code. However, only 11 Google apps that matched with the selected cloud learning services which are Google translate, Google search, Google Syn, Google drive, Google site, Google group, Google+, Blogger, Gmail, Youtube for education, Google code. Google+ could match both social network and internet telephony. Gmail could provide services for e-mail, sms and chat.

Designing Web Portal

The result prioritization can design the webpage. The services in critical and standard group are applied in the learning support web portal as shown in Fig. 4. Five categories of services are also rank from the most important category to the least important category, personal learning, collaborative learning, social community learning, communicative learning, and multimedia-rich material respectively.

Fig. 4. The Overview of Learning Support Web Portal

The services in each category is classify into two groups which are the critical and standard. Firstly, personal learning has synchronization in standard group and other services in critical group. Secondly, collaborative learning has file storage &sharing in critical group and other services in standard group. Thirdly, social community learning has social network service in critical group and blogging service in standard group. Fourthly, communicative learning has internet telephony in critical group and other services in standard group. Lastly, multi-media rich material has only services in the standard group. 5. illustrates the example of the personalized learning page which arranges the services in the critical group to the left hand-side and arrange the services in the standard group in the right-hand. Moreover, the priority of each page are also prioritized from the amount of money that the students spent from the service analysis phase with the 100-dollar technique.

Applying Learning Cloud Services in the Class

For the feedback, 41 students from 42 students respond about the software as a service that the researchers select for them in Table1. The researchers provide an online questionnaire to let them evaluate each selected service about the supportedusability for their learning. Fig. 6. The Feedback from the students Fig. 6. demonstrate the feedback of students about the each selected service's supported-usability for their learning. The horizontal axis shown the priority of each services, while the vertical axis represent the feedback score of each service. The priority ranked from the heist priority to the lowest priority which are number 1 and number 14 respectively. The average score of these services are 3.81.

The trend of the students' feedback is the same direction like the prioritization services. Most of SaaS from the critical group got the score over average, except Google+. In this case, the integration approach is recommended; the organization can apply another SaaS to the social network approach because there are many social networking. Whereas Most of SaaS from the standard group got the score under average, only Youtube for education that slightly greater that the average. In this case, we found learning by video on demand is suitable for the software engineering class of CAMT. .

Conclusion

Applying SaaS in the classroom could enhances the learning method of the students. This paper tries to prove that the selected learning cloud services from the previous work could make the student satisfying. The process starts from matching the prioritized services with the selected cloud provider, and then designing the web portal and letting the students try before give the feedback. The limitation of this research is applying the services from only one cloud provider. For the future work, the researchers can use the integration approach with various sources of services.

Fig. 1 .

 1 Fig. 1. The Common Software Selection Process

Fig. 5 .

 5 Fig.5. The Personalized Learning Page Fig.5. illustrates the example of the personalized learning page which arranges the services in the critical group to the left hand-side and arrange the services in the standard group in the right-hand. Moreover, the priority of each page are also prioritized from the amount of money that the students spent from the service analysis phase with the 100-dollar technique.