平面应力问题7

平面应力问题

平面应力问题的平面应力应变关系

{σxxσyyτxy}=E1γ2[1γ0γ1000(1γ)2]{εxxεyyγxz}⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪σxxσyyτxy⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪=E1γ2⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢1γ0γ1000(1γ)2⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪εxxεyyγxz⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪

平面应变问题的应力-应变关系

{σxxσyyτxy}=E(1+v)(12v)[1vv0v1v000(12v)2]{ϵxxϵyyγxy}⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪σxxσyyτxy⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪=E(1+v)(12v)⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢1vv0v1v000(12v)2⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ϵxxϵyyγxy⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪

小变形下的应力-位移关系

εxx=uxεyy=νyγxy=uy+νxεxx=uxεyy=νyγxy=uy+νx

写成矩阵形式

{εxxεyyγxy}=[x00yyx]{uν}⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪εxxεyyγxy⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢x00yyx⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥uν

简写

{ϵ}=[L]U{ϵ}=[L]U

对于n节点的单元,未知位移的节点近似插值函数为

u=N1u1+N2u2++Nnunν=N1ν1+N2ν2++Nnνn

矩阵形式

{uv}=[N10|N20||Nn00N1|0N2||0Nn]{u1ν1ν2u2ν2unνn}

简写

{U}=[N]{a}

用上式中U代替 {ϵ}=[L]U中的U

{ϵ}=[L]U=[L][N]{a}=[B]{a}

其中

[B]=[N1x0|N2x0||Nnx00N1y|0N2y||0NnyN1yN1x|N2yN2x||NnyNnx]

矩阵[B] 被称为应变矩阵,可以通过插值函数Ni(x,y)求导得到

为了得到单元单元上的载荷与节点位移之间的关系,我们用虚功原理

Veδ{ϵ}T{σ}dV=Veδ{U}T{b}dV+Γeδ{U}T{t}dΓ+iδ{U}T({x}={¯x}){P}i

式中:

{ϵ}为应变矢量;{σ}为应力矢量;{U}为位移矢量

{b}为体力矢量;{t}为表面力矢量

{P}i为作用在 {x}={¯x} 处的集中力矢量

dV为单元体积;dΓ 为表面力{t}所作用单元的单元边界

应变的变分 δ{ε}和位移的变分 δ{U} 可以分别表示为

δ{ϵ}=δ([B]{a})=[B]δ{a}δ{U}=δ([N]{a})=[N]δ{a}

应力应变关系

{σ}=[D]{ϵ}=[D][B]{a}

代入虚功原理方程

Veδ{a}T[B]T[D][B]{a}dV=Veδ{a}T[N]T{b}dV+Γeδ{a}T[N]T{t}dΓ+iδ{a}T[N((x)=(ˉx))]T{P}iVeδ{ϵ}T{σ}dV=Veδ{U}T{b}dV+Γeδ{U}T{t}dΓ+iδ{U}T({x}={¯x}){P}i

注意,对于平面单元,单元体积 dV和单元边界 dГ 可以分别写成 dν=tdA和 dΓ=tdl ,式中t
为单元的厚度,dA为单元的无限小面积,dl为单元的无限小边界长度。
由于δ{a} 是节点位移的变分,因此关于坐标独立,可以把它从积分符号中拿出来并消除,方程就变为

[Ac[B]T[D][B]tdA]{a}=Ac[N]T{b}tdA+Lc[N]T{t}tdl+i[N({x}={¯x})]T{P}i

写成矩阵形式

[Ke]{a}=fe

[Ke]=[Ae[B]T[D][B]tdA]

上式为刚度矩阵,其中

{fe}=Ae[N]T{b}tdA+Le[N]T{t}tdl+i[N((x)={¯x})]T{P}i

上式为单位力矢量

空间离散

常应变三角形(CST)单元

单元的插值函数

N1(x,y)=m11+m12x+m13yN2(x,y)=m21+m22x+m23yN3(x,y)=m31+m32x+m33y

其中

m11=x2y3x3y22Am12=y2y32Am13=x3x22Am21=x3y1x1y32Am22=y3y12Am23=x1x32Am31=x1y2x2y12Am32=y1y22Am33=x2x12A

A=12det[1x1y11x2y21x3y3]

image-20240714191212161

位移场

三角形单元的位移场可以表示为

u=N1u1+N2u2+N3u3ν=N1ν1+N2ν2+N3ν3

矩阵表示

{uν}=[N10|N20|N300N1|0N2|0N3]{u1ν1u2ν2u3ν3}

简写为

{U}=[N]{a}

应变矩阵

{ϵ}=[B]{a}

其中

=[N1x0|N2x0|N3x00N1y|0N2y|0N3yN1yN1x|N2yN2x|N3yN3x]

代入得

[B]=[m120|m220|m3200m13|0m23|0m33m13m12|m23m22|m33m32]

备注:矩阵[B]与笛卡儿坐标系的x和 y 坐标无关,它只是节点坐标的函数,并且在整个单元内为常量,因此在整个单元内应变矢量也是常量。这就是为什么这种单元被称为“常应变三角形单元”的原因,本书中将常应变三角形单元简称为 CST(Constant Strain Triangle)单元。

由于矩阵[B]和[D]都是常数矩阵,因此刚度矩阵可表示为

[Ke]=[B]T[D][B]tAe

其中 Ae为单元的面积

单位力矢量

体力

考虑体力{b}是重力引起的,因此

Aϵ[N]T{b}tdA=tAϵ[N100N1N200N2N300N3]{00ρg}dA=t[0AϵN1ρgdA0AϵN2ρgdA0AϵN3ρgdA]

这里采用式(8.18)和式(8.19)给出的三角形积分公式,计算在整个单元上的插值函数的积分。运用这些公式,可将积分结果整理如下:

AeN1ρgdA=ρgAeN11N02N03dA=ρg1!0!0!(1+0+0+2)!2Ae=ρgAe3AeN2ρgdA=ρgAeN01N12N03dA=ρg0!1!0!(0+1+0+2)!2Ae=ρgAe3AeN3ρgdA=ρgAeN01N02N13dA=ρg0!0!1!(0+0+1+2)!2Ae=ρgAe3

代入得

Ae[N]T{b}tdA=t3{0ρgAe0ρgAe0ρgAe}

可以发现,单元自重在各节点是平均分配的。

表面力

如图 所示的单元作用有大小为q的均布载荷,该均布载荷作用在单元的边 2-3 上,并
且与 x 轴的夹角为θ 。因此,表面力矢量可以表示为{t}={qcosθ,qsinθ}T

image-20240714192828717

则表面力可以表示为

Ln[N]T{t}tdl=L23[0000N200N2N300N3]{qcosθqsinθqsinθ}tdl=tL23{00N2qcosθN2qsinθN3qcosθN3qsinθ}dl

注意,在边 2-3 上 N1=0

采用式(8.18)给出的沿三角形边长的积分公式来计算沿长度方向的积分。运用这个公式上式可以表示为

Le[N]T{t}tdl=t{00qcosθL23/2qsinθL23/2qcosθL23/2qsinθL23/2}

从上式可以看出,节点2和节点3平均分配了作用在它们之间的均布载荷qL23

集中力

% 文件:CST_COARSE_MESH_DATA.m
%
% 下列变量被声明为全局变量,以便被构成程序的所有函数(M文件)使用
%
global nnd nel nne nodof eldof n     % 全局变量声明
global geom connec dee nf Nodal_loads  % 包括节点数、元素数、节点每元素数、每节点自由度数等

format short e  % 设置MATLAB的输出格式为科学计数法

% 定义模型参数
nnd = 21;              % 节点数
nel = 24;              % 元素数
nne = 3;               % 每个元素的节点数
nodof = 2;             % 每个节点的自由度数
eldof = nne * nodof;   % 每个元素的自由度数

% 节点坐标x和y
geom = zeros(nnd, 3);  % 初始化几何矩阵
% ... 此处省略了具体的节点坐标赋值 ...

% 元素连接信息
connec = zeros(nel, 3);  % 初始化连接矩阵
# ... 此处省略了具体的元素连接信息 ...

% 材料属性
E = 200000.;  % 弹性模量,单位MPa
vu = 0.3;      % 泊松比
thick = 5.;    % 梁厚度,单位mm

% 形成平面应力的弹性矩阵
dee = formdsig(E, vu);  % 调用formdsig函数,该函数未在代码中给出

% 边界条件
nf = ones(nnd, nodof);  % 初始化自由度矩阵
% ... 此处省略了具体的边界条件赋值 ...

% 计算自由度的总数
n = 0;
for i = 1:nnd
    for j = 1:nodof
        if nf(i, j) ~= 0
            n = n + 1;
            nf(i, j) = n;
        end
    end
end

% 载荷
Nodal_loads = zeros(nnd, 2);  % 初始化节点载荷矩阵
% ... 此处省略了具体的载荷赋值 ...

%%%%%%%%%%%%   输入结束   %%%%%%%%%%%%%
posted @   redufa  阅读(56)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· [翻译] 为什么 Tracebit 用 C# 开发
· 腾讯ima接入deepseek-r1,借用别人脑子用用成真了~
· Deepseek官网太卡,教你白嫖阿里云的Deepseek-R1满血版
· DeepSeek崛起:程序员“饭碗”被抢,还是职业进化新起点?
· 深度对比:PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
点击右上角即可分享
微信分享提示