用于声学模拟的传输线矩阵
5.00/5 (27投票s)
声波传播TLM建模的实现和理论，以及2D和3D视图。还包括雨滴和船尾波模拟。
介绍
声波模拟在许多情况下是一个备受关注的话题，但它们难以在最少计算机能力的情况下实现。有几种候选方法，但正如我们将看到的，它们都有优缺点。最精确的模拟技术之一是传输线矩阵，简称TLM，因为它能够在给定正确前提条件的情况下模拟任何情况。一些最著名的声学技术包括：
射线追踪通常在音乐厅声学中大量使用，因为其主要领域（其有效频率）通常落在可以应用射线追踪的范围。一般来说，这适用于较高频率领域，其中足够大的部分可以是
所谓的有限差分时域技术（简称FDTD）实际上与声学完全相同，因为FD-TD方法求解波方程，而TLM提供了一个物理传播模型，这意味着两者的局限性和优点完全相同。需要提及的是，TLM和FDTD都可以用于声学和电动力学，因为它们都是可以通过机电波中的电压和电流，以及声学中的压力和粒子速度建模的波。
FEM方法是许多商业产品中使用的非常流行的算法，特别适用于模拟低频域的波方程，因为该方法必须根据其模拟的频率设置点。FEM技术求解亥姆霍兹波方程，其中需要首先输入要求解的频率。
BEM方法在生成有限问题的精确解方面非常流行，例如声波通过狭缝和其他开口的传播。如果几何形状变得更加多样化，它很快就会变得复杂，通常用于简单情况，本质上是一种分析技术。
本文的目标是展示波模拟可以多么简单地完成，以及设计过程实际上是多么直观。还有一些更近期的研究涉及格林-玻尔兹曼模型，这也可以用于线性波方程以及纳维-斯托克斯方程。
如果您想了解这些技术背后的更多细节，我建议您阅读我的一位数值声学教授撰写的这篇PDF文档。
背景
TLM方法的原理最初由惠更斯描述，他是一位荷兰科学家，在数学和物理学领域取得了许多发展（他还根据连分数建造了第一个天象仪），他将波传播简化为一系列球体的弹性碰撞，这实际上在物理上几乎完全正确，因为声波是借助分子碰撞传播的，分子周围环绕着电磁场（围绕原子核旋转的电子）以及分子之间的空间。
利用惠更斯原理计算波方程的想法最初由彼得·约翰斯提出，他也创造了“传输线矩阵”（简称TLM）这个名称，因为模拟电路类似于矩阵网格上的一系列相互连接的点。这个想法非常简单，接收来自其中一个节点的入射波（维基百科插图）：
第一张图片（最左边）显示了图中心节点处的入射波，从中心点南部的点向北传播。
压力随后被散射（中间的图片），但总能量保持不变（所有散射的总和仍然是1），向南的波是反射，因为它沿着一条线传播，并且在节点处，它看到3条线可以进一步传播，这导致了由波进入节点时看到的特性阻抗计算出的反射和散射系数。 用数学术语可以写成如下（再次来自维基百科），对于传输：
以及相应的反射系数：
最后一张图片只是展示了模式如何一遍又一遍地重复。这三张图片显示的是电路的传输线，而声学中的TLM模型则使用一系列管道，如下图所示，节点显示为灰色
这些信息实际上就是您在计算机上模拟简单自由空间传播所需的一切，但我们对更复杂的场景感兴趣，因此了解数学理论和其含义很重要。
声学管道需要一些额外的解释，因为它们也必须具有一定的尺寸才能正常工作。这并没有明确设置，因为在TLM模型中，假设只有平面波会在管道内从一个节点传播到下一个节点，这实际上可以用一个方程矩阵来描述：
角度通常不包含在此模型中，因为假设只有正入射平面波在管道内传播。还可以注意到，长度（或距离）包含在上述矩阵中，但未包含在模型中。为了使模型有效，需要选择足够小的距离，以便压力在管道上不会发生变化。 上述方程常用于计算吸声材料的行为，以及声波穿过墙壁等的传播。它还意味着存在一个截止频率，超过此频率，我们的模型不再可靠。
在TLM模型中，以不同的形式查看传播矩阵更有趣，因为它从一个节点到另一个节点。方程因此取决于4条可能的传播线，所以从西向东传播的波的表达式变为：
上式中描述的角度是节点之间的长度 l 和模拟波长，其关系如下：
这个原理可以用来解决计算声波穿过薄壁传播的问题，因为薄壁以产生近场而闻名，这意味着一些声波不会从墙壁传播到房间，而是垂直于墙壁传播，从而抵消了效果。这在房屋中常见的轻质墙体结构中，对于较高频率的情况最为典型。我意识到最后这三个方程是为特别感兴趣的人准备的，但如果您想在实际世界中使用它，了解这些非常有用。这仅仅表明，无论您的计算机程序有多好，结果都取决于用户所拥有的知识水平。
代码中的基本TLM传播模型
基本原理现已确立，但毕竟这是CodeProject，所以我们需要看看实际计算波传播需要存储什么。通过对模型的仔细审查，我们可以确定在每个离散时间步中，以下事件可能按此顺序发生：
- 如果需要，通过设置入射压力，为每个声源创建声压
- 计算每个节点所有入射声压的散射
- 将声压从入射节点传播到相对的出射节点
- 传入压力分量的总和给出每个节点当前的压力
这种澄清是绘制每个节点所需变量方案的全部需求，其中每个节点在每个方向上连接到一个相似的节点：
现在，这完全是计算节点的输入和输出数量的问题，每个四个方向都有两个，以及当前压力。缩写词表示连接的位置以及它是向内还是向外传播。总而言之，对于一个没有损耗的简单传播模型，我们需要9个矩阵来存储所有变量；像这样：
Private Dimensions As Integer = 100
Private _
IncommingEast(Dimensions, Dimensions), _
InncommingNorth(Dimensions, Dimensions), _
IncommingWest(Dimensions, Dimensions), _
IncommingSouth(Dimensions, Dimensions), _
ScatteredEast(Dimensions, Dimensions), _
ScatteredNorth(Dimensions, Dimensions), _
ScatteredWest(Dimensions, Dimensions), _
ScatteredSouth(Dimensions, Dimensions) _
As Decimal
Private CurrentPressure(Dimensions, Dimensions) As Decimal
该算法是迭代的，因为需要遍历矩阵中的所有点，我首先在给定节点设置压力，然后计算散射压力。重要的是要在添加当前压力之前，在一个单独的for循环中计算所有散射脉冲。这实际上是符合逻辑的，因为所有元素的散射必须在您能够判断下一个压力是什么之前发生。
For x As Integer = 0 To Dimensions - 1
For y As Integer = 0 To Dimensions - 1
If Source.Contains(New Point(x, y)) Then
IncommingEast(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + IncommingEast(x, y)
IncommingWest(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + IncommingWest(x, y)
IncommingSouth(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + IncommingSouth(x, y)
InncommingNorth(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + InncommingNorth(x, y)
End If
ScatteredNorth(x, y) = 0.5 * (+IncommingEast(x, y) - InncommingNorth(x, y) + _
IncommingWest(x, y) + IncommingSouth(x, y))
ScatteredEast(x, y) = 0.5 * (-IncommingEast(x, y) + InncommingNorth(x, y) + _
IncommingWest(x, y) + IncommingSouth(x, y))
ScatteredWest(x, y) = 0.5 * (+IncommingEast(x, y) + InncommingNorth(x, y) - _
IncommingWest(x, y) + IncommingSouth(x, y))
ScatteredSouth(x, y) = 0.5 * (+IncommingEast(x, y) + InncommingNorth(x, y) + _
IncommingWest(x, y) - IncommingSouth(x, y))
Next
Next
散射压力是向外的，我需要考虑到每个节点可以从所有4个不同的连接点（北、西、东、南）接收传入压力：
For x As Integer = 0 To Dimensions - 1
For y As Integer = 0 To Dimensions - 1
'Free field propagation, no obstacles in the current cell
'and all surrounding cells gives a scattered field
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
CurrentPressure(x, y) = (0.5 * (ScatteredEast(x, y) + _
ScatteredNorth(x, y) + ScatteredWest(x, y) + ScatteredSouth(x, y)))
Next
Next
这两个循环都在每个离散时间步中发生，但如果没有完美匹配层（PML），此代码将无法正常工作，因为我们尚未添加末端修正。上面的代码只有在实际点四周都被空气（或所有点中的相同物体）包围时才能工作。现在是引入边界和具有反射系数的物体的时候了。
网格的反射和终止
反射已经使用并描述过，但在那种情况下相当简单，因为所有的声学管都是相同的。反射和传输实际上只是说明当波到达节点时，它看到多少阻力。如果阻力高，它将主要被反射；如果阻力低，它将主要被传输。
在声学中，反射由介质的特征阻抗和两个系数之间的差异决定。反射系数由关系式定义
其中Za是墙壁或边界的特征阻抗，而Zt是TLM模型中自由场的特征阻抗。Zt值实际上不是空气中使用的正常值，其原因实际上是网格的离散化，因为它被平方了。这具有深远的影响，因为2D中传播速度（在空气中定义为C0，值为343 m/s）实际上是c0 / Sqrt(2)（在3D中是c0 / Sqrt(3)）。这只是对角传播的勾股定理。
要在TLM结构中使用此信息，我们需要通过此方法计算反射
这有一些有趣的效果，例如；如果你想让你的网格具有无限大TLM矩阵的边界条件，你需要将Za与Zt匹配（相等的值，通常称为PML），这样你就会得到以下非反射边界条件：
直观上这实际上是合理的，因为人们会期望它不为零，因为这将导致无限吸收，而是粒子实际上会碰撞，并且一些能量会传回网格。
使PML反射可行的代码如下所示：
For x As Integer = 0 To Dimensions - 1
For y As Integer = 0 To Dimensions - 1
If x = 0 And y = 0 Then
'Compared to the free field propagation x-1 and y-1 is illegal
'This is the top left corner
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncommingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = 0 And y = Dimensions - 1 Then
'Compared to the free field propagation x-1 and y+1 is illegal
'This i sthe top right side
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncommingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf x = Dimensions - 1 And y = 0 Then
'Compared to the free field propagation x+1 and y-1 is illegal
'This is the bottom left corner
IncommingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = Dimensions - 1 And y = Dimensions - 1 Then
'Compared to the free field propagation x+1 and y+1 is illegal
'This is the Bottom right corner
IncommingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncommingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf x = 0 Then
'Compared to the free field propagation x-1 is illegal
'This is the top line
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf y = 0 Then
'Compared to the free field propagation y-1 is illegal
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf y = Dimensions - 1 Then
'Compared to the free field propagation y+1 is illegal
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncommingWest(x, y) = ScatteredWest(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf x = Dimensions - 1 Then
'Compared to the free field propagation x+1 is illegal
IncommingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
'Free field propagation, no obstacles in the current cell
'and all surrounding cells gives a scattered field
If Wall.Contains(New Point(x + 1, y)) Then
IncommingEast(x, y) = 0.5 * ScatteredEast(x, y)
Else
IncommingEast(x, y) = ScatteredWest(x + 1, y)
End If
If Wall.Contains(New Point(x, y + 1)) Then
InncommingNorth(x, y) = 0.5 * ScatteredNorth(x, y)
Else
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
End If
If Wall.Contains(New Point(x - 1, y)) Then
IncommingWest(x, y) = 0.5 * ScatteredWest(x, y)
Else
IncommingWest(x, y) = ScatteredEast(x - 1, y)
End If
If Wall.Contains(New Point(x, y - 1)) Then
IncommingSouth(x, y) = 0.5 * ScatteredSouth(x, y)
Else
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
End If
End If
CurrentPressure(x, y) = (0.5 * (ScatteredEast(x, y) + _
ScatteredNorth(x, y) + ScatteredWest(x, y) + ScatteredSouth(x, y)))
Next
Next
TLM 是一个非常直观的模型，人们可以很容易地理解实际发生的事情，而无需付出太多努力（至少在简单情况下）。 如果以这里为例
InncommingNorth(x, y) = 0.5 * ScatteredNorth(x, y)
我只是将所有墙壁（非反射边界除外）设置为0.5，我只是取向北散射的能量，并将其以循环方式发送回传入的北方，只是它被反射系数衰减了。在TLM节点的图片中，它对应于通过电阻器将模型在北端短路。侧壁设置为非反射边界条件（PML的实际定义），其中R的实际值约为-0.17。
文章中显示的代码片段在实际程序中有所改动，因为有更快的方法来检查它是否是自由场节点。想法是设置一个与节点大小相等的布尔数组。在这里，我们指定任何节点是否有一堵墙，或者在北、南、东或西方向上有一堵相邻的墙。这是在添加墙壁元素时完成的
Public Sub SetWalls(ByVal WallPoints As List(Of Point), ByVal ReflectionFactor As Decimal)
Dim R As Decimal = ((1 + ReflectionFactor) - SquareRoot(2) * (1 - ReflectionFactor)) / _
((1 + ReflectionFactor) + SquareRoot(2) * (1 - ReflectionFactor))
For Each p As Point In WallPoints
Walls(p.X, p.Y) = 0.5
WallNeigbor(p.X, p.Y) = True
If p.X <> 0 Then
WallNeigbor(p.X - 1, p.Y) = True
End If
If p.X <> Dimensions - 1 Then
WallNeigbor(p.X + 1, p.Y) = True
End If
If p.Y <> 0 Then
WallNeigbor(p.X, p.Y - 1) = True
End If
If p.Y <> 0 Then
WallNeigbor(p.X, p.Y + 1) = True
End If
Next
End Sub
该代码还存储了每个墙元素在添加时指定的反射系数。这也存储在一个数组中，并且可以在每次遇到墙元素时调用。
现在可以通过简单地检查Wall Neighbor是否为真来确定它是否是自由传播的场
If Not (WallNeigbor(x, y)) Then
'Its a free field computation
IncomingEast(x, y) = ScatteredWest(x + 1, y)
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
IncomingWest(x, y) = ScatteredEast(x - 1, y)
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
If WallNeigbor(x + 1, y) Then
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
Else
IncomingEast(x, y) = ScatteredWest(x + 1, y)
End If
If WallNeigbor(x, y + 1) Then
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
Else
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
End If
If WallNeigbor(x - 1, y) Then
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
Else
IncomingWest(x, y) = ScatteredEast(x - 1, y)
End If
If WallNeigbor(x, y - 1) Then
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
Else
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
End If
End If
Else
'Its on the borders
If x = 0 And y = 0 Then
'Compared to the free field propagation x-1 and y-1 is illegal
'This is the top left corner
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
IncomingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncomingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = 0 And y = Dimensions - 1 Then
'Compared to the free field propagation x-1 and y+1 is illegal
'This i sthe top right side
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
InncomingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncomingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
ElseIf x = Dimensions - 1 And y = 0 Then
'Compared to the free field propagation x+1 and y-1 is illegal
'This is the bottom left corner
IncomingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
IncomingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = Dimensions - 1 And y = Dimensions - 1 Then
'Compared to the free field propagation x+1 and y+1 is illegal
'This is the Bottom right corner
IncomingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncomingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
If Walls(x - 1, y) = 0 Then
IncomingWest(x - 1, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x - 1, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
ElseIf x = 0 Then
'Compared to the free field propagation x-1 is illegal
'This is the top line
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
IncomingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
ElseIf y = 0 Then
'Compared to the free field propagation y-1 is illegal
IncomingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
ElseIf y = Dimensions - 1 Then
'Compared to the free field propagation y+1 is illegal
InncomingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
ElseIf x = Dimensions - 1 Then
'Compared to the free field propagation x+1 is illegal
IncomingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
End If
End If
由于假设大多数节点都在自由场中，因此这是一个相当大的速度优势，尽管因此整体代码变得复杂得多。
关于端节点也有一些特殊之处，因为墙也可以出现在模型的末端，这意味着顶部的节点，在北方没有反射，可以在东方或西方（或两者兼有）具有可能的墙元素的反射。
我应该提及，这种模拟每个反射表面的特定方式不允许它们在不同频率下具有不同的吸收和反射。如果要使用脉冲来查找两点之间的传递函数，则必须一次完成一个频率，并组合结果。还有另一种可能性，即使用FIR滤波器作为吸收和传输系数，因为现实生活中普通的墙壁和其他边界通常是低通滤波器，但这需要一些DSP（数字信号处理）知识，并且还需要知道每种材料的吸收系数。然而，这远远超出了解释简单TLM模型如何实现的范畴！
一种不同类型的墙壁
细心的读者可能已经注意到，目前使用的墙壁不能将波穿过墙壁，未被反射的能量简单地消失了，这当然不能在现实世界中发生。然而，如果您想在TLM模型中使用均匀网格，在实现这样的墙壁时存在一个问题，但幸运的是，我们可以通过两种辅助方法来解决这个问题。
我们不能在TLM模型中增加声速，使其超过原始值，但我们可以减小它，从而产生反射和传输，例如两种气体。 实现这一点的通常方法是在每个节点中使用两个特殊分支，它们同时处理声速变化和声波衰减。
声速的可能变化也使得非线性传播的模拟成为可能，其中速度根据节点处的当前压力而变化。当然，这不包含在视图中，因为它是一个高度专业化的任务，但如果设置得当，TLMGrid中的代码是完全可用的。
我们还可以根据两种不同材料的重量设置反射和透射，如果传播波在从一个节点到下一个节点时看到两种不同的重量，这将产生反射和透射，这也会在网格中导致真实的反射和透射情况。但是，请不要尝试将这种类型的墙壁与全吸收墙壁结合使用，因为这会完全混淆模型，结果会变得一团糟。一些专家在此PowerPoint演示文稿中进行了快速介绍（Pdf文档）。
最后，我应该提到，TLM模型目前只适用于纵波，并且由薄壁（具有惯性矩）产生的波不能在此模拟中结合。您将不得不通过一些花哨的FIR滤波器方法或通过在FDTD模型中求解方程来解决这个问题。
将模型中的频率与现实世界关联起来
将TLM模型与现实世界关联起来存在许多问题，正如您已经看到的模型中的声速，在每个方向上都是C0/Sqrt(2)。因此，要求真实速度C0出现在一个等边三角形的斜边上，传播角度为45度，如下图所示
45度角的原因是这是模型正确处理的最差角度，因为它必须在x方向和y方向上移动才能到达顶点。如果将模型设置为将任何非0的绝对压力绘制成红色，您将看到模型在直线x和y方向上传播得相当快，但圆形45度半径长度之外的点的振幅非常小。波浪将显示为源点周围的一个正方形。
有两种类似的方式来看待这一点，并得出关于它们如何关联的结论。一种看待方式是通过推理，压力需要两个离散时间脉冲才能到达45度角处的点，因为它必须先沿x轴移动，然后再沿y轴移动（反之亦然），这意味着它将以Ct的速度移动2 Delta L。
将传播介质中驻波的波长或频率（对应于介质的波长）关联起来的标准公式是：
然而，TLM模型中的正弦源只需要时间步长t作为输入（但也可以通过恒定相移进行调整，这样你就可以将正弦波的起始延迟一个固定的量。这将允许你使用两个不同相位的源来模拟偶极子源）。因此，给定频率的生成正弦函数如下所示：
为了避免混淆，我必须强调是t减半而不是l，因为在上面的方程中写的方式可能会造成混淆。结果是当它放入模型时，频率实际上减半了。最后，TLM模型中实际正弦源的代码：
If Sources(x, y) Then
Dim Sine As Decimal = Amplitude * Math.Sin(Math.PI * CurrentTimeStep * _
DeltaLength * SourceFrequency / CT)
IncomingEast(x, y) = Sine + IncomingEast(x, y)
IncomingWest(x, y) = Sine + IncomingWest(x, y)
IncomingSouth(x, y) = Sine + IncomingSouth(x, y)
InncomingNorth(x, y) = Sine + InncomingNorth(x, y)
End If
除了正弦函数之外，还有其他生成源的方法。最常用，也许最有价值的是狄拉克脉冲。从理论上讲，这是一个在无限小时间间隔内的无限高点源，在数学理论中，这意味着源处具有完全平坦的频率响应，或者换句话说，如果对源信号进行傅里叶变换，所有频率都具有相同的振幅。这编码起来极其简单，因为当时间步长等于t=0时，它将只有一个振幅。
If DeltaSources(x, y) Then
IncomingEast(x, y) = Amplitude + IncomingEast(x, y)
IncomingWest(x, y) = Amplitude + IncomingWest(x, y)
IncomingSouth(x, y) = Amplitude + IncomingSouth(x, y)
InncomingNorth(x, y) = Amplitude + InncomingNorth(x, y)
End If
最后我提到了一些关于截止频率的问题，超过这个限制，TLM就不再有效了。截止频率定义为dl/Lambda = 0.25，但推荐的最大值是dl/Lambda = 0.1，如果您感兴趣的模拟频率高于此值，您不能假定结果会准确。然而，即使您超过此限制，由于抗混叠效应，您仍有可能获得稳定的输出。
这导致在一个完整的周期内对正弦波进行10个或更多离散采样，其思想是正弦波的每个步骤都可以很好地建模，使其近似为一条直线。顺便说一下，这并非TLM模型所独有，因为FDTD和TLM都依赖于仅使用泰勒级数的一阶导数，如果我们需要模拟非线性效应，我们还需要进行调整，以便每个离散样本之间的差异可以由一条直线精确建模。
截止频率实际上是不使用狄拉克δ脉冲的一个原因，因为它必然会包含高于截止频率的频率分量，从而破坏测量点处良好的传递函数效应，并且传播在模型中也会变得“粗糙”，因为它似乎在图中包含了一些随机噪声。因此，使用所谓的高斯脉冲作为源而不是狄拉克δ脉冲更为方便，其定义为：
当a趋于0时，此表达式将生成狄拉克δ脉冲。请记住，我们不希望有任何频率分量高于TLM截止频率，因此我们还需要知道截止频率（在频率频谱中幅度减半的地方），这可以通过高斯滤波器的理论轻松找到，并且可以在这里找到。然而，高斯源有点麻烦，因为频率响应退化速度相当慢。实际上，这意味着脉冲对于非常高的频率会变得极其长，因此不同类型的源可能更好。使用巴特沃斯或切比雪夫滤波器特性的源可能是更好的方法，但在本文中我将忽略这一点。
构建3D查看器并提高速度
3D 查看器的设计深受 Sean Sexton 的雨滴动画（你可以在这里下载他的原始 C# 代码）的影响。然而，他做了一些简化，因为他不需要担心墙壁和其他物体会阻碍波浪传播。他的第一个实现直接使用了狄拉克δ脉冲作为源，而他的第二个版本使用了近似于简单高斯滤波器的平滑滤波器。（高斯滤波器基本上是一个去除高频分量的低通滤波器。） 通过直接使用高斯源，可以获得更平滑的表面
我的实现中实际上并没有剩下他代码的太多内容，除了这篇MSDN文章中提供的一些技巧，它基本上阐述了如何在WPF中以最快的方式初始化和更新3D模型。基本思想是只初始化所有3D数组一次，并将一个3D点列表连接到视图模型，当您更新计算时，您应该更新当前未链接到视图模型的数组，然后交换这两个数组。
Dim RunTimeInMS As Double = 1
Dim LastTimeRendered As Double = 0
Dim IsSimulationRunning As Boolean = False
Dim RunSimulationEvent As New EventHandler(AddressOf RunSimulation)
Private Sub btnStart_Click(sender As System.Object, _
e As System.Windows.RoutedEventArgs) Handles btnStart.Click
If IsSimulationRunning Then
RemoveHandler CompositionTarget.Rendering, RunSimulationEvent
btnStart.Content = "Resume"
IsSimulationRunning = False
Else
AddHandler CompositionTarget.Rendering, RunSimulationEvent
btnStart.Content = "Pause"
IsSimulationRunning = True
End If
End Sub
有一点有点令人费解，那就是您在VB中需要EventHandler，但在C#中不需要。如果您像下面这样编写代码，您会收到一个不会消失的警告。我猜这与VB不能在没有参数的情况下提供委托，而C#可以有关。
AddHandler CompositionTarget.Rendering, AddressOf RunSimulation
在VB中会收到一个意义不大的警告，通过EventHandler委托来解决。现在可以通过如下编写子程序RunSimulation在特定的时间间隔设置渲染
Sub RunSimulation(ByVal sender As Object, ByVal e As RenderingEventArgs)
Dim rargs = DirectCast(e, RenderingEventArgs)
If ((rargs.RenderingTime.TotalMilliseconds - LastTimeRendered) > RunTimeInMS) Then
...
...
...
LastTimeRendered = rargs.RenderingTime.TotalMilliseconds
End If
End Sub
但是，如果您只想尽快运行它，您可以删除上一节中的几乎所有内容（这就是我在代码中所做的，因为我没有以最快速度绘制它的问题）
Sub RunSimulation()
...
...
...
End Sub
也可以使用DispatcherTimer，但它可能会稍微滞后一些，但那样它也会按顺序更新。
TLM方法也是并行处理的完美候选，因为它在一个时间步内对所有节点循环两次。将其更改为并行实际上非常容易，您只需阅读此MSDN文档即可。显然，如果您的计算机只有一个处理器，它不会更快，但我有6个，它大致将速度提高到以前代码的两倍（其他处理器甚至可能更快，例如8个CPU）。其思想是更改类TLMGrid中名为CalculateNextTimeStep的函数中的代码。没有并行处理的清除函数如下所示：
Public Function CalculateNextTimeStep() As Decimal(,)
For x as Integer = 0 to Dimensions - 1
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
Next
For x as Integer = 0 to Dimensions - 1
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
Next
CurrentTimeStep += 1
Return CurrentPressure
End Function
要并行实现代码，首先添加Imports System.Threading.Tasks，然后简单地将两个for循环更改为：
Public Function CalculateNextTimeStep() As Decimal(,)
Parallel.For(0, Dimensions, Sub(x)
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
End Sub)
Parallel.For(0, Dimensions, Sub(x)
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
End Sub)
CurrentTimeStep += 1
Return CurrentPressure
End Function
就是这样，您应该注意到在原始的For循环中，它一直到Dimensions - 1，而在Parallel.For中，它一直到x >= Dimensions，所以您不包括-1。遗憾的是，我无法以这种方式处理3D查看器的代码，因为其他进程占用了所使用的变量。
据说还有一件事可以改善 WPF 的性能，那就是降低屏幕更新频率，称为帧率。这在应用程序启动时设置，如下所示
Imports System.Windows.Media.Animation
Class Application
' Application-level events, such as Startup, Exit, and DispatcherUnhandledException
' can be handled in this file.
Protected Overrides Sub OnStartup(e As System.Windows.StartupEventArgs)
MyBase.OnStartup(e)
Timeline.DesiredFrameRateProperty.OverrideMetadata(GetType(Timeline), _
New FrameworkPropertyMetadata() With {.DefaultValue = 60})
End Sub
End Class
例如，将其从60降低到20，应该可以改善CPU使用率。我没有看到任何改进，所以将其保留为默认的60帧率，但您可以尝试将其降低以查看会发生什么。我怀疑这适用于故事线动画，但我在这里不确切知道会发生什么以及它是如何工作的。
有了这个新获得的速度，现在可以开始尝试不同类型的波浪，即船尾波（图片来自维基百科）
这可以用切伦科夫辐射来建模，正如谢尔盖·克留科夫在这个问题中首次向我描述的。 下面的动画图片展示了原理
它只是一系列比波传播速度更快的脉冲。这在YouTube的视频中也解释得很好。这意味着子弹和船波在数学上是相当相似的。程序中使用一系列高斯脉冲产生的动画如下所示
使用TLM模型创建它真是令人惊叹的简单。
在3D世界中移动
还有一点是关于我如何使用鼠标实现非常基本的相机视角方向更改。具体做法如下：
Dim StartPos, EndPos As New Point
Dim OriginalLookdirection As New Vector3D
Dim look As Boolean = False
Private Sub mainViewport_MouseDown(sender As System.Object, _
e As System.Windows.Input.MouseButtonEventArgs) Handles Me.MouseDown
StartPos = e.GetPosition(viewport3D1)
OriginalLookdirection = camMain.LookDirection
look = True
End Sub
Private Sub mainViewport_MouseUp(sender As System.Object, _
e As System.Windows.Input.MouseButtonEventArgs) Handles Me.MouseUp
look = False
End Sub
Private Sub mainViewport_MouseMove(sender As System.Object, _
e As System.Windows.Input.MouseEventArgs) Handles Me.MouseMove
If look Then
EndPos = e.GetPosition(viewport3D1)
Dim v As New Vector3D
v = OriginalLookdirection
v.X -= (StartPos.X - EndPos.X) * 0.1
v.Y += (StartPos.Y - EndPos.Y) * 0.1
camMain.LookDirection = v
End If
End Sub
正如您所看到的，我只是改变了X和Y的位置，同时对视角方向的X和Y进行了固定量的改变。这似乎提供了一种非常自然的方式来调整视角方向。
通过一些向量微积分的知识，您甚至可以相当容易地移动相机位置
Private Sub mainViewport_KeyDown(sender As System.Object, _
e As System.Windows.Input.KeyEventArgs) Handles Me.KeyDown
If e.Key = Key.Up Then
Dim p As New Point3D
p = camMain.LookDirection
p = Vector3D.Multiply(p, 0.01)
camMain.Position += p
End If
If e.Key = Key.Down Then
Dim p As New Point3D
p = camMain.LookDirection
p = Vector3D.Multiply(p, 0.01)
camMain.Position -= p
End If
If e.Key = Key.Left Then
Dim r As New Vector3D
r = Vector3D.CrossProduct(Normalize(camMain.LookDirection), camMain.UpDirection)
r *= 0.2
camMain.Position -= r
End If
If e.Key = Key.Right Then
Dim r As New Vector3D
r = Vector3D.CrossProduct(Normalize(camMain.LookDirection), camMain.UpDirection)
r *= 0.2
camMain.Position += r
End If
End Sub
Private Function Normalize(ByVal v As Vector3D) As Vector3D
Dim result As New Vector3D
Dim length As Double = Math.Sqrt(v.X ^ 2 + v.Y ^ 2 + v.Z ^ 2)
result = v
result /= length
Return result
End Function
然而，这会稍微搞乱简单的视角方向，因为视角方向的向量被改变了。如果你想让它正常工作，你必须调整向量的长度。
然后是缩放。
Private zoomDelta As Vector3D
Private Sub Window_MouseWheel(sender As Object, e As MouseWheelEventArgs)
If e.Delta > 0 Then
' Zoom in
camMain.Position = Point3D.Add(camMain.Position, zoomDelta)
Else
' Zoom out
camMain.Position = Point3D.Subtract(camMain.Position, zoomDelta)
End If
End Sub
其中ZoomDelta取决于视角方向，并在Window_Loaded子程序中设置。
Private Sub Window_Loaded(sender As System.Object, e As System.Windows.RoutedEventArgs) Handles MyBase.Loaded
' On each WheelMouse change, we zoom in/out a particular % of the original distance
Const ZoomPctEachWheelChange As Double = 0.02
zoomDelta = Vector3D.Multiply(ZoomPctEachWheelChange, camMain.LookDirection)
End Sub
这是一个非常快速的概述，说明了我在模拟运行时如何实现移动。尽管在模拟运行时它相当滞后，但只要视角方向的向量尺寸合理，它的功能就相当不错。我几乎可以肯定有更好的方法可以做到这一点，但目前对我来说已经足够了。
从矢量到像素
使用画布和背景图像非常好，因为我可以将输入、源墙放置与实际计算算法分离，从而形成两个独立的层。这尤其好，因为您无需担心每个像素的句柄，而是将所需数据存储在数组中。
从画布上的矢量放置到栅格的转换对于一个点来说是相当直接的，因为您可以简单地将点按宽度和高度进行缩放，然后乘以矩阵的栅格长度。将画布上的点转换为栅格中单元格的代码是相当直接的
Dim PointOnCanvas As New Point
PointOnCanvas = e.GetPosition(cnvTLM)
Dim RasterPointInMatrix As New Point
RasterPointInMatrix.Y = CInt(Math.Truncate(PointOnCanvas.X * Dimensions / cnvTLM.Width))
RasterPointInMatrix.X = CInt(Math.Truncate(PointOnCanvas.Y * Dimensions / cnvTLM.Height))
哦，X和Y的值被调换了，因为左上角点是0,0而不是0, Height。当我们要定义两个点，并希望它们在矩阵中（或在位图图像上）生成一条连续的线时，麻烦就开始了，人们通常会求助于Bresenham的画线算法。还有其他可用的算法，但它被广泛认为是标准实现。
鉴于它的广泛使用，人们可能会认为在网上很容易找到一个可用的示例，但我可能错了，因为我找不到一个在所有情况下都能正常工作的现成代码。（一些VB中无法工作的示例可以在这里和这里找到）。我没有再继续寻找，而是自己实现了一个简化版的Bresenham算法，它并不完美，因为我没有处理舍入误差控制，但修复它应该相当容易！
''' <summary>
''' Simplyfied Bresenham algorithm
''' </summary>
''' <param name="P1"></param>
''' <param name="P2"></param>
''' <returns>A list of continuos points on raster data</returns>
''' <remarks></remarks> <remarks />
Public Function GetBresenham_ishLine(ByVal P1 As Point, ByVal P2 As Point) As List(Of Point)
' Get the numbers of steps in each direction
Dim DeltaX, DeltaY As Integer
DeltaX = Math.Abs(CInt(P1.X) - CInt(P2.X))
DeltaY = Math.Abs(CInt(P1.Y) - CInt(P2.Y))
' Figure out if its a positive or negative step
Dim PositiveX, PositiveY As Double
If P1.X > P2.X Then
PositiveX = -1
Else
PositiveX = 1
End If
If P1.Y > P2.Y Then
PositiveY = -1
Else
PositiveY = 1
End If
Dim LineList As New List(Of Point)
'Figure out which one has the most steps and use it
If DeltaX > DeltaY Then
For i As Integer = 0 To DeltaX
Dim TempX, TempY As Integer
TempX = P1.X + PositiveX * i
TempY = CInt(P1.Y + DeltaY / DeltaX * PositiveY * i)
LineList.Add(New Point(TempX, TempY))
Next
Else
For i As Integer = 0 To DeltaY
Dim TempX, TempY As Integer
TempX = CInt(P1.X + DeltaX / DeltaY * PositiveX * i)
TempY = CInt(P1.Y + PositiveY * i)
LineList.Add(New Point(TempX, TempY))
Next
End If
Return LineList
End Function
对于这个程序的目的，它工作得足够好，因为它只是为了展示TLM代码的工作原理。有了这个简单的事情，你就可以观察不同设置的行为。
问题
这实际上是您生成无损TLM代码所需的全部内容，这意味着压力随距离下降的唯一原因是能量分散到更大的区域。这实际上是一种相当正常的看待短距离声传播的方式，然而，在多次反射或长距离传播时，这可能无法正确反映现实。
TLM 技术具有一些特性，这些特性对于抛物线方程而言优于 FDTD，因为 FDTD 在中心差分处不稳定，并且仅在后向差分处条件稳定，而 TLM 是无条件稳定的。然而，TLM 确实存在不同的问题。
- 点源发出的波不会以均匀模式扩散，因为它是一种离散技术
- 传播速度取决于
1/Sqrt(2)
所提供的代码不包括传播损耗，为此需要另外两个矩阵，因为每个节点都会有一个非反射管，消耗能量。还需要实现一个简单的高斯滤波器来改进结果，以便进行任何专业用途。
其他问题是这是一个时间相关的算法，因此为了利用结果，通常需要在接收端执行傅里叶变换，而这并未包含在现有代码中。
创建水效应的更快方法
这里使用的TLM方法旨在提供对波的精确数学描述，如果您只对创建雨滴和船尾波等涟漪效应感兴趣，则有更快的方法可用。
肖恩提供的一个链接中很好地解释了这种方法，其目标只是创建一个看起来足够好，而不是精确数学准确的东西。主要思想是使用当前单元格周围单元格的区域平均技术。
历史
本文的主要目标是帮助您理解如何合理地编写TLM算法，并为您提供一些关于如何在WPF引擎上实现真实世界物理模拟的想法。
高斯脉冲仍需改进，目前它不能很好地对应变化的频率，因为截止频率到滤波器的实现不正确。还有一个问题是如何在系统中实现阻尼，在声学中表示为奈培系数的衰减。进一步开发不仅能阻尼波浪，还能实现传输的元素也很有用，但这需要在墙壁部分进行一些重写，因为我们将需要为每面墙壁提供4个数组，而不是当前只有一个。在这种2D模拟使用的均匀网格中，声速也无法改变，因此除了如何在VB/C#中使用WPF相对快速地实现TLM方法之外，它的应用有限。
如果您想将3D TLM模型的结果用于创建雨滴效果以外的用途，那么它会更有趣，但同样，要了解发生了什么并不那么容易，因为您需要查看每个横截面。不过，该模型并不难扩展，您只需实现上下方向，以及两个相应的传入和传出压力支路。然后应该提到的是，构建网格有比正方形更好的方法。
参考文献
- 传输线矩阵（TLM）-扩散应用技术，作者：Donard de Cogan
- 声波传播的离散惠更斯模型方法，作者：Y. Kagawa et. al。
- 传输线矩阵 – 多极展开方法，一份博士论文的pdf文档，作者：Petr Lorenz
- 波和散射方法在数值模拟中的应用，作者：Stefan Bilbao
- 一个简单的TLM Matlab程序，由NTNU, Trondheim创建。