Adjungierten-basierte Optimierung von akustischen Berandungen
im Zeitbereich mittels Volumenpenalisierung
Arne H¨
olter1, Emanuele Porcinai1, Mathias Lemke2, Stefan Weinzierl1
1TU Berlin, Fachgebiet Audiokommunikation, Email: hoelter@campus.tu-berlin.de
2TU Berlin, Fachgebiet Numerische Fluiddynamik
Einleitung
F¨
ur akustische Simulationen werden in der Raum-
akustik meist strahlen- oder frequenzbasierte Metho-
den eingesetzt. Diese ben¨
otigen vergleichsweise gerin-
ge Rechenleistungen, sind jedoch in ihrer Genauig-
keit begrenzt. Wellen- bzw. zeitbasierte Methoden wie
die finite Differenzenmethode im Zeitbereich (FDTD)
k¨
onnen eine deutlich h¨
ohere Genauigkeit erreichen,
allerdings stellt insbesondere die Beschreibung von
akustischen Berandungen eine Herausforderung dar.
Ein vielversprechender Ansatz hierf¨
ur ist die Verwen-
dung von Volumenpenalisierungsmethoden, wie sie in
der Str¨
omungsmechanik Anwendung finden [1]. Hierbei
wird der Str¨
omungswiderstand (Darcy Term) und die
Ver¨
anderung des effektiven Volumens in die beschreiben-
den Gleichungen eingebracht [2; 3; 4].
In diesem Beitrag wird die Optimierung des Darcy Terms
mittels adjungierten-basierter 2D-Simulationen an zwei
Beispielen untersucht. Im ersten Beispiel wird ein Ab-
sorptionselement auf Basis diskreter Mikrofonmessungen
identifiziert. Im zweiten Beispiel wird eine gegebene ener-
gy decay curve (EDC) in einem Rechteckraum durch ge-
eignete Berandung erzeugt.
Beschreibende Gleichungen
Als beschreibende Gleichungen f¨
ur die pr¨
asentierten Bei-
spiele dienen die nicht-linearen Euler-Gleichungen
∂t(ϱ) + ∂xj(φϱui)
φ= 0 (1)
∂t(ϱuj) + ∂xj(φϱuiuj)
φ+∂xj(pδij) = −χuj(2)
∂tp
γ−1+
∂xjφuipγ
γ−1−ui∂xj(φp)
φ= 0,(3)
mit der Dichte ϱ, der Geschwindigkeit ujin xj-
Richtung, dem Druck pund den Anpassungen durch
ein effektives Volumen φ=Vfluid/Vtotal, sowie dem
Str¨
omungswiderstand χ. Eine detaillierte Beschreibung
von φund χist in [3] gegeben. Es sei erw¨
ahnt, dass die
gezeigten Beispiele auch mit linearen akustischen Glei-
chungen erzeugt werden k¨
onnen
∂t(p′) + ϱ0c2
0∂xj(φu′
j)
φ= 0 (4)
∂t(u′
j) + ∂xj(p′)
ϱ0
=−χu′
j
ϱ0
.(5)
Adjungierten-basierter Ansatz
Die Optimierung der Volumenpenalisierung erfolgt mit-
hilfe eines adjungierten-basierten Verfahrens. Dabei wird
das Optimierungsproblem als Minimierungsaufgabe for-
muliert. Das Zielfunktional
J=gTq(6)
mit dem Zustandsvektor q= [ϱ, ui, p] und einem Gewich-
tungsfaktor gbeschreibt die N¨
ahe der L¨
osung zu einem
gew¨
unschten Zielwert. Der Faktor gwird durch die ad-
jungierte Gleichung definiert
ATq∗=g, (7)
wobei ATdie transponierte Systemmatrix der beschrei-
benden Gleichungen (1)-(3), vereinfacht als Aq =χ, und
q∗die adjungierte Zustandsvariable darstellt. Mittels der
Beziehung
J=gTq=ATq∗Tq=q∗TAq =q∗Tχ(8)
und L¨
osung der adjungierten Gleichung l¨
asst sich ein
einfacher Zusammenhang zwischen dem Zielfunktional
Jund dem Str¨
omungswiderstand χherstellen. Der ad-
jungierte Zustandsvektor q∗kann dementsprechend als
Gradient interpretiert werden. Der Str¨
omungswiderstand
wird damit iterativ im Sinne eines einfachen Gradienten-
verfahrens optimiert
χn+1 =χn+α∇χJ, (9)
wobei der Parameter αdie Optimierungsschrittweite dar-
stellt. F¨
ur Details zu den adjungierten Euler-Gleichungen
siehe [5; 6].
Simulationseinstellungen
Zur Simulation der Testf¨
alle wird als r¨
aumliche Ab-
leitung ein implizites Verfahren 4. Ordnung verwen-
det. Die zeitliche Integration wird mit einem Runge-
Kutta-Verfahren 4. Ordnung durchgef¨
uhrt. Zus¨
atzlich
wird nach jedem Zeitschritt ein impliziter, r¨
aumlicher
Filter 6. Ordnung zur Vermeidung numerischer Insta-
bilit¨
aten angewendet. Eine Analyse der Verfahren ist
in [7] zu finden. Als Randbedingungen werden charak-
teristische, nicht-reflektierende Formulierungen mit ei-
nem zus¨
atzlichen sponge-Bereich im Fall des einzelnen
Str¨
omungswiderstands verwendet. Zur Anregung des Re-
chengebiets dienen Gauss-Pulse als Anfangsbedingung.
Das Rechengebiet besteht aus gleichm¨
aßig verteilten Git-
terpunkten mit einem Abstand von ∆x= 0,01 m zwi-
schen den Gitterpunkten. Die Samplingfrequenz betr¨
agt
50 kHz, wodurch eine CFL-Zahl von ≈0,68 resultiert.
DAGA 2023 Hamburg
1070
0.0 0.5 1.0 1.5 2.0 2.5 3.0
x1in m
0.0
0.5
1.0
1.5
2.0
2.5
3.0
x2in m
mic 1
mic 2
source
sponge
absorber
Abb. 1: Schematischer Aufbau des Rechengebiets zur Opti-
mierung des Str¨
omungswiderstands. Im grauen Bereich wird
der sponge angewendet. Der gr¨
une Punkt stellt die Quellposi-
tion, die roten Punkte die Empf¨
angerpositionen und der blaue
Bereich den Str¨
omungswiderstand dar.
Identifikation eines einzelnen Str¨
omungs-
widerstands χ
Das Rechengebiet dieses Testfalls erstreckt sich ¨
uber
3 m ×3 m und es wird als Referenz ein Str¨
omungs-
widerstand bei x1= [1,85; 2,05] m und x2= [1,0; 2,0] m
mit einer Amplitude von χ= 5000 verwendet. In
Abb. 1 ist der Aufbau des Rechengebiets dargestellt.
In der Optimierung ist die Amplitude des Referenz-
str¨
omungswiderstandes unbekannt, lediglich die Druck-
signale an den Empf¨
angerpunkten sind gegeben. In
Abb. 4a ist der Referenzstr¨
omungswiderstand dargestellt
und in Abb. 4b das Ergebnis der Optimierung. Es ist
zu erkennen, dass die reflektierten Wellen durch einen
gleichm¨
aßigen Hauptstr¨
omungswiderstand reproduziert
werden. Der Wert der Referenz von χ= 5000 wird
nur in den vorderen 10 cm des optimierten Absor-
bers ben¨
otigt. Im hinteren Teil des Absorbers f¨
allt der
Str¨
omungswiderstand auf etwa χ= 4000 ab. Abb. 2 zeigt
eine sehr gute ¨
Ubereinstimmung der Rekonstruktion der
Referenzsignale. Sowohl die Reflexion (mic 1), als auch
die Transmission (mic 2) k¨
onnen exakt rekonstruiert wer-
den. Die G¨
ute der Rekonstruktion der Referenzsignale
bleibt von χ= 0 (keine Reflexion und Absorption) bis
zu sehr hohen Werten χ > 105(schallhart) konstant mit
sehr guter ¨
Ubereinstimmung.
Optimierung mit FEM Referenz
Zur weiteren Validierung werden im Folgenden die Re-
ferenzsignale mithilfe der kommerziellen Software COM-
SOL Multiphysics®erzeugt. Dabei wird eine Simulation
mit der finiten Elemente Methode (FEM) im Frequenz-
bereich durchgef¨
uhrt und die Signale an den Aufpunkten
anschließend in den Zeitbereich transformiert. Der Ab-
sorber wurde mit dem Miki-Modell [8] modelliert. Darin
wird eine verringerte, effektive Schallgeschwindigkeit in-
nerhalb des Absorbervolumens ber¨
ucksichtigt. Ansonsten
0.000 0.002 0.004 0.006 0.008 0.010
tin s
−0.02
0.00
0.02
0.04
0.06
p’ in Pa
mic 1opt
mic 2opt
mic 1ref
mic 2ref
Abb. 2: Druckverlauf ¨
uber Zeit an den Empf¨
angerpositionen
xmic1 und xmic2 (siehe Abb. 1) des Referenzsignals und des
optimierten Signals.
0.00 0.25 0.50 0.75 1.00 1.25
xin m
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
yin m
Quelle 1
Quelle 2
0.2
0.4
0.6
0.8
1.0
φ
Abb. 3: Schematischer Aufbau des Rechengebiets zur Opti-
mierung des Str¨
omungswiderstands. Die Werte von φlaufen
von dunkelblau bis gelb. Die gr¨
unen Punkte stellt die Quellpo-
sitionen dar und der (hell)blaue Bereich ist der Auswertungs-
bzw. Empf¨
angerbereich.
werden die Einstellungen des vorigen Abschnitts behal-
ten.
In Abb. 5 sind die Schalldrucksignale im Zeitbereich dar-
gestellt. In Abb. 5a wird der Defekt der Schallgeschwin-
digkeit nicht ber¨
ucksichtigt. Durch die Anpassung der
Dichte ϱ0in Abb. 5b im Gebiet des Absorbers konnte der
Schallgeschwindigkeitsdefekt kompensiert werden und ei-
ne deutliche Verbesserung der Reproduktion erreicht
werden. Zur korrekten Abbildung der Reflexion ist die
Oberfl¨
achenimpedanz des Absorbers entscheidend und
die Modellierung der Schallgeschwindigkeit im Absor-
ber nicht zwingend n¨
otig. Letztlich ist die Ver¨
anderung
der Schallgeschwindigkeit im por¨
osen Medium frequenz-
abh¨
angig und erfordert weitere Untersuchungen, um die
Transmission durch einen por¨
osen Absorber mithilfe den
hier vorgestellten Gleichungen zu modellieren.
DAGA 2023 Hamburg
1071
0.0 0.5 1.0 1.5 2.0 2.5 3.0
x2in m
0.0
0.5
1.0
1.5
2.0
2.5
3.0
x1in m
(a) Referenz des Str¨
omungswiderstand.
0.0 0.5 1.0 1.5 2.0 2.5 3.0
x2in m
0.0
0.5
1.0
1.5
2.0
2.5
3.0
x1in m
(b) Optimierter Str¨
omungswiderstand.
0
1000
2000
3000
4000
5000
χ
Abb. 4: R¨
aumliche Verteilung des Str¨
omungswiderstandes bzw. von χ.
0.002 0.004 0.006 0.008 0.010
tin s
−0.2
−0.1
0.0
0.1
0.2
0.3
p’ in Pa
mic 1opt
mic 2opt
mic 1ref
mic 2ref
(a) ϱ0,absorber =ϱ0,air.
0.002 0.004 0.006 0.008 0.010
tin s
−0.2
−0.1
0.0
0.1
0.2
0.3
p’ in Pa
mic 1opt
mic 2opt
mic 1ref
mic 2ref
(b) ϱ0,absorber = 2,4 kg/m3.
Abb. 5: Optimierung mit Comsol Referenz mit verwen-
detem Miki-Modell mit Str¨
omungswiderstand χComsol =
5000 Pa s/m2.
Optimierung des Str¨
omungswiderstands
im 2D Raum
Das in der Abb. 3 dargestellte Rechengebiet dieses
Testfalls erstreckt sich ¨
uber 1,4 m ×1,4 m. An den
R¨
andern wirkt das effektive Volumen bei φ= 10−3
¨
ahnlich zu einer schallharten Wand [3]. Im orange-
farbigen Bereich wird innerhalb der Optimierung der
Str¨
omungswiderstand χangepasst. Zur Ermittlung der
Referenz wird zun¨
achst eine Simulation im schallharten
Zustand gerechnet und die mittlere energy decay curve
(EDC) ermittelt. Diese kann anschließend durch einen
gew¨
unschten, k¨
unstlichen Abfall angepasst werden. In
diesem Fall wird der Abfall an jedem Gitterpunkt im
Empf¨
angerbereich durch
Lp,decay =[t1, t2, ..., tN]
tN
(−9) dB (10)
pmod =prigid ·10
Lp,decay
20 Pa (11)
modifiziert und dient anschließend als Optimierungsziel.
Da die lokale Anpassung von χnicht glatte Verl¨
aufe er-
zeugt (Abb. 7a), die in der Praxis nur schwer umsetz-
bar sind, wird der Einfluss von ¨
ortlich konstanten χ-
Bereichen untersucht (Abb. 7b).
Das Ergebnis der Optimierung des Darcy Term ist in
Abb. 7 f¨
ur (a) ohne und (b) mit Kassetierung gezeigt. Die
Ergebnisse der Optimierung hinsichtlich der manipulier-
ten EDC sind in Abb. 6 dargestellt. Es ist zu erkennen,
dass die mit den Gl. (10) und Gl. (11) erstellten Refe-
renzkurven durch die Optimierung sehr gut reproduziert
werden k¨
onnen.
Fazit
In diesem Beitrag wurde die prinzipielle Eignung der
Optimierung des Str¨
omungswiderstandes mittels ei-
nes adjungierten-basierten Verfahren in einem FDTD-
Framework untersucht. Weiterhin wurde die Optimie-
rung hinsichtlich eines gegebenen EDC-Verlaufs unter-
sucht. In beiden F¨
allen wurde festgestellt, dass das Ver-
fahren zur Optimierung geeignet ist und die Ergeb-
nisse der Optimierung mit den Referenzl¨
osungen sehr
gut ¨
ubereinstimmen. F¨
ur den Fall der Reflexion am
por¨
osen Absorbern wurde eine sehr gute Reproduktion
DAGA 2023 Hamburg
1072
der Referenzsignale erreicht. Es wird allerdings deutlich,
dass f¨
ur eine korrekte Beschreibung der Transmission ei-
ne Erweiterung der Modellierung notwendig ist. Auch
die Optimierung hinsichtlich raumakustischer Parameter
ist durch die erfolgreiche Reproduktion der EDC-Kurve
m¨
oglich. Es muss dabei aber der Zusammenhang zwi-
schen der EDC und den raumakustischen Parametern
beachtet werden, da raumakustische Parameter gene-
rell keine Eindeutigkeit bzgl. des zeitlichen Schalldruck-
signals liefern.
Danksagung
Die Autoren bedanken sich f¨
ur die finanzielle Un-
terst¨
utzung durch die Deutsche Forschungsgemeinschaft
(DFG) innerhalb der Projekte LE 3888/2 und WE
4057/16.
Literatur
[1] Liu, Q.; Vasilyev, O.V. (2007): “A brinkman penalization
method for compressible flows in complex geometries.” In:
Journal of Computational Physics,227(2):946–966.
[2] Reiss, J. (2022): “Pressure-tight and non-stiff volume
penalization for compressible flows.” In: Journal of Scien-
tific Computing,90(86):1–29.
[3] Lemke, M.; Reiss, J. (2023): “Approximate acoustic boun-
dary conditions in the time-domain using volume penali-
zation.” In: The Journal of the Acoustical Society of Ame-
rica,153(2):1219–1228.
0.00 0.01 0.02 0.03 0.04 0.05 0.06
tin s
−25
−20
−15
−10
−5
0
Amplitude in dB
schallhart
Referenz/Ziel
optimiert
(a) EDC der Optimierung.
0.00 0.01 0.02 0.03 0.04 0.05 0.06
tin s
−25
−20
−15
−10
−5
0
Amplitude in dB
schallhart
Referenz/Ziel
optimiert
(b) EDC der Optimierung mit 10 cm Bereichen mit kon-
stanten χ.
Abb. 6: Energy Decay Curve (EDC) der Optimierung.
[4] Bilbao, S. (2022): “Immersed boundary methods in wave-
based virtual acoustics.” In: The Journal of the Acoustical
Society of America,151(3):1627–1638.
[5] Lemke, M. (2015): Adjoint based data assimilation in com-
pressible flows with application to pressure determination
from PIV data. Ph.D. thesis, Technische Universit¨
at Ber-
lin.
[6] Stein, L.; Straube, F.; Sesterhenn, J.; Weinzierl, S.; Lem-
ke, M. (2019): “Adjoint-based optimization of sound re-
inforcement including non-uniform flow.” In: The Journal
of the Acoustical Society of America,146(3):1774–1785.
[7] H¨
olter, A.; Lemke, M.; Weinzierl, S. (2022): “Analysis and
comparison of FDTD discretization procedures for room
acoustical simulations.” In: Fortschritte der Akustik: Ta-
gungsband d. 48. DAGA, 903–906.
[8] Miki, Y. (1990): “Acoustical properties of porous
materials-modifications of delany-bazley models.” In: The
Journal of the Acoustical Society of Japan,11(1):19–24.
0.00 0.25 0.50 0.75 1.00 1.25
x2in m
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
x1in m
0
20
40
60
80
100
120
140
χ
(a) Optimierter Str¨
omungswiderstand.
0.00 0.25 0.50 0.75 1.00 1.25
x2in m
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
x1in m
0
80
90
100
110
χ
(b) Optimierter Str¨
omungswiderstand mit ¨
ortlicher Mittelung
in 10 cm Bereichen.
Abb. 7: R¨
aumliche Verteilung des Str¨
omungswiderstandes
im 2D Raum.
DAGA 2023 Hamburg
1073
