首页/文章/ 详情

形象理解弹性力学中的张量

2小时前浏览10

概述

  弹性力学中有很多张量,尽管这些张量都能展开放到直角坐标系中理解,但是要想在弹性力学的基础上深入学习力学,张量的深刻理解和熟练运用是必须的。

        帖子转载自知乎,内容深入浅出,从矢量的角度引入张量概念,不让人觉得突兀,值得一读。已征得贴主同意转发,点击文末“阅读原文”即可跳转原文。

正文

  我开始学习张量的时候认为,二阶张量    矩阵    二维数组。这实际上是不准确的,因为张量实际上是矢量在维度上的扩展,它不仅有大小。还有方向。我们通常说的矢量    中,    和    一般指的是分量,表示的是当前坐标系下矢量在各个方向上的分量,通常取的是笛卡尔坐标系,但实际上也可以换成其他坐标系。在二维的空间中极坐标系中的向量    和直角坐标系中的    实际上表示的是两个相同的矢量。  那么二阶张量其实也是一样,在不同的坐标系下它的表达形式会有所不同,矩阵只是它的分量表达形式,矩阵的规模取决于张量的维度,一维空间就是1×1,二维空间就是2×2,三维空间就是3×3,而矩阵的各个数值代表了它的各个分量的大小,分量结合基矢量才是这个二阶张量本身,也就是说,对于不同的坐标系,这个用于表示张量分量的矩阵的数值会发生变化。
  下面以应力为例,说明一下二阶应力张量表示的含义。对于一个一维问题,应力和中学时候学到的压强的概念是类似的,就是在这个方向上单位面积受到的力的大小。对于多维问题,在笛卡尔坐标系下,应力张量可以写成    ,其中    表示张量的分量,其中的 ij分别独立地表示不同的维度,在三维笛卡尔坐标系中分别取    。于是应力分量    可以写成矩阵的形式

 

  而    和    分别表示基矢量,在这里可以通俗地理解为当前坐标系下的各个方向。上式中主对角元上的三个分量可以理解为在分别在三个坐标轴对应的法平面上的单位面积沿着三个坐标轴方向上的力,如下图所示,它们会使得单元体在对应的方向上受到拉伸或者压缩,即正对该方向的变形,因此被称为正应力;而非对角元上的分量则表示在法平面上的单位面积沿着其他方向的力,它会使得单元体在这个方向上发生错位,因此被称为切应力。  由于这个单元体的体积假设为无限小,它可以看做是一个点,在这个点上受力的状态是确定的,但是将坐标系的方向旋转一下,或者取其它的坐标系,那么它的各个分量就会发生变化,例如每个应力状态都有三个主应力,对应三个主方向,在三个方向组成的坐标系上,它们的切应力都为0,只有主应力。同样地,在球坐标系或柱坐标系中,它们的分量也会有所不同,但是实际他们的应力状态是客观存在的,不会发生任何变化。应力张量描述的也是这个点本身的应力状态,和选取的坐标系没有关系。也就是说,在坐标转换的过程中,实际的应力的状态也就是说应力张量并没改变,改变的只是我们的度量方式,或者说基矢量的选取。因此这也是用张量计算的意义:张量不随坐标变换而发生改变,在任何坐标系下相应的张量表示的公式都是成立的。当然,系统地学习张量还需要讲解张量的运算、斜角坐标系以及各种曲线坐标系以及各种度量张量,这里只是为了便于理解,不做叙述了。
  此外,应力张量的各个分量就表示在第一个角标表示的平面内单位面积沿着第二个角标表示的方向上的力,它一共包含了两个“方向”,因此是二阶的张量。类似的,力矢量中的各个分量表示的是沿着该方向上的力的大小,只有一个“方向”,因此是一阶张量。这里的“方向”其实就是基矢量,也就是上文说的    或者在正交的笛卡尔坐标系中,基矢量分别为沿着三个坐标轴正方向上的单位矢量。
  和应力张量类似,应变张量也是个二阶张量,可以写成    。正应变,例如    表示的是单元体沿着x方向上单位长度在 xy方向的形变量,可以理解为单元体在某个方向上拉伸或者压缩变形相对于初始尺寸的比例;剪应变,例如表示的是单元体沿着    方向上单位长度在    方向的形变量(这里的剪应变和材料力学中的剪应变    的关系为    ,可以理解为单元体在某个方向上沿着另一个发生错位的变形程度,由于它表示的是变形的尺寸和物体本身尺寸的比值,因此没有单位。和应力张量一样,应变张量的描述也涉及到了两个“方向”,因此也是二阶张量。
  对于弹性张量,或者非弹性材料中的本构张量,它是一个四阶张量可以写成    ,乍一看很抽象,完全脑补不出来。实际上它可以理解为:材料发生变形时会产生应力,应力的大小可以由变形的程度确定,例如一维的问题中应变和应力的关系就可以类比弹簧中的    ,有  ,    为杨氏模量。弹性张量就是建立应力张量和应变张量之间的联系的,其中弹性张量的前两个指标由应力确定,后两个指标由应变确定,可以看做为:应变    对应力    的影响程度,那么在笛卡尔坐标系下的    的计算公式为

 

  其他的应力分量也可以以此类推,最后的形式为(这里为了便于书写把3×3的应力和应变张量写成了1×9的形式):

 

  又因为剪应力互等定理,    ,    ,    。所以为了便于程序编写和表示方便,将对称部分合并,于是应变和应力分量变成了    个,弹性矩阵变成了    的形式,也就是    标记

 

  这也刚好解释了为什么要定义前文提到的    。再将分量

 

分别当做第一到第六分量,就可以得到:

 

  当然需要说明的是,弹性张量的严格推导应该从应变能密度出发,是应变能密度对应变求二阶导数得到的,即

  因为求导顺序不影响求导结果,因此


 

  即上式的弹性矩阵是一个对称矩阵。

来源:有限元先生

材料
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-11-03
最近编辑:2小时前
外太空土豆儿
硕士 我们穷极一生,究竟在追寻什么?
获赞 2粉丝 1文章 45课程 0
点赞
收藏
作者推荐

newmark法求解运动方程(附matlab代码)

问题引入在数值仿真中,常常涉及到静力和动力的问题。静力问题通常忽略物体的惯性力和阻尼效应。也就是说,系统中的力只与物体的静态形变和外部载荷有关,而不考虑动态过程中速度和加速度带来的影响,因此静力问题可以用方程AX=B表示,这种方程用线性代数里面的内容就可以求解,主要是A矩阵求逆的计算速度决定了方程的求解速度,这部分的资料很多,在此不多说,这个帖子重点讲解动力问题。动力问题极为广泛,动力问题(DynamicProblem)是指在力学和物理学中研究物体或结构在外力作用下随时间变化的运动、变形和应力响应。与静力问题不同,动力问题考虑了系统的惯性力和阻尼效应,这使得其分析需要解决随时间变化的运动方程。动力问题广泛应用于涉及冲击、振动、地震、波动等现象的工程和物理系统中。动力问题的特点考虑惯性力和加速度:在动力问题中,物体或结构的运动和变形不仅受外力和内力的影响,还受到惯性力的作用。根据牛顿第二定律,惯性力与物体的加速度成正比,因此在动力分析中必须考虑加速度的影响。时间依赖性:动力问题的解通常是时间的函数,要求在整个时间域内求解物体或结构的运动状态。随着时间的推移,外力、位移、速度和应力等物理量会发生变化。系统响应的类型:自由振动:不受外部力作用的振动,仅依赖于系统的初始状态和自然特性。强迫振动:在外力作用下,系统的响应随外力变化而变化,通常伴随频率和振幅的变化。冲击和瞬态响应:系统受到瞬时力(如冲击或爆炸)时的快速响应,通常表现为瞬态振动。阻尼效应:动力系统中还可能存在阻尼,阻尼力通常与速度成正比,并会导致系统的运动逐渐衰减。动力问题中必须考虑阻尼效应对系统振动和响应的影响。动力问题的基本方程:动力问题的基本方程是牛顿第二定律或达朗贝尔原理的推广,称为运动方程(EquationofMotion)。一个简单的动力问题可以用以下形式的运动方程来描述:这个方程说明了外力F(t)如何引起系统的位移、速度和加速度的时间响应。这是一个二阶微分方程,求解这个方程的方法大致分为解析法和数值解法。解析解当然是最精确的,但只有简单的动力系统或边界条件较为规则的问题,动力问题可以通过解析方法求解。例如:单自由度系统的自由振动:对于简单的弹簧-质量系统,运动方程可以通过特征值问题解析求解,得到系统的固有频率和振动模态。正弦激励下的强迫振动:当系统受到周期性外力(如正弦波)作用时,可以通过解析方法求解系统的稳态响应。数值解法是对于复杂的几何形状、边界条件或载荷,解析解往往难以得到,因此通常采用数值方法来求解动力问题。常见的数值方法包括:有限元法(FEM):将连续体离散为有限个单元,通过求解离散的动力方程来求解整个结构的动态响应。Abaqus等软件广泛使用有限元方法进行复杂动力学分析。有限差分法(FDM):将动力方程的时间和空间导数进行离散化,用数值方法逐步推进时间步长,求解各个时间点的响应。显式和隐式时间积分法:显式方法:如中心差分法,适合处理瞬态问题和高频振动问题。隐式方法:如Newmark方法,适合求解稳定性要求较高的动力问题。这里就重点介绍newmark法求解运动方程,包括固定位移边界条件的施加。newmark理论及程序newmark法是一种常用的时间积分方法,用于求解结构动力学中的运动方程,尤其是在有限元分析中求解动力问题(如振动、冲击、地震响应等)。该方法由NathanM.Newmark于1959年提出,它通过离散时间步长,将动力方程转化为一系列可以逐步求解的代数方程,从而计算出位移、速度和加速度的时间响应。newmark法的优点在于它既适用于求解线性问题,也适用于求解非线性问题,并且它在保持数值稳定性和精度之间提供了一种灵活的平衡。根据参数的选择,Newmark法可以实现不同的时间积分方案,如条件稳定和无条件稳定。newmark公式资料非常的多,基本随便拎起来一本涉及到动力学的书籍都会给出该方法的详细公式,如下图所示上面的资料来自于朱伯芳院士的《有限单元法原理与应用》第四版的第499页。朱伯芳院士在北京于2024年2月4日逝世,每次我想起朱院士序言中的内容,我都会对朱院士充满敬意。朱院士为我国大江大河治理和水电开发作出了卓著贡献,我在此借这个帖子沉痛悼念朱伯芳院士!深切缅怀朱伯芳院士!朱院士的书是我学习有限元过程中推公式用的最多的一本书,当时我是看了一些国内其他书籍,觉得完蛋了,学不会有限元了,感觉山重水复疑无路了,但是看到朱院士的这本书,顿时豁然开朗,柳暗花明。这里贴一下书的封面,有需要的可以私信我。书归正传,有了newmark公式之后,先不慌着写程序,我们通过自己的理论求解了刚度矩阵、质量矩阵和阻尼矩阵之后,还要添加边界条件,才能消除刚体位移,至于添加边界条件又是一个很大的领域,我涉及不多,这里重点讲解newmark,就以最简单的固定位移边界条件为例。添加固定边界的方法有乘大数法、划零置一法和划行划列法等等,我在这里采用的是行划列法,即假如模型的第7个节点是全固定的,就把刚度矩阵等的第3*7-2:3*7(三维问题中)的行和列都直接删除,然后将矩阵传进newmark法中进行求解,求解完成后,再将自由度还原,方便后处理,这里贴一下主要的matlab程序。function[U,V,A]=newmark(deltat,beta,kk,mm,cc,ndim,fbc,kinc,ndofel,nodeu)%结果初始化uu=zeros(ndofel-ndim*length(nodeu),kinc);%didpvv=zeros(ndofel-ndim*length(nodeu),kinc);%velaa=zeros(ndofel-ndim*length(nodeu),kinc);%acc%%基本参数计算alpha=0.25*(0.5+beta)*(0.5+beta);%newamrk计算参数a0=1/(alpha*deltat*deltat);a1=beta/(alpha*deltat);a2=1/(alpha*deltat);a3=1/2/alpha-1;a4=beta/alpha-1;a5=deltat/2*(beta/alpha-2);a6=deltat*(1-beta);a7=beta*deltat;%%迭代计算ek=a0*mm+a1*cc+kk;%有效刚度矩阵ekk=inv(ek);%有效刚度矩阵求逆,只进行一次fori=2:kincfprintf("kinc:%dtime:%f\n",i,deltat*(i-1));ft=fbc(:,i);u=uu(:,i-1);v=vv(:,i-1);a=aa(:,i-1);%计算i+1时刻的有效荷载ft=fbc(:,i);f=ft+mm*(a0.*u+a2.*v+a3.*a)+cc*(a1.*u+a4.*v+a5.*a);%计算i+1时刻的位移uu(:,i)=ekk*f;%计算i+1时刻的速度和加速度aa(:,i)=a0.*(uu(:,i)-u)-a2.*v-a3.*a;vv(:,i)=v+a6.*a+a7.*(aa(:,i));end%%还原自由度U=zeros(ndofel,kinc);V=U;A=U;fori=1:length(nodeu)node=nodeu(i);U(ndim*node-(ndim-1):ndim*node,:)=4e20;%在位移向量中标记固定边界的位置V(ndim*node-(ndim-1):ndim*node,:)=4e20;%在速度向量中标记固定边界的位置A(ndim*node-(ndim-1):ndim*node,:)=4e20;%在加速度向量中标记固定边界的位置end%还原自由度index=0;fori=1:ndofel/ndimifU(ndim*i-(ndim-1),1)==4e20index=index+1;elseU(ndim*i-(ndim-1):ndim*i,:)=uu(ndim*(i-index)-(ndim-1):ndim*(i-index),:);V(ndim*i-(ndim-1):ndim*i,:)=vv(ndim*(i-index)-(ndim-1):ndim*(i-index),:);A(ndim*i-(ndim-1):ndim*i,:)=aa(ndim*(i-index)-(ndim-1):ndim*(i-index),:);endendend注意传进newmark函数的刚度矩阵等都是已经添加过边界条件的。数值算例设计悬臂梁受动荷载算例,悬臂梁一端固定,另一端受动荷载,统计悬臂端下边界的重点数据作对比,示意图为荷载曲线为计算总时长为10s,增量步长为0.01s,采用六面体网格离散,离散后的网格图为采用C3D8单元,共计5120个六面体单元,节点总数为6273。设置两种工况工况一:采用abaqus计算全部的结果工况二:采用自编newmark程序计算统计了悬臂端下边界中点的位移曲线如下图速度曲线如下图加速度曲线如下图自编newmark程序求解运动方程,与abaqus的计算结果相对比,最终结果一致(也有可能是碰巧对准了),但是依然不影响在这过程中对newmark加深了理解。详细的数据和参数不在此罗列了,有兴趣的可以私信我获取全部数据。点击卡片关注我们来源:有限元先生

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈