In this section the finite element solution of the example problem is followed step-by-step to-gether with the correspondingMatlabcodes.
D.2.1 Parametrization
First, the global parameters are set, i.e. the properties of the material and the source frequency.
1% Material and excitation properties
2rho = 1.2; % Average density of air [kg/m^3]
3c = 343; % Speed of sound [m/s]
4f = 400; % Frequency of excitation [Hz]
5om = 2*pi*f; % Angular frequency [rad/s]
Next, the parameters of the geometry and the discretized model are configured. The element size is chosen asLe= 0.05 m, which roughly corresponds to16elements per wavelength. Using this configuration we haveNPML= 8elements along the thickness of the PML.
6% Geometry
7R_cyl = 0.5; % Radius of the cylinder [m]
8N_cyl = 60; % Number of elements on cylinder
9L_slab = 3.5; % Side length of computational domain [m]
10 Le = 0.05; % Length of elements [m]
11 N_slab = ceil(L_slab / Le); % Number of elements along the edge
D.2. FINITE ELEMENT SOLUTION CODE 143
12 tol = 1e-3; % Relative tolerance
13 % PML parameters
14 N_PML = 8; % Number of PML elements in the layer
15 L_PML = 0.4; % Thickness of PML layer [m]
16 L_sum = L_slab + 2*L_PML; % Total side length [m]
D.2.2 Mesh generation
Our mesh consists of two parts, i.e. the pure FE mesh, and the attached PML. The inner mesh is generated first, by creating its contour and using automatic triangular mesh generation. The latter is provided by themesh2dfunction of theMESH2D2toolbox.
17 % Create contour
18 cylinder = create_circle_boundary(R_cyl, N_cyl/2);
19 slab_in = translate_mesh( ...
20 create_slab_boundary(L_slab, N_slab/2), -1/2*L_slab*[1, 1, 0]);
21 contour = join_meshes(cylinder, slab_in);
22 % Create triangle mesh
23 [p,t] = mesh2d(contour.Nodes(:,2:3), contour.Elements(:,5:6), ...
24 struct(’hmax’, Le), struct(’output’, false));
25 mesh = pt2nihumesh(p, t); % NiHu mesh from points and triangles
The functionpt2nihumeshconverts the generated nodes and triangles into the mesh format ofNiHu. Note: Themesh2dfunction halves all the edges on the boundary, thus, in lines 18 and 20 the parametersNcylandNslabare halved.
The PML mesh is generated by creating a square with a side length ofLslab+ 2LPML, and cutting out its middle part. The latter is performed in lines 29–30 using the mesh_section function. After the middle elements are cut out, the corresponding nodes are eliminated by calling the functiondrop_unused_nodes.
26 % Create the PML mesh
27 pml_mesh = translate_mesh( ...
28 create_slab(L_sum, N_slab+2*N_PML), -1/2*(L_sum)*[1, 1, 0]);
29 pml_mesh = mesh_section(pml_mesh, ...
30 1/2*L_slab*[-1+tol, -1+tol, -inf; 1-tol, 1-tol, inf], ’none’);
31 pml_mesh = drop_unused_nodes(pml_mesh);
32 % Assign PML properties
33 pml_mesh.Elements(:,2) = 224; % PML elem identifier
34 pml_mesh.Properties = [1, 3, 1, L_PML, 0, 0];
The PML properties are assigned to the mesh in lines 33–34.
When both parts of the mesh are constructed, they are joined together and their coincident nodes are merged. The material properties are configured in line 38. Finally, the absorbing func-tion of the PML is configured in lines 39–40. The funcfunc-tionpml_sigma_globimplements the unbounded absorbing function as defined by equation (4.49).
35 % Join the meshes
36 model = merge_coincident_nodes(join_meshes(mesh, pml_mesh), 1e-3);
37 model = drop_mesh_IDs(model);
38 model.Materials(1, 3:4) = [rho, c];
39 model.PMLData = {1, @pml_sigma_glob, @pml_sigma_glob, [], ...
40 [-L_slab/2 L_slab/2], [-L_slab/2 L_slab/2], []};
2c Darren Engwirda, 2009. See: http://www.mathworks.com/matlabcentral/fileexchange/25555-mesh2d
The generated mesh is displayed in Figure D.1(b), with the PML part shown in transparent green for clarity. The numerical model has 6 541nodes and 10 182elements comprising 7 686 triangles in the inner part and2 496quadrangles in the absorbing layer.
D.2.3 System matrix assembly
In the next step, the system matrices of the FE–PML model are assembled. The matricesKand Mfrom equation (4.13) are constructed by calling themodel_mkfunction. The matricesK(PML) andM(PML)are assembled separately by callingpml_mk. Then, the two subsystems are joined by the functionmodel_system.
41 % System matrix assembly
42 [I, J, M, K, DOF] = model_mk(model, ’ind’); % Pure FEM
43 [Ip, Jp, Mp, Kp, DOFp] = pml_mk(model, om, ’ind’); % Pure PML
44 % Join subsystems
45 [M, K, DOF] = model_system({I J M K DOF}, {Ip Jp Mp Kp DOFp});
46 A = model_a(get_boundary(model)); % Excitation matrix
The reason for assembling the pure FE and the PML system matrices in two subsequent steps is that the former are independent of the frequency, while the latter are not. If we would solve the same problem with a number of different testing frequencies, the pure FE matrices would need to be assembled only once. Finally, the matrixAfrom equation (4.13) is assembled. The structure of the sparse system matricesK andMis displayed in Figure D.2(a). The size of the matrices and hence the number of DOFs of the system correspond the number of nodes in the mesh.
D.2.4 Excitation and boundary conditions
Before the definition of the excitatation the nodes on the boundary are selected in line 48 by the mesh_selectfunction.
47 % Find excitation indices
48 i_exc = mesh_select(model, sprintf(’r < %f’, R_cyl*(1+tol)), ’ind’);
The incident field generated by a plane wave is obtained by calling the functionincident with the parameter ’plane’. The function calculates the incident pressure and its normal derivative. For the latter, we make use of the fact that the normal vector onx ∈ Γis given as n(x) =−x/Rcyl. The normal particle velocity is attained using relation (4.3).
49 % Calculate incident field
50 [p_inc, q_inc] = incident(’plane’, [1 0 0], model.Nodes(: ,2:4), ...
51 -model.Nodes(:, 2:4) / R_cyl, om/c);
52 % Calculate the reflected normal velocity
53 v_refl = zeros(size(model.Nodes, 1), 1);
54 v_refl(i_exc) = -1./(1i*om*rho)*q_inc(i_exc);
D.2.5 Solution
The numerical solution is performed usingMatlab’s direct solver. In line 56 equation (4.12) is solved. After executing lines 56–57 the vectorpsolcontains the complex sound pressure values in the nodal coordinates of the mesh.
55 % Numerical solution
56 p_refl = (K - om^2*M)\(1i*om*A*v_refl);
57 p_sol = p_inc + p_refl;
D.2. FINITE ELEMENT SOLUTION CODE 145
(a) Structure of the system matrices (b) The numerical solution
Figure D.2.System matrix structure and the resulting acoustic field
The analytical solution is attained by evaluating equations (D.1)–(D.3). The infinite sum in (D.1) is truncated tommax= 100and evaluated by the functionplanewave_cyl2din line 59.
58 % Analytical solution
59 p_ana_refl = planewave_cyl2d(model.Nodes(:,2:3), R_cyl, om/c, 100);
60 p_ana = p_inc + p_ana_refl;
The numerical result is depicted in Figure D.2(b). As it is seen, the reflection from the cylin-drical obstacle perturbs the incident planar wavefield to a great extent. The thin white rectangle indicates the transition interface of standard finite elements and the perfectly matched layer. Go-ing outwards in the absorbGo-ing layer, the amplitude converges to the amplitude of the incident wave, as expected.
D.2.6 Validation
Finally, the error of the solution is computed by comparing the numerical and analytical results.
This is done by evaluating
ε= kpsol−panak
kpanak , (D.4)
with k·kdenoting the 2-norm. Since the amplitude of the reflected field decays exponentially inside the PML, the DOFs corresponding to the nodes inside the PML are excluded from the evaluation of (D.4). The nodes of the pure FE region are selected in lines 62–63.
61 % Find fem indices
62 i_fem = mesh_select(model, ...
63 sprintf(’abs(x) < %f & abs(y) < %f’, L_slab/2,L_slab/2), ’ind’);
64 % Calculate error
65 err = norm(p_ana(i_fem) - p_sol(i_fem))/norm(p_ana(i_fem));
66 fprintf(’Log10 L2 error of solution: %f\n’, log10(err));
The output reads as:
Log10 L2 error of solution: -1.704463
This corresponds to a relative error ofε ≈ 1.97%, which proves that our numerical solution is correct. The execution time of the whole example code is less than3 son an average desktop computer.
The above example demonstrated that a relatively complex acoustical problem involving fi-nite elements and perfectly matched layers can be modeled and solved from scratch in about 60 lines of code using theMatlabroutines of theNiHutoolbox. Besides the simplicity of the code, it was also seen that the toolbox also provides a considerable amount of insight into the internals of the finite element method, that are hidden by commercial FE software. The latter feature is espe-cially useful in education, where understanding how the solution process works is as important as the actual solution of the specific problem.