计算力学|采用python进行有限元模拟

从abaqus输出的inp文件中读取节点和单元信息

import meshio

mesh = meshio.read('Job-3.inp')

coords = mesh.points###coords即为各个节点的坐标

Edof = mesh.cells_dict['triangle']#Edof为三角形单元的节点号

1.单元刚度矩阵

def element_stiffness(n1,coords,E,v,t):

    node1 = coords[n1[0]]

    node2 = coords[n1[1]]

    node3 = coords[n1[2]]

    xi = node1[0]

    yi = node1[1]

    xj = node2[0]

    yj = node2[1]

    xm = node3[0]

    ym = node3[1]

    A=abs((xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))*0.5)

    bi=yj-ym

    bj=ym-yi

    bm=yi-yj

    ci=-xi+xm

    cj=-xm+xi

    cm=-xi+xj

B=np.array([[bi,0,bj,0,bm,0],[0,ci,0,cj,0,cm],[ci,bi,cj,bj,cm,bm]])/(2*A)

    BT=B.transpose()

    D=np.array([[1,v,0],[v,1,0],[0,0,(1-v)/2]])*E/(1-v*v)

    Ke=BT.dot(D).dot(B)*t*A

    return Ke

2.总体刚度矩阵

def assemble_stiffness(coords,edof,E,v,t):

    NoN = np.size(coords, 0)

    NoE = np.size(edof, 0)

    PD = np.size(coords, 1)

    NPE = np.size(edof, 1)

    K=np.zeros([NoN*PD,NoN*PD])

    for i in range(len(edof)):

        n1 = edof[i, 0:(NPE)]

        Ke=element_stiffness(i,n1,coords, E,v,t)

        for r in range(NPE):

            for p in range(NPE):

                K[2*n1[r] , 2*n1[p] ] = K[2*n1[r] - 0, 2*n1[p]-0] + Ke[2*r, 2*p]

                K[2*n1[r] , 2*n1[p]+1] = K[2*n1[r] , 2*n1[p]+1] + Ke[2*r, 2*p+1]

                K[2*n1[r] +1, 2*n1[p] ] = K[2*n1[r]+1, 2*n1[p] ] + Ke[2*r+1, 2*p]

                K[2*n1[r] +1, 2*n1[p]+1] = K[2*n1[r]+1 , 2*n1[p]+1] + Ke[2*r+1, 2*p+1]

    return K

3.施加载荷

在薄板的左右两端施加均布载荷,首先根据单元节点横坐标判断在左右两端,判断在左右两端的单元个数,然后在单元上设置近似的均布载荷。

def F(coords,NoN,PD,b,c):

    G = np.zeros([NoN , PD])

    U = np.zeros([NoN , PD])

    # 施加左侧载荷

    m=0

    for i in range(NoN):

        if int(coords[i,:][0]*10**2) == -(b/2)*10**2+1.:

            m=m+1

    p=c/m*689.5

    k = 0

    for i in range(NoN):

        if int(coords[i][0]*10**2) == -(b/2)*10**2+1.:

            print(i)

            k = k + 1           

            G[i][0]= (-p)

    # 施加右侧载荷

    n=0

    for i in range(NoN):

        if int(coords[i,0:][0]*10**2) == (b/2)*10**2-1.:

            n = n + 1

    p = c/ n*689.5

    for i in range(NoN):

        if int(coords[i,0:][0]*10**2) == (b/2)*10**2-1.:

            G[i][0] = p

    return G

4.求解位移

U = np.zeros([2*NoN , PD-1])

T= np.linalg.inv(K)#T为K(总体刚度矩阵)逆矩阵

U=np.dot(T,GT)#U=[K]-1[G]

U=U.reshape(2*NoN,PD-1)#把位移转化为一维数组

5.计算应力应变

epsilonx = np.zeros(NoE)

epsilony = np.zeros(NoE)

epsilonxy = np.zeros(NoE)

A = np.zeros([3,3])

for i in range(NoE):

      dot_id = Edof[i]

      dot1 = np.array(coords[int(Edof[i, 0])])

      dot2 = np.array(coords[int(Edof[i, 1])])

      dot3 = np.array(coords[int(Edof[i, 2])])

      beta1 = dot2[1] - dot3[1]

      beta2 = dot3[1] - dot1[1]

      beta3 = dot1[1] - dot2[1]

      gamma1 = -dot2[0] + dot3[0]

      gamma2 = -dot3[0] + dot1[0]

      gamma3 = -dot1[0] + dot2[0]

      B=np.array([[1, dot1[0], dot1[1]], [1, dot2[0], dot2[1]],[1, dot3[0], dot3[1]]]) / 2

      A = np.linalg.det(B)

      epsilonx[i - 1] = (beta1 * U[dot_id[0]] + beta2 * U[dot_id[1]] + beta3 * U[dot_id[2]]) / (2 * A)

      epsilony[i - 1] = (gamma1 * U[NoN + dot_id[0]] +

            gamma2 * U[NoN + dot_id[1]] + gamma3 * U[NoN + dot_id[2]]) / (2 * A)

      epsilonxy[i - 1] = (beta1 * U[NoN + dot_id[0]] +

             beta2 * U[NoN + dot_id[1]] + beta3 * U[NoN + dot_id[2]]

             + gamma1 * U[dot_id[0]] + gamma2 * U[dot_id[1]] + gamma3 * U[dot_id[2]]) / (2 *A)

      sigmax = (E / (1 - v ** 2)) * (epsilonx + v * epsilony)

      sigmay = (E / (1 - v ** 2)) * (epsilony + v * epsilonx)

      sigmaxy = (E / (2 * (1 + v))) * epsilonxy

五、结果及分析

编程结果:

X方向位移(Abaqus)

y方向位移(Abaqus)

应力:

数值比较:最大x方向位移

编程

Abaqus

1.79345518e-03

1.183e-03

数值比较:最大y方向位移

编程

Abaqus

8.84718366e-04

1.838e-04

数值比较:最大x方向应力

编程

Abaqus

2.82313695e+04

1.028e+04

分析:在用Abaqus求解时,将中间整个圆环边界位移设置为零,相当于减少了薄板的位移量,因此出现了编程位移数值大于Abaqus数值的问题。而在使用相同模型时,位移的大小决定了应变的大小,应力呈现与应变相同的变化趋势,所以编程的应力大于Abaqus的应力,但结果在同一个数量级,可以基本确定是以上原因。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/894140.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

目标检测——Cascade R-CNN算法解读

论文: Cascade R-CNN: Delving into High Quality Object Detection (2017.12.3) 链接:https://arxiv.org/abs/1712.00726 Cascade R-CNN: High Quality Object Detection and Instance Segmentation (2019.6.24) 链接:https://arxiv.org/abs…

ubuntu22.04下GStreamer源码编译单步调试

前言 本文会通过介绍在linux平台下的GStreamer的源码编译和单步调试example实例。官网介绍直接通过命令行来安装gstreamer可以参考链接:Installing on Linux。 这种方法安装后,基于gstreamer的程序,单步调试的时候并不会进入到gstreamer源码…

LSTM预测:糖尿病的发生情况

本文为为🔗365天深度学习训练营内部文章 原作者:K同学啊 本期,做个二维结构化数据的分类预测。提到结构化数据,一般的分类算法常用有:逻辑回归(二分类)、KNN、SVM、决策树、贝叶斯、随机森林、X…

Jenkins配置流水线任务-实践操作(Pipeline-script)

Jenkins配置流水线任务-实践操作(Pipeline-script) 1、新增jenkins 任务,选择流水线 2、参数化 3、流水线配置 pipeline {agent anystages {stage(aoePlugin_mysql) {steps {echo "xxx,数据库:Mysql"echo "${HOST},${USER_NAME}"b…

王爽汇编语言第三版实验1

前言 本系列的文章是对王爽老师的汇编语言中的实验的解答记录,原书一共有17个实验,由于学校的教学流程只做到了第14个实验,因此本文章只会有前十四个实验的解答记录,还有个比较重要的是,文章中会有原书实验中没有的题目&#xff…

C语言 | Leetcode C语言题解之第477题汉明距离总和

题目&#xff1a; 题解&#xff1a; int totalHammingDistance(int* nums, int numsSize) {int ans 0;for (int i 0; i < 30; i) {int c 0;for (int j 0; j < numsSize; j) {c (nums[j] >> i) & 1;}ans c * (numsSize - c);}return ans; }

element plus的el-select分页

摘要&#xff1a; el-select的数据比较多的时候&#xff0c;必须要分页&#xff0c;处理方案有全部数据回来&#xff0c;或者添加搜索功能&#xff0c;但是就有个问题就是编辑的时候回显问题&#xff0c;必须要保证select的数据有对应的id与name匹配回显&#xff01; <el-fo…

如何用pyhton修改1000+图片的名字?

import os oldpath input("请输入文件路径&#xff08;在windows中复制那个图片文件夹的路径就可以):") #注意window系统中的路径用这个‘\分割&#xff0c;但是编程语言中一般都是正斜杠也就是’/‘ #这里写一个代码&#xff0c;将 \ > / path "" fo…

数字图像处理:图像复原应用

数字图像处理&#xff1a;图像复原应用 1.1 什么是图像复原&#xff1f; 图像复原是图像处理中的一个重要领域&#xff0c;旨在从退化&#xff08;例如噪声、模糊等&#xff09;图像中恢复出尽可能接近原始图像的结果。图像复原与图像增强不同&#xff0c;复原更多地依赖于图…

服务器数据恢复—服务器硬盘指示灯亮黄灯,raid崩溃的数据恢复案例

服务器数据恢复环境&#xff1a; 一台浪潮服务器中有一组由6块SAS硬盘组建的RAID。服务器上划分了1个卷&#xff0c;存放Oracle数据库文件。 服务器故障&检测&#xff1a; 服务器上有两个硬盘指示灯亮黄灯&#xff0c;RAID崩溃&#xff0c;服务器不可用。 将故障服务器中所…

LLM:deepspeed zero-2时模型训练所占显存分析

前置&#xff1a; fp16占2字节&#xff0c;fp32占4字节。换算就是1B的参数量&#xff0c;以fp16表示&#xff0c;占2G的内存。 模型参数为32B 全量微调&#xff1a; 模型参数&#xff1a;fp16的模型前向传播副本。fp32的模型的优化参数副本。这就是322324192G 梯度&#xff…

Jmeter简介

基础介绍 Jmeter录制脚本的原始是配置一个HTTP代理&#xff0c;然后浏览器通过这个代理访问测试页面从而完成脚本录制。 一、下载安装 jmeter本身不需要安装&#xff0c;需要配置环境变量JDK&#xff0c;然后打开bin文件夹中的jmeter.vbs即可。建议jdk 1.7及以上版本。 基本祖…

CVE-2024-22120:Zabbix低权限SQL注入至RCE+权限绕过

所有利用代码&#xff1a; GitHub - W01fh4cker/CVE-2024-22120-RCE: Time Based SQL Injection in Zabbix Server Audit Log --> RCE 一、漏洞环境搭建 1.1 下载vmware镜像并设置 直接懒人一键搭建&#xff1a; https://cdn.zabbix.com/zabbix/appliances/stable/6.0/6.0…

得物App3D创新应用引关注,世界设计之都大会启幕

近日&#xff0c;2024世界设计之都大会&#xff08;WDCC&#xff09;在上海盛大启幕。此次大会以“设计无界 新质生长”为主题&#xff0c;汇聚了全球设计领域的精英与前沿成果&#xff0c;展现了设计作为新质生产力的巨大潜力。主场展览占据了整整3个楼面&#xff0c;总面积达…

k8s-对命名空间资源配额

对k8s命名空间限制的方法有很多种&#xff0c;今天来演示一下很常用的一种 用的k8s对象就是ResourceQuota 一&#xff1a;创建命名空间 kubectl create ns test #namespace命名空间可以简写成ns 二&#xff1a; 对命名空间进行限制 创建resourcequota vim resourcequ…

基于Javaweb的医院挂号预约管理系统

系统展示 用户前台界面 管理员后台界面 医生后台界面 系统背景 在现代社会&#xff0c;随着医疗需求的不断增长&#xff0c;病患挂号成为医院面临的一大挑战。传统的挂号方式不仅耗时耗力&#xff0c;还容易引发混乱和不满。病患需要排队等候&#xff0c;挂号过程繁琐&#xff…

Nginx(Linux):启动停止Nginx

目录 1、理解Nginx后台进程2、停止Nginx(方式一&#xff1a;使用信号源)2.1 获取master进程号2.1 设置信号源 3、停止Nginx(方式二&#xff1a;使用命令行) 1、理解Nginx后台进程 Nginx后台进程包含master和worker两类进程。 master进程&#xff1a;主要用来管理worker进程&am…

Docker 教程四 (Docker 镜像加速)

Docker 镜像加速 国内从 DockerHub 拉取镜像有时会遇到困难&#xff0c;此时可以配置镜像加速器。 目前国内 Docker 镜像源出现了一些问题&#xff0c;基本不能用了&#xff0c;后期能用我再更新下。* Docker 官方和国内很多云服务商都提供了国内加速器服务&#xff0c;例如…

C++ | Leetcode C++题解之第479题最大回文数乘积

题目&#xff1a; 题解&#xff1a; class Solution { public:int largestPalindrome(int n) {if (n 1) {return 9;}int upper pow(10, n) - 1;for (int left upper;; --left) { // 枚举回文数的左半部分long p left;for (int x left; x > 0; x / 10) {p p * 10 x %…

Maxwell 底层原理 详解

Maxwell 是一个 MySQL 数据库的增量数据捕获&#xff08;CDC, Change Data Capture&#xff09;工具&#xff0c;它通过读取 MySQL 的 binlog&#xff08;Binary Log&#xff09;来捕获数据变化&#xff0c;并将这些变化实时地发送到如 Kafka、Kinesis、RabbitMQ 或其他输出端。…