Modeling Acoustics with Differential Equations

Yi-Lin Chiu,
Software Developer, WRI.

Join the Conversation #WolframTechConf

Partial differential equations (PDEs) Model

Acoustics Modeling with Wolfram Language

Frequency Domain Analysis

Define the domain

Ω=RegionUnion[chamber,pipeIn,pipeOut];

Specify boundary conditions.

(Subscript[Γ, RBC])=NeumannValue[(-I ω p[x,y,z]+2I ω pIn)/(ρ c),z==inlet];                   (*Radiation BC at inlet*)
(Subscript[Γ, ABC])=NeumannValue[-I ω p[x,y,z]/(ρ c),z==outlet]                                             (*Absorbing BC at outlet*)

Set the PDE model

FrequencyAcousticsModel3D=AcousticModel[p[x,y,z],{x,y,z},ρ,{0,{0,0,0}}]-(ω^2) p[x,y,z]/(ρ  (c^2));
pde ={FrequencyAcousticsModel3D==(Subscript[Γ, ABC])+(Subscript[Γ, RBC])}/.{c->(Subscript[c, air]),ρ->(Subscript[ρ, air])};

Solve

pfun=ParametricNDSolveValue[pde,p,{x,y,z}∈mesh,{ω},Method->{"PDEDiscretization"->{"FiniteElement"}}];
pfunTable=Table[pfun[ω]/.{ω-> 2π f},{f,50,1500,25}];
{vals,modes}=NDEigensystem[{AcousticModel[p[x,y],{x,y},(Subscript[ρ, air]),{0,{0,0}}]},p,{x,y}∈chamberRegion,5];

Post-processing and Visualization

(Legended[, (BarLegend[({{"sound amplitude (Pa)"}, {}})])])

Time Domain Analysis

Define the domain

Ω=RegionUnion[Disk[{0,cl},radius,{0,π}],Rectangle[{-d/2-cw/2,bottom},{-d/2+cw/2,cl}],Rectangle[{d/2-cw/2,bottom},{d/2+cw/2,cl}]];

Specify initial and boundary conditions.

ics={p[0,x,y]==0,Derivative[1,0,0][p][0,x,y]==0};
(Subscript[Γ, RBC])=NeumannValue[1/(ρ c)*(2D[pIn,t]-Derivative[1,0,0][p][t,x,y]),y==bottom];
(Subscript[Γ, ABC])=NeumannValue[-1/(ρ c)*Derivative[1,0,0][p][t,x,y],y⩾cl];

Set the PDE model

TimeAcousticsModel2D=(1/(ρ (c^2)))*D[p[t,x,y],{t,2}]+AcousticModel[p[t,x,y],{x,y},ρ,{0,{0,0}}];
pde ={TimeAcousticsModel2D==(Subscript[Γ, ABC])+(Subscript[Γ, RBC]),ics}/.{c->(Subscript[c, air]),ρ->(Subscript[ρ, air])};

Solve

pfun=NDSolveValue[pde,p,{t,0,tend},{x,y}∈Ω,AcousticsOptions[(Subscript[c, air]),ω,"Dimension"->2]]

Post-processing and Visualization

nframes=50;
frames=Table[ContourPlot[pfun[t,x,y],{x,y}∈Ω,AspectRatio->Automatic,Evaluate[opts],FrameLabel->{"x","y"},PlotLabel->Style["Diffraction of Sound",Bold]],{t,0,tend,tend/nframes}];
rasterizedFrames=Map[Rasterize[#,"Image",ImageResolution->150]&,frames];
ListAnimate[rasterizedFrames,SaveDefinitions->True]

Acoustics PDE setup

Acoustics PDE setup

Time domain Analysis --- Wave equation

(1/(ρ (c^2)))((∂^2)/(∂(t^2)))p(t,X)+∇·{-(1/ρ)[∇p(t,X)+F]}=Q

TimeAcousticsModel=(1/(ρ (c^2)))*D[p[t,x],{t,2}]+AcousticModel[p[t,x],{x},ρ,{0,{0}}];

Frequency domain analysis --- Helmholtz equation

-((ω^2)/(ρ (c^2)))p(X)+∇·{-(1/ρ)[∇p(X)+F]}=Q

FrequencyAcousticsModel=(-(ω^2)/(ρ (c^2)))*p[x]+AcousticModel[p[x],{x},ρ,{0,{0}}];

p: Sound pressure
ρ: Density
c: Speed of sound
t:  Time
X: Position vector
ω: Angular frequency (=2π f)
Q,F: Monopole and Dipole sound source

Sound sources

Example: Monopole source

({{\[Piecewise], ({{Time Domain PDE:, (1/(ρ (c^2)))((∂^2)/(∂(t^2)))p(t,X)+∇·{-(1/ρ)[∇p(t,X)+F]}=Q}, {Frequency Domain PDE :     , -((ω^2)/(ρ (c^2)))p(X)+∇·{-(1/ρ)[∇p(X)+F]}=Q}})}})

Monopole source:Q=(Subscript[\[LeftBracketingBar]Q\[RightBracketingBar], amp])·δ(X-(Subscript[X, s]))

Monopole Source Setup

(Subscript[Q, amp])=1;xs=2.5;ys=2.5;f=200; 
Q2D[t_,x_,y_,xs_,ys_]:=RegularizedDelta2D[h,x,y,xs,ys]*(Subscript[Q, amp])*Cos[ω t]/.{ω-> 2π f};
TimeAcousticsModel2D=(1/(ρ (c^2)))*D[p[t,x,y],{t,2}]+AcousticModel[p[t,x,y],{x,y},ρ,{Q2D[t,x,y,xs,ys],{0}}];

Example: Dipole source

({{\[Piecewise], ({{Time Domain PDE:, (1/(ρ (c^2)))((∂^2)/(∂(t^2)))p(t,X)+∇·{-(1/ρ)[∇p(t,X)+F]}=Q}, {Frequency Domain PDE:     , -((ω^2)/(ρ (c^2)))p(X)+∇·{-(1/ρ)[∇p(X)+F]}=Q}})}})

Dipole source:(Overscript[F, →])=(Subscript[\[LeftBracketingBar]F\[RightBracketingBar], amp])·δ(X-(Subscript[X, s]))·(cos θ,sin θ)

Dipole Source Setup

(Subscript[F, amp])=1;xs=2.5;ys=2.5;f=200;θ=π/4;
F2D[t_,x_,y_,xs_,ys_]:=RegularizedDelta2D[h,x,y,xs,ys]*(Subscript[F, amp])*{Cos[θ],Sin[θ]}*Cos[ω t]/.{ω-> 2π f};
TimeAcousticsModel2D=(1/(ρ (c^2)))*D[p[t,x,y],{t,2}]+AcousticModel[p[t,x,y],{x,y},ρ,{0,F2D[t,x,y,xs,ys]}];

Sound medium

Example: Lossy media

(Subscript[R, f]): flow resistivity, ϕ: porosity                        ({{\[Piecewise], ({{Time Domain PDE:, ((ϕ (Subscript[R, f]))/((ρ^2) (c^2)))(∂/(∂t))p(t,X)+(1/(ρ (c^2)))((∂^2)/(∂(t^2)))p(t,X)+∇·{-(1/ρ)[∇p(t,X)+F]}=Q}, {Frequency Domain PDE:     , I((ϕ (Subscript[R, f])ω)/((ρ^2) (c^2)))p(X)-((ω^2)/(ρ (c^2)))p(X)+∇·{-(1/ρ)[∇p(X)+F]}=Q}})}})

Setup

Ω=Line[{{left},{right}}];
TimeAcousticsModelLossy=D[p[t,x],{t,2}]/(ρ (c^2))+AcousticModel[p[t,x],{x},ρ,{0,{0,0}}]+ϕ Rf*D[p[t,x],{t,1}]/((ρ^2)(c^2));
(Subscript[Γ, D])=DirichletCondition[p[t,x]==Sin[ω t],x==left];

Solve & Visualize

pde ={TimeAcousticsModelLossy==0,silentIC,(Subscript[Γ, D])}/.{c->(Subscript[c, air]),ρ->(Subscript[ρ, air]),ϕ-> 0.617,Rf-> 5000};
pfun=NDSolveValue[pde,p,{t,0,tend},x∈Ω,AcousticsOptions[(Subscript[c, air]),ω]];

nframes=50;
frames=Table[Plot[{pfun[t,x],pLosslessAna[t,x]},{x}∈Ω,PlotLegends->{Style["Lossy media",18],Style["Lossless media",18]},PlotRange->{-1.2,1.2},AxesLabel->{Style["x",18],Style["sound pressure",18]},PlotStyle->{{},LightGray},ImageSize->Large],{t,0,tend,tend/nframes}];
ListAnimate[frames,SaveDefinitions->True,ImageSize->Automatic]

Boundary Conditions

Boundary Conditions in Acoustics

Example: Impedance boundary conditions (IBC)

(Subscript[Z, b]):Impedance at the boundary                                                         ({{\[Piecewise], ({{Time Domain IBC:, (OverVector[n])·( (1/ρ)∇p(t,X))=-(1/(Subscript[Z, b]))(( ∂p)/(∂t))}, {Frequency Domain IBC:     , (OverVector[n])·((1/ρ)∇p(X)) =-I ω ((p(X))/(Subscript[Z, b]))}})}})          on ∂Ω

Setup

Ω=Line[{{0},{0.5}}];
(Subscript[Z, b])=(Subscript[c, 2])*(Subscript[ρ, 2]);
(Subscript[Γ, IBC])=NeumannValue[-1/(Subscript[Z, b])*Derivative[1,0][p][t,x],x==0.5];
ics={p[0,x]==2Exp[-(x^2)/(0.05^2)],Derivative[1,0][p][0,x]==0};

Solve & Visualize

pde ={TimeAcousticsModel==(Subscript[Γ, IBC]),ics}/.{c->(Subscript[c, 1]),ρ->(Subscript[ρ, 1])};
pfun=NDSolveValue[pde,p,{t,0,tend},x∈Ω,AcousticsOptions[(Subscript[c, air]),0,"WaveType"-> "TimeInharmonic","StandardDeviation"-> β]];

Perfectly Matched Layer (PML)

Example: Perfectly Matched Layer (PML)

(∂/(∂x))(-(1/ρ)(((Subscript[s, y])(Subscript[s, z]))/(Subscript[s, x]))((∂p)/(∂x)))+(∂/(∂y))(((-1)/ρ)(((Subscript[s, x])(Subscript[s, z]))/(Subscript[s, y]))((∂p)/(∂y)))+(∂/(∂z))(((-1)/ρ)(((Subscript[s, x])(Subscript[s, y]))/(Subscript[s, z]))((∂p)/(∂y)))-(((Subscript[s, x])(Subscript[s, y])(Subscript[s, z])(ω^2))/(ρ (c^2)))p+∇·F=Q

(Subscript[σ, x, y,  z]):Absorbing coefficients of the PML                (Subscript[s, x, y,  z])=1-I ((c (Subscript[σ, x, y,  z]))/ω): Auxiliary variables

Setup

pmlWidth=0.5;σMax=35;
(Subscript[Ω, extended])= Line[{{left},{right+pmlWidth}}];
FrequencyAcousticsModelPML=AcousticModel[p[x],{x},ρ*sx,{0,{0}}]-sx (ω^2)p[x]/(ρ (c^2));
(Subscript[Γ, RBC])=NeumannValue[(-I ω p[x]+2I ω pIn)/(ρ c),x==left];

Solve & Visualize

pde ={FrequencyAcousticsModelPML==(Subscript[Γ, RBC])}//.{c-> (Subscript[c, air]),ρ->(Subscript[ρ, air]),sx->1-I c σx[σMax]/ω};
pfun=ParametricNDSolveValue[pde,p,{x}∈(Subscript[Ω, extended]),{ω},AcousticsOptions[(Subscript[c, air]),ω]];

f=1500;
pfunFrequencyDomain=pfun[ω]/.{ω->2π f};
pfunTimeDomain[t_,x_]:=Re[pfunFrequencyDomain[x]*Exp[I ω t]/.{ω->2π f}]

Application: Modeling an Automotive Muffler

({{\[Piecewise], ({{Radius of the chamber:, R=75mm}, {Radius of inlet/outlet pipes:     , r=40 mm                                           }, {Length of the chamber, L=600mm}, {Length of inlet/outlet pipes, l=40 mm}})}})

Define the domain

Ω=RegionUnion[chamber,pipeIn,pipeOut];

Specify boundary conditions.

pIn=1;
(Subscript[Γ, RBC])=NeumannValue[(-I ω p[x,y,z]+2I ω pIn)/(ρ c),z==inlet];                   (*Radiation BC at inlet*)
(Subscript[Γ, ABC])=NeumannValue[-I ω p[x,y,z]/(ρ c),z==outlet]                                             (*Absorbing BC at outlet*)

Set the PDE model

FrequencyAcousticsModel3D=AcousticModel[p[x,y,z],{x,y,z},ρ,{0,{0,0,0}}]-(ω^2) p[x,y,z]/(ρ  (c^2));
pde ={FrequencyAcousticsModel3D==(Subscript[Γ, ABC])+(Subscript[Γ, RBC])}/.{c->(Subscript[c, air]),ρ->(Subscript[ρ, air])};

Solve

pfun=ParametricNDSolveValue[pde,p,{x,y,z}∈Ω,{ω},Method->{"PDEDiscretization"->{"FiniteElement"}}];
pfunTable=Table[pfun[ω]/.{ω-> 2π f},{f,50,1500,25}];

Application: Modeling an Automotive Muffler

Post-processing and Visualization

Sound pressure field

Sound amplitude:         \[LeftBracketingBar]p(X)\[RightBracketingBar]=(Sqrt[p(X)(OverBar[p])(X)])

fRange=Table[f,{f,50,1500,25}];
pAmplitude[i_]:=Sqrt[pfunTable[[i]][x,y,z]*Conjugate[pfunTable[[i]][x,y,z]]];

(Legended[, (BarLegend[({{"sound pressure amplitude (Pa)"}, {}})])])

Application: Modeling an Automotive Muffler

Transmission loss of the muffler

({{\[Piecewise], ({{Time-averaged sound intensity:, I=(((\[LeftBracketingBar]p\[RightBracketingBar])^2)/(2ρ c))}, {Sound power:    , W=(∫_( ∂Ω)) I\[DifferentialD]A}, {Transmission-loss:, TL=10(Subscript[log, 10])(((Subscript[W, in])/(Subscript[W, out]))).}})}})

Iin=(pIn^2)/(2(Subscript[ρ, air]) (Subscript[c, air]));
Iout=Table[(((pAmplitude[i])^2)/(2(Subscript[ρ, air]) (Subscript[c, air])))/.{x-> 0,y-> 0,z-> outlet},{i,fRange//Length}];

Win=Iin*pipeArea;
Wout=Iout*pipeArea;

TL=10*Log[10,Win/Wout];

ListLinePlot[Transpose[{fRange,TL}],PlotLabel->Style["Transmission loss",Bold,24], AxesLabel->{Style["frequency",18],Style["db",18]},ImageSize-> Large]

Application: Modeling an Automotive Muffler

Eigenfrequency Analysis

Eigenfrequency:         (Subscript[f, i])=((c(Sqrt[ρ (Subscript[λ, i])]))/(2π ))

{vals,modes}=NDEigensystem[{AcousticModel[p[x,y],{x,y},(Subscript[ρ, air]),{0,{0,0}}]},p,{x,y}∈chamberRegion,5];
EigenfrequencyList=(c/(2π))*Sqrt[ρ*vals]/.{c-> (Subscript[c, air]),ρ-> (Subscript[ρ, air])}

modeAmplitude[i_]:=Sqrt[modes[[i]][x,y]*Conjugate[modes[[i]][x,y]]];

bl=BarLegend[{"TemperatureMap",{0,9}},50,LegendLabel->Style["sound amplitude (Pa)",18]];
opts={PlotRange->{0,9},ColorFunction->ColorData[{"TemperatureMap",{0,9}}],ContourStyle-> None,ColorFunctionScaling->False,AspectRatio->Automatic,Contours->50,PlotPoints->41,ImageSize-> Large,PlotLabel->Style["Eigenmode: f = 1241.9 Hz",Bold,24],FrameLabel->{Style["x",18],Style["y",18]}};
Legended[ContourPlot[modeAmplitude[4],{x,y}∈chamberRegion,Evaluate[opts]],bl]

Future goals

Thank you for your attention!

Which field should be the next to tackle?

Questions? Comments? Suggestions?

yilinc@wolfram.com

Init

Needs["NDSolve`FEM`"]

AcousticModel[p_,X:{x__},ρ_,{Q_,F_}]:=Module[{a,factor},
factor=-1/ρ*IdentityMatrix[Length[X]];
a=PiecewiseExpand[Piecewise[{{factor,True}}]];
Inactive[Div][a.Inactive[Grad][p,X],X]-Inactive[Div][F/ρ,X]-Q]

MaxGridSize[c_,ω_,Resolution_]:=
Module[{λ},
λ=2π c/ω;
λ/Resolution]

Options[AcousticsOptions]={"Resolution"-> 20,"Dimension"-> 1};
AcousticsOptions[c_,ω_,opts:OptionsPattern[]]:=Sequence[Method->{"PDEDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->MaxGridSize[c,ω,OptionValue["Resolution"]]^OptionValue["Dimension"]}}}]

(Subscript[c, air])=QuantityMagnitude[ThermodynamicData["Air","SoundSpeed"]];
(Subscript[ρ, air])=QuantityMagnitude[ThermodynamicData["Air","Density"]];
tend=0.01;

(*wave diffraction example*)
bottom=0;radius=0.6;cl=1/4;cw=3/100;d=1/5;

(*muffler example*)
ratio=1000;
R=75/ratio;
L=600/ratio;
r=40/ratio;
l=150/ratio;
inlet=-l; outlet=L+l;
pipeArea=π (r^2);
offset=10/ratio;
chamberRegion=BoundaryDiscretizeRegion[RegionUnion[Disk[{-R,0},R],Rectangle[{-R,-R},{R,R}],Disk[{R,0},R]]];
chamber=RegionBoundary[RegionProduct[chamberRegion,BoundaryDiscretizeRegion@Line[{{0},{L}}]]];
pipeIn=RegionBoundary[RegionProduct[BoundaryDiscretizeRegion[Disk[{0,0},r],MaxCellMeasure->Infinity],BoundaryDiscretizeRegion[Line[{{inlet},{offset}}]]]];
pipeOut=RegionBoundary[RegionProduct[BoundaryDiscretizeRegion[Disk[{0,0},r],MaxCellMeasure->Infinity],BoundaryDiscretizeRegion@Line[{{L-offset},{outlet}}]]];

bl=BarLegend[{"TemperatureMap",{0,1.65}},50,LegendLabel->"sound pressure amplitude (Pa)"];
opts={PlotRange->{0,1.65},ColorFunction->ColorData[{"TemperatureMap",{0,1.65}}],ImageSize-> Large,ContourStyle-> None,ColorFunctionScaling->False,Contours->30,PlotPoints->41,Boxed-> False,Axes->False,BoxRatios-> Automatic,PlotLabel->"f = 1500 Hz" };
Legended[SliceContourPlot3D[pAmplitude[-1],{x,y,z}∈mesh ,Evaluate[opts]],bl]

(*Time Domain Analysis*)

ClearAll[MaxGridSize]
MaxGridSize[c_,ω_,WaveType_,Resolution_,StandardDeviation_]:=
Module[{λ,β},
β=StandardDeviation;
λ=If[WaveType=="TimeHarmonic",2π c/ω,2β(Sqrt[2Log[10]])];
λ/Resolution]

ClearAll[AcousticsOptions]
Options[AcousticsOptions]={"WaveType"-> "TimeHarmonic","Resolution"-> 20,"StandardDeviation"-> 0.1,"Dimension"-> 1};
AcousticsOptions[c_,ω_,opts:OptionsPattern[]]:=Sequence[AccuracyGoal->5,PrecisionGoal->5,Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->MaxGridSize[c,ω,OptionValue["WaveType"],OptionValue["Resolution"],OptionValue["StandardDeviation"]]^OptionValue["Dimension"]}}}}]

(*Lossy media example*)
left=0;right=1;ω=5025;tend=0.003;
silentIC = {p[0,x]==0,Derivative[1,0][p][0,x]==0};
pLosslessAna[t_,x_]:=If[x⩽(Subscript[c, air])t,Re[Exp[I(-(ω/(Subscript[c, air]))x+ω t-π/2)]],0]

(*Impedance BC example*)
left=0;right=1;ω=5025;tend=0.0025;
β=0.05;
p0[x_]:=2*Exp[-(x^2)/(β^2)];
ics={p[0,x]==p0[x],Derivative[1,0][p][0,x]==D[p0[x],t]};
(Subscript[c, 1])=QuantityMagnitude[ThermodynamicData["Air","SoundSpeed",{"Pressure"->Quantity[101325,"Pascals"], "Temperature"->Quantity[250,"Kelvins"]}]];
(Subscript[ρ, 1])=QuantityMagnitude[ThermodynamicData["Air","Density",{"Pressure"->Quantity[101325,"Pascals"], "Temperature"->Quantity[250,"Kelvins"]}]];
(Subscript[c, 2])=QuantityMagnitude[ThermodynamicData["Air","SoundSpeed",{"Pressure"->Quantity[101325,"Pascals"], "Temperature"->Quantity[600,"Kelvins"]}]];
(Subscript[ρ, 2])=QuantityMagnitude[ThermodynamicData["Air","Density",{"Pressure"->Quantity[101325,"Pascals"], "Temperature"->Quantity[600,"Kelvins"]}]];

Ω= Line[{{left},{right}}];
TimeAcousticsModel=(1/(ρ (c^2)))*D[p[t,x],{t,2}]+AcousticModel[p[t,x],{x},ρ,{0,{0}}];
(Subscript[Γ, ABC])=NeumannValue[-1/(ρ c)*Derivative[1,0][p][t,x],x==right];
pde ={TimeAcousticsModel==(Subscript[Γ, ABC]),ics}/.{c->If[x≤1/2,(Subscript[c, 1]),(Subscript[c, 2])],ρ->If[x≤1/2,(Subscript[ρ, 1]),(Subscript[ρ, 2])]};
pfun=NDSolveValue[pde,p,{t,0,tend},x∈Ω,AcousticsOptions[(Subscript[c, air]),0,"WaveType"-> "TimeInharmonic","StandardDeviation"-> β]];

nframes=50;
frames=Table[Show[Plot[pfun[t,x],{x,0,1},ImageSize-> {500,250},PlotRange->{-0.5,1.5}],Graphics[{{Dashed,LightGray,Line[{{1/2,0},{1/2,3}}]},Text[Style["(Subscript[ρ, 1]), (Subscript[c, 1])",20],{0.25,1.2}],Text[Style["(Subscript[ρ, 2]), (Subscript[c, 2])",20],{0.75,1.2}],Text[Style["Medium 1",Bold,20],{0.25,1.4}],Text[Style["Medium 2",Bold,20],{0.75,1.4}]}]],{t,0,tend,tend/nframes}];
ListAnimate[frames,SaveDefinitions->True,ImageSize->Automatic]

nframes=50;
frames=Table[Show[Plot[pfun[t,x],{x}∈Ω,ImageSize-> {500,250},PlotRange->{-0.5,1.5}],Graphics[{{Dashed,LightGray,Line[{{1/2,0},{1/2,3}}]},Text[Style["(Subscript[ρ, 1]), (Subscript[c, 1])",20],{0.25,1.2}],Text[Style["Medium 1",Bold,20],{0.25,1.4}]}]],{t,0,tend,tend/nframes}];
ListAnimate[frames,SaveDefinitions->True,ImageSize->Automatic]

(*Impedance BC example*)
left=0;right=1;ω=2π*1500;tend=0.0025;
(OverHat[p])=Sin[ω t];
(Subscript[Γ, RBC])=NeumannValue[1/(ρ c)*(2D[(OverHat[p]),t]-Derivative[1,0][p][t,x]),x==left];
silentIC = {p[0,x]==0,Derivative[1,0][p][0,x]==0};
TimeAcousticsModel=(1/(ρ (c^2)))*D[p[t,x],{t,2}]+AcousticModel[p[t,x],{x},ρ,{0,{0}}];
pde ={TimeAcousticsModel==(Subscript[Γ, RBC]),silentIC}/.{c->(Subscript[c, air]),ρ->(Subscript[ρ, air])};
pfun=NDSolveValue[pde,p,{t,0,2tend},{x,left,right},AcousticsOptions[(Subscript[c, air]),ω,"Resolution"-> 50]];

nframes=50;
frames=Table[Show[Plot[pfun[t,x],{x,left,right},ImageSize-> {500,250},PlotRange->{-2.2,2.2}],Graphics[{{Dashed,Gray,Line[{{1,-2.2},{1,2.2}}]},Text[Style["Simulation Region",Bold,18],{0.5,1.8}]}]],{t,0.5tend,1.75tend,tend/nframes}];
ListAnimate[frames,SaveDefinitions->True,ImageSize->Automatic]

σx[σMax_]:=With[{right=right,pmlWidth=pmlWidth},If[x-right≥0,((x-right)/pmlWidth)^2*σMax,0]];

tend=0.0025

nframes=50;
frames=Table[Show[Plot[pfunTimeDomain[t,x],{x}∈(Subscript[Ω, extended]),ImageSize-> {500,250},PlotRange->{-1.5,1.5}],Graphics[{{Dashed,Gray,Line[{{1,-1.5},{1,1.5}}]},Text[Style["Simulation Region",Bold,18],{0.5,1.3}],Text[Style["PML Region",Bold,18],{1.3,1.3}]}]],{t,0,tend,tend/nframes}];
ListAnimate[frames,SaveDefinitions->True,ImageSize->Automatic]