清华大学计算固体力学第七次课件ALE公式市公开课一等奖百校联赛特等奖课件.pptx

非线性有限元非线性有限元第第7 7章章 任意任意LagrangianLagrangian和和EulerianEulerian公式公式 计算固体力学计算固体力学第1页第第7 7章章 任意任意LagrangianLagrangian和和EulerianEulerian公式公式 1 1引言引言2 2ALEALE连续介质力学连续介质力学3 3ALEALE守恒规则守恒规则4 4ALEALE控制方程控制方程5 5弱形式弱形式6 6网格更新算法网格更新算法7 7Petrov-GalerkinPetrov-Galerkin方法方法第2页1 引言引言处处理理：在在发发生生严严重重大大变变形形模模拟拟中中，重重新新划划分分网网格格是是不不可可防止，工作量大，而且因为网格投影引入了误差。防止，工作量大，而且因为网格投影引入了误差。提出提出：许多问题应用：许多问题应用Lagrangian网格不能有效地处理。网格不能有效地处理。问问题题：当当材材料料严严重重变变形形时时，LagrangianLagrangian单单元元一一样样发发生生严严重重扭扭曲曲，因因为为它它们们随随材材料料一一起起变变形形，从从而而恶恶化化了了这这些些单单元元近近似似精精度度，尤尤其其是是对对于于高高阶阶单单元元。所所以以，在在积积分分点点JacobianJacobian行行列列式式可可能能成成为为负负值值，从从而而使使计计算算中中止止或或者者引引发发严严重重局局部部误误差差。另另外外，也也恶恶化化了了线线性性化化牛牛顿顿方方程程条条件件，而且显式稳定时间步长显著地下降。而且显式稳定时间步长显著地下降。第3页 一一个个Lagrangian网网格格像像在在材材料料上上蚀蚀刻刻：当当材材料料变变形形时时，蚀蚀刻刻(和和单单元元)伴伴随随变形。变形。一一个个Eulerian网网格格像像放放在在材材料料前前面面一一薄薄片片玻玻璃璃上上蚀蚀刻刻：当当材材料料变变形形时时，蚀蚀刻刻不不变变形形，而而材材料料横穿过网格。横穿过网格。1 引言引言第4页Lagrangian网格，材料点与网网格，材料点与网格点保持重合，单元随材料变格点保持重合，单元随材料变形，适合描述固体与结构变形，形，适合描述固体与结构变形，但轻易严重扭曲。但轻易严重扭曲。处理方法：处理方法：ALE网格网格(Arbitrary Lagrangian Eulerian)节点能够有序地任意运动，在节点能够有序地任意运动，在边界上节点保持在边界上运动，边界上节点保持在边界上运动，内部节点运动使网格扭曲最小内部节点运动使网格扭曲最小化。化。1 引言引言第5页1 引言引言Mesh adaptivity is based on solution variables as well as minimum element distortion Elements concentrate in areas where they are needed Adaptation is based on boundary curvature Deformation of a rubber sealInitial configurationUniform adaptivitySolution-dependent adaptivity第6页1 引言引言 在在一一些些问问题题中中，LagrangianLagrangian方方法法是是根根本本不不适适用用。比比如如，对对于于高高速速流流动动流流体体力力学学问问题题，如如围围绕绕机机翼翼区区域域，喷射等。喷射等。在在EulerianEulerian有有限限元元中中，网网格格与与物物质质是是相相互互独独立立，网网格格在在空空间间上上是是固固定定，材材料料从从网网格格中中流流过过。这这么么EulerianEulerian有有限限元元不不会会伴伴随随材材料料运运动动而而扭扭曲曲；不不过过，因因为材料经过单元对流，本构方程处理和更新是复杂。为材料经过单元对流，本构方程处理和更新是复杂。应用应用EulerianEulerian单元单元处理移动边界和相互作用问题处理移动边界和相互作用问题是困难，所以，发展了是困难，所以，发展了ALEALE。第7页2 ALE连续介质力学连续介质力学材料坐标与空间坐标材料坐标与空间坐标 空间坐标与空间坐标与ALE坐标坐标 在在Lagrangian、Eulerian和和ALE域之间映射域之间映射 ALE坐标（参考）坐标（参考）ALE坐标与材料坐标坐标与材料坐标 相对运动关系相对运动关系第8页2 ALE连续介质力学连续介质力学在在ALE算法中，网格运动是预先设置或者是由计算得到。算法中，网格运动是预先设置或者是由计算得到。网格位移网格位移 网格速度网格速度 网格加速度网格加速度 ALE ALE网格加速度和速度没有任何物理意义。当网格是网格加速度和速度没有任何物理意义。当网格是Lagrangian Lagrangian 时，它们对应于材料速度和加速度。时，它们对应于材料速度和加速度。定义定义传递速度传递速度 c，作为材料速度和网格速度之间差，作为材料速度和网格速度之间差 c0，为，为L格式；格式；cv，()为为E格式。格式。第9页2 ALE连续介质力学连续介质力学考虑一个指定函数考虑一个指定函数 为为ALE坐标坐标 和时间和时间t 函数函数 参考质点速度参考质点速度w 材料速度和网格速度差材料速度和网格速度差 对于材料速度对于材料速度第10页2 ALE连续介质力学连续介质力学利用空间梯度建立材料时间导数表示式利用空间梯度建立材料时间导数表示式 代入代入f 若代表是速度，上式为加速度若代表是速度，上式为加速度 坐标之间转换关系见例坐标之间转换关系见例7.1。(7.2.17)第11页3 ALE守恒规则守恒规则 守守恒恒规规则则，在在形形式式上上与与在在第第3 3章章EulerianEulerian描描述述中中那那些些几几乎乎相相同同，唯唯一一修修改改是是用用材材料料时时间间导导数数ALEALE形形式式(7.2.17)(7.2.17)代代替替全全部部材材料料时时间间导导数数，其其结结果果是是在在更更新新L L格格式式中中EulerianEulerian描描述述和和ALEALE描描述述之之间间唯一唯一区分是材料时间导数项。区分是材料时间导数项。与与在在第第4 4章章中中建建立立LagrangianLagrangian格格式式主主要要区区分分是是，现现在在需需要要以以偏偏微微分分方方程程(即即连连续续方方程程)形形式式考考虑虑质质量量守守恒恒方方程程，因因为为域域随随时时间改变，质量亦随时间改变间改变，质量亦随时间改变。所所以以，我我们们几几乎乎总总是是在在处处理理两两个个系系统统偏偏微微分分方方程程：标标量量连连续续方方程程和和向向量量动动量量方方程程。当当它它们们与与热热交交换换或或者者其其它它能能量量转转换换耦耦合时，还必须包含合时，还必须包含能量方程能量方程。第12页4 ALE控制方程控制方程连续方程连续方程(质量守恒质量守恒)或者或者动量方程动量方程能量方程能量方程自然边界条件自然边界条件在在 上上在在 上上基本位移边界基本位移边界在在 上上在在 上上初始条件初始条件第13页5 弱形式弱形式有限元近似有限元近似 对于单元对于单元e，ALE坐标给出为坐标给出为 单元单元e坐标坐标 网格运动给出为网格运动给出为 节点运动节点运动 这代表两个映射复合：从母单元到这代表两个映射复合：从母单元到ALE映射和网格运动映射映射和网格运动映射 网格速度为网格速度为 节点节点I 网格速度网格速度 第14页Lagrangian、Eulerian、ALE和自然坐标域之间映射和自然坐标域之间映射 5 弱形式弱形式有限元近似有限元近似 即单元坐标、网格坐标、空间坐标和材料坐标之间映射即单元坐标、网格坐标、空间坐标和材料坐标之间映射 第15页5 弱形式弱形式有限元近似有限元近似 在在ALE格式中，密度也是一个非独立变量。密度被近似为格式中，密度也是一个非独立变量。密度被近似为密度形状函数，可能不一样于网格运动形状函数密度形状函数，可能不一样于网格运动形状函数 速度材料时间导数速度材料时间导数 在离散运动方程中动力学项将以网格加速度在离散运动方程中动力学项将以网格加速度v,t表示，表示，经过积分和插值，给出经过积分和插值，给出材料速度材料速度为为应用类似插值，得到应用类似插值，得到传递速度传递速度为为 第16页5 弱形式弱形式有限元近似有限元近似 由传递速度公式由传递速度公式 得到得到 材料速度材料速度ALE时间导数时间导数 结论是，在弱形式中材料速度时间导数为结论是，在弱形式中材料速度时间导数为 对于密度材料时间导数，应用一样过程，给出对于密度材料时间导数，应用一样过程，给出第17页5 弱形式弱形式有限元矩阵有限元矩阵 连续方程连续方程 容量、转换和散度矩阵分别为容量、转换和散度矩阵分别为 动量方程动量方程 M和和L分分别别是广是广义质义质量和量和传递传递矩矩阵阵，对应对应于在参考构形下速度于在参考构形下速度 第18页5 弱形式弱形式有限元矩阵有限元矩阵 动量方程动量方程 注注意意到到除除了了它它们们是是以以变变分分形形状状函函数数形形式式定定义义之之外外，内内部部和和外外部部节节点点力力与与更更新新Lagrangian格格式式（框框4.3）中中对对应应项项是是一一致致。质质量量矩矩阵不是时间常量，因为密度和域随时间改变阵不是时间常量，因为密度和域随时间改变。第19页6 网格更新算法网格更新算法 在在ALEALE中中，网网格格能能够够任任意意移移动动给给出出了了大大变变形形可可能能性性。经经过过ALEALE移移动动边边界界（指指物物理理表表面面）能能够够利利用用LagrangianLagrangian准准确确特特征征来来循循迹迹，内内部部网网格格也也能能够够移移动动以以防防止止过过渡渡单单元元扭扭曲曲。然然而而这这需需要要一一个个有有效效算算法法来来更更新新网网格格，即即网网格格速速度度 必必须须给给定定，以以防防止止网网格格扭扭曲曲和和确保边界和接触面最少局部地保持确保边界和接触面最少局部地保持LagrangianLagrangian 。对对应应于于域域边边界界在在每每一一时时刻刻均均已已知知分分析析，预预先先给给定定网网格格运运动动。当域边界有一个已知运动时，网格随这一边界运动能够预先给定。当域边界有一个已知运动时，网格随这一边界运动能够预先给定。第20页6 网格更新算法网格更新算法建立材料和网格速度关系，只有给定其中一个，自动建立材料和网格速度关系，只有给定其中一个，自动确定确定另一个另一个1 假如给定假如给定，能够计算位移，能够计算位移 和加速度和加速度 。应用显式积分算法和中心差分算法，而不需要计算相对速度应用显式积分算法和中心差分算法，而不需要计算相对速度w。2 假如假如 未知，给定了未知，给定了 w，在更新网格前求解上式计算，在更新网格前求解上式计算 。发展发展ALE关键问题是给定这些速度最正确选择和更新网格算关键问题是给定这些速度最正确选择和更新网格算法。法。3 给定给定 和和 w 分量形式，混合算法分量形式，混合算法。第21页1.当当域域边边界界有有一一个个已已知知运运动动时时，网网格格随随这这一一边边界界运运动动能能够够预先给定。预先给定。6 网格更新算法网格更新算法2.Lagrange-Euler 矩阵方法，任意定义相对速度（参考）矩阵方法，任意定义相对速度（参考）Lagrange-Euler 参数矩阵参数矩阵假如假如 ，则，则 ，能够是时间和空间函数。能够是时间和空间函数。由上式，相对速度是材料速度线性函数，由上式，相对速度是材料速度线性函数，假如假如 ，则，则 ，Lagrangian网格描述。网格描述。假如假如 ，则，则 ，Eulerian网格描述。网格描述。第22页6 网格更新算法网格更新算法由传递速度和相对速度由传递速度和相对速度得到得到 在参考域边界上必须满足后一个方程。得到一个网格再分区在参考域边界上必须满足后一个方程。得到一个网格再分区基本方程：基本方程：在二维情况下显式格式以下：在二维情况下显式格式以下：问问：或或分别对应什么格式？分别对应什么格式？第23页6 网格更新算法网格更新算法 带有网格更新带有网格更新ALE技术基于技术基于L-E参数，对参数，对表面波动表面波动问题是非问题是非常有用。我们假设自由表面相对于总体坐标是有导向，曲面方常有用。我们假设自由表面相对于总体坐标是有导向，曲面方程能够写为程能够写为 欧拉坐标用于欧拉坐标用于 方向，方向，自由表面经过一个空间坐标定义，它对其余两个空间坐标和自由表面经过一个空间坐标定义，它对其余两个空间坐标和时间是连续可微函数，时间是连续可微函数，L-E矩阵只有一个非零项矩阵只有一个非零项 普通等于普通等于1称为累计率函数，表示在自由表面得到或失去质量。称为累计率函数，表示在自由表面得到或失去质量。自由表面是物质表面；沿着表面累计率必须为零；所以自由表面是物质表面；沿着表面累计率必须为零；所以 等于等于1第24页6 网格更新算法网格更新算法在贮箱内液体晃动在贮箱内液体晃动自由表面液体晃动自由表面液体晃动第25页Coupled equationsCoupled equations Free surface Free surface Fluid-structure interfaceFluid-structure interface S Sf fS Sww 1 On S1 On Sw w:the geometrical the geometrical compatibility conditionscompatibility conditions(slipping slipping boundary conditionboundary condition)are applied)are applied:2 On S2 On Sw w:the equilibrium conditionsthe equilibrium conditions are appliedare applied:6 网格更新算法网格更新算法第26页Elephant foot bulging(EFB)Diamond shape bulging(DSB)流体流体/结构耦合分析结构耦合分析6 网格更新算法网格更新算法第27页流体流体/结构耦合分析结构耦合分析程序平台：程序平台：ABAQUS/用户单元附加质量用户单元附加质量例：动力作用下液体贮箱例：动力作用下液体贮箱Elephant foot bulging Diamond shape bulgingEFB and DSBEFB and DSBEFB and DSBEFB and DSB6 网格更新算法网格更新算法第28页程序平台：程序平台：ABAQUS/ALE流体单元流体单元例例2：动力作用下液体贮箱：动力作用下液体贮箱Elephant foot bulgingElephant foot bulgingDeformation between test data and FEM by added mass and ALE流体流体/结构耦合分析结构耦合分析6 网格更新算法网格更新算法第29页6 网格更新算法网格更新算法 因为因为 方法对表面单元跟踪很好，不过极难确保流体域内部方法对表面单元跟踪很好，不过极难确保流体域内部单元扭曲。因为这个缺点，引入了一个混和方法，变形梯度法，一旦单元扭曲。因为这个缺点，引入了一个混和方法，变形梯度法，一旦边界已知，经过给定网格位移或者速度，防止单元缠结和扭曲。边界已知，经过给定网格位移或者速度，防止单元缠结和扭曲。3.3.变形梯度法混合算法变形梯度法混合算法 因为网格位移和速度直接控制单元形状，因为网格位移和速度直接控制单元形状，沿着域边界给定沿着域边界给定 ，在内部给定节点位移或者速度在内部给定节点位移或者速度，混合算法。，混合算法。ALE网格，边界处，网格，边界处，Lagrangian化；化；在内部，在内部，Eulerian化。化。第30页 作作为为一一个个使使用用修修正正弹弹性性方方程程ALEALE网网格格更更新新例例子子，考考虑虑一一个个位位于于一一个个矩矩形形流流体体域域中中圆圆柱柱绕绕流流有有限限元元网网格格。在在圆圆柱柱沿沿着着y y方方向向移移动动了了一一个个位位移移0.25w0.25w后后网网格格见见图图b b。矩矩形形域域边边界界保保持持固固定定。由由图图可可见见，圆圆柱柱附附近近网网格格精精度度在在更更新新后后网网格格中中得得以以保保持持，而而且且无显著单元扭曲发生。无显著单元扭曲发生。网格更新例子网格更新例子 6 网格更新算法网格更新算法第31页7 Petrov-Galerkin方法方法 伽辽金伽辽金(Galerkin)(Galerkin)方法，是利用满足位移和应力边界条件方法，是利用满足位移和应力边界条件函数寻求积分方程解答。函数寻求积分方程解答。关键问题是关键问题是增加粘性增加粘性，消除不稳定项，确保数值稳定性。，消除不稳定项，确保数值稳定性。建建立立Petrov-Galerkin方方法法迎迎风风流流线线(Streamline Upwind Petrov-GalerkinSUPG)公公式式。对对流流扩扩散散方方程程是是一一个个有有用用方方法法，它它对对应应于于动动量量方方程程线线性性化化。对对于于离离散散稳稳态态对对流流扩扩散散方方程程，将将得得到到闭闭合合解解答答。将将证证实实当当网网格格参参数数(已已知知Peclet数数)超超出出临临界界值值时时，这这个个解解答答在在空空间间是是振振荡荡。经经过过建建立立P-G方方法法以以消除这些振荡，即纠正不稳定性。消除这些振荡，即纠正不稳定性。在一维中，离散方程类似于标准迎风方程。然而在多维中，在一维中，离散方程类似于标准迎风方程。然而在多维中，它们提供了沿着流线引导迎风项一致理论框架。它们提供了沿着流线引导迎风项一致理论框架。第32页 迎迎格格调调式式基基本本思思想想：当当用用差差分分方方程程求求解解偏偏微微分分方方程程时时，利利用用特特征征线线方方向向一一侧侧单单边边差差商商来来代代替替空空间间偏偏导导数数。目目标标是是确确保保数数值稳定性。如一阶线性常系数双曲型方程：值稳定性。如一阶线性常系数双曲型方程：7 Petrov-Galerkin方法方法在特征线方向一侧单边差商来代替偏导数，迎格调式为：在特征线方向一侧单边差商来代替偏导数，迎格调式为：j-1n+1a0njj+1j-1n+1a0njj+1第33页 中中心心差差分分格格式式需需要要加加入入非非线线性性人人工工粘粘性性项项以以去去掉掉激激波波前前振振荡，同时加入荡，同时加入线性数值粘性项线性数值粘性项以确保计算稳定。以确保计算稳定。迎迎格格调调式式依依据据EulerEuler方方程程中中波波传传输输信信息息结结构构格格式式，无无需需添添加加人人工工粘粘性性项项，迎迎格格调调式式为为空空间间差差分分格格式式(偏偏心心差差分分格格式式)，理理论论上比中心差分格式扎实。上比中心差分格式扎实。7 Petrov-Galerkin方法方法 假假如如差差分分格格式式(所所用用网网格格点点)与与微微分分方方程程特特征征线线方方向向一一致致，那那么么网网格格比比 在在满满足足一一定定条条件件下下是是稳稳定定，不不然然，差差分分格格式式是不稳定。是不稳定。双曲线型微分方程特征线双曲线型微分方程特征线存在交叉，导数不连续。存在交叉，导数不连续。第34页7 Petrov-Galerkin方法方法稳态线性对流扩散方程为稳态线性对流扩散方程为运动粘度运动粘度 给定速度给定速度 对于一维问题，偏微分方程成为常微分方程对于一维问题，偏微分方程成为常微分方程 应用边界条件应用边界条件这是在这是在0 xL 域上两点边值问题，轻易证实公式准确解答是域上两点边值问题，轻易证实公式准确解答是空间非独立变量空间非独立变量 第35页7 Petrov-Galerkin方法方法应用线性形状函数建立应用线性形状函数建立Galerkin离散，并在全域上积分离散，并在全域上积分 变分函数变分函数 分部积分而且应用散度原理，对流扩散方程弱形式是分部积分而且应用散度原理，对流扩散方程弱形式是将域（将域（0，L）划分成相同尺寸单元）划分成相同尺寸单元 在每个单元上离散方程给出为在每个单元上离散方程给出为 有限元形状函数有限元形状函数 对于第对于第j个节点内部方程为个节点内部方程为 第36页7 Petrov-Galerkin方法方法上式上式恰好是中心差分方程，能够方便地重写成为恰好是中心差分方程，能够方便地重写成为1 假如假如Peclet数小于数小于1，则离散解答是类似于准确解答；则离散解答是类似于准确解答；2 假如假如Peclet数大于数大于1，则离散解答是正或者负而振荡。则离散解答是正或者负而振荡。这种不稳定是数值离散空间不稳定。这种不稳定是数值离散空间不稳定。Peclet数数 准确解准确解离散解离散解第37页将上式中变分函数将上式中变分函数重写为重写为式中式中选择选择 是为了消除振荡，期望得到准确解。是为了消除振荡，期望得到准确解。如如简化成为中心差分方法简化成为中心差分方法 是一个完全迎风公式是一个完全迎风公式 7 Petrov-Galerkin方法方法应用线性形状函数建立应用线性形状函数建立Galerkin离散，并在全域上积分离散，并在全域上积分 逆逆P-G项项G项项第38页总结总结 ALE网格，调整网格运动相对速度。网格，调整网格运动相对速度。边界处，边界处，Lagrangian化；化；在内部，在内部，Eulerian化；化；稳定性，增加粘性。稳定性，增加粘性。第39页