Tutorial 06 Second Order Elements
So far we have deployed linear elements for the FEM discretization of the Helmholtz equation. But WavesAndEigenvalues.Helmholtz can do more. Second order elements are also available. All you have to do is to set the order keyword in the call to order=:quad.
Model set-up
We demonstrate things with the usual Rijke tube set-up from Tutorial 01.
using WavesAndEigenvalues.Helmholtz
mesh=Mesh("Rijke_mm.msh",scale=0.001) #load mesh
dscrp=Dict() #initialize model descriptor
dscrp["Interior"]=(:interior, ()) #define resonant cavity
dscrp["Outlet"]=(:admittance, (:Y,1E15)) #specify outlet BC
γ=1.4 #ratio of specific heats
ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
Tu=300.0 #K unburnt gas temperature
Tb=1200.0 #K burnt gas temperature
P0=101325.0 # ambient pressure in Pa
A=pi*0.025^2 # cross sectional area of the tube
Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0
x_ref=[0.0; 0.0; -0.00101] #reference point
n_ref=[0.0; 0.0; 1.00] #directional unit vector of reference velocity
n=0.01 #interaction index
τ=0.001 #time delay
dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,:n,:τ,n,τ)) #flame dynamics
R=287.05 # J/(kg*K) specific gas constant (air)
speedofsound(x,y,z) = z<0. ? sqrt(γ*R*Tu) : sqrt(γ*R*Tb)
c=generate_field(mesh,speedofsound)3380-element Array{Float64,1}:
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
694.4379021914054
⋮
347.2189510957027
347.2189510957027
347.2189510957027
347.2189510957027
347.2189510957027
347.2189510957027
347.2189510957027
347.2189510957027
347.2189510957027The key difference is that we opt for second order elements during discretization
L=discretize(mesh,dscrp,c, order=:quad)6172×6172-dimensional operator family: ω^2*M+K+n*exp(-iωτ)*Q+ω*Y*C Parameters ---------- n 0.01 + 0.0im λ Inf + 0.0im ω 0.0 + 0.0im τ 0.001 + 0.0im Y 1.0e15 + 0.0im
All subsequent steps are the same. For instance, solving for an eigenvalue is done by
sol,nn,flag=mslp(L,340*2*pi,maxiter=20,tol=1E-11,output=true)
#Launching MSLP solver... Iter dz: z: ---------------------------------- 0 Inf 2136.2830044410593 1 382.9705879256757 1753.4014913373253 + 8.259427312643453im 2 42.502458544562025 1710.9235160596425 + 9.70185534924655im 3 0.5342818649578048 1710.3901346228909 + 9.732861060842836im 4 8.4413480680881e-5 1710.3900506855896 + 9.732870014343618im 5 1.5462854887275517e-11 1710.390050685605 + 9.73287001434383im 6 9.322369603784163e-12 1710.3900506856144 + 9.73287001434386im Solution has converged! ...finished MSLP! ##################### Results ##################### Number of steps: 6 Last step parameter variation:9.322369603784163e-12 Auxiliary eigenvalue λ residual (rhs):3.155516543983256e-8 Eigenvalue:1710.3900506856144 + 9.73287001434386im
And writing output to paraview by
data=Dict()
data["abs mode"]=abs.(sol.v)
vtk_write("tutorial06",mesh, data)Note However, that the vtk-file name will be tutorial06_quad.vtu.
Summary
Quadratic elements are an easy matter. Just set order=:quad when calling discretize.
This page was generated using Literate.jl.