📄 elmt_shell_4n.c
字号:
/*
* =============================================================================
* ALADDIN Version 1.0 :
* elmt_shell_4n.c : STATIC/DYNAMIC SHELL_4N Element
*
* Copyright (C) 1995 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin
* Institute for Systems Research,
* University of Maryland, College Park, MD 20742
*
* This software is provided "as is" without express or implied warranty.
* Permission is granted to use this software for any on any computer system
* and to redistribute it freely, subject to the following restrictions:
*
* 1. The authors are not responsible for the consequences of use of
* this software, even if they arise from defects in the software.
* 2. The origin of this software must not be misrepresented, either
* by explicit claim or by omission.
* 3. Altered versions must be plainly marked as such, and must not
* be misrepresented as being the original software.
* 4. This notice is to remain intact.
*
* -------------------------------------------------------------------
*
* The Algorithms are rederived in order to suit for Implicit Algorithms.
* For implicit integration: linear strain and displacement relation is
* employed
*
* -------------------------------------------------------------------
*
* Written by: Xiaoguang Chen January 1995
* =============================================================================
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "defs.h"
#include "units.h"
#include "matrix.h"
#include "vector.h"
#include "fe_database.h"
#include "symbol.h"
#include "fe_functions.h"
#include "elmt.h"
#include "miscellaneous.h"
/*
#define DEBUG
#define DEBUG1
#define DEBUG3
*/
/* ============================================================== */
/* Element SHELL_FOUR_NODES */
/* Implicit code */
/* Shell Element: */
/* material properties array */
/* Input Properties: */
/* ============================================================== */
/* p->work_material[0] = E;
p->work_material[1] = G;
p->work_material[2] = fy;
p->work_material[3] = ET;
p->work_material[4] = nu;
p->work_material[5] = density;
p->work_material[6] = fu;
p->work_section[0] = Ixx;
p->work_section[1] = Iyy;
p->work_section[2] = Izz;
p->work_section[3] = Ixy;
p->work_section[4] = Ixz;
p->work_section[5] = Iyz;
p->work_section[6] = weight;
p->work_section[7] = bf;
p->work_section[8] = tf;
p->work_section[9] = depth;
p->work_section[10] = area;
p->work_section[11] = thickness; */
/* ============================================================== */
ARRAY *elmt_shell(p, isw)
ARRAY *p;
int isw;
{
#ifdef DEBUG
printf(" enter elmt_shell() \n");
#endif
p = elmt_shell_4nodes_implicit(p, isw);
#ifdef DEBUG
printf(" Leaving elmt_shell() \n");
#endif
return(p);
}
ARRAY *elmt_shell_4nodes_implicit(p, isw)
ARRAY *p;
int isw;
{
static double nu;
static QUANTITY fy, G, E, ET, density, lambda;
static double Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
static double bf, tf, A, depth, h, EA, EIzz;
static DIMENSIONS *dp_length, *dp_force, *dp_moment;
static DIMENSIONS *dp_stress, *dp_degree, *dp_temperature;
double x21, x31, x42;
double y21, y31, y42;
double z21, z31, z42;
double co_x31, co_y31, co_z31;
double co_x42, co_y42, co_z42;
double EE, Jac,jacobian, k_bar = 0.833; /* k_bar is shear factor */
double elmt_length, aspect_ratio;
double **co_coord; /* co_coord is a coordinate rotated with material deformation */
double **T_matrix, **T_Transpose;
double **Direction_Matrix;
static double *integ_coord, *weight;
double **stress, **displ;
double **nodal_load, **nodal_load_temp;
double diff, sum, temp;
double **Cep, **stiff, **mass;
int i, j, k, ii, jj, kk, n, n1, n2, k1, k2;
int surface_pts, surf_integ_pts, no_integ_pts;
int size, dof, length, length1, length2;
int UNITS_SWITCH;
#ifdef DEBUG
printf("*** Enter elmt_shell_4nodes_implicit() : isw = %4d\n", isw);
#endif
/* =======================================================*/
/* Calculation of Direction Matrix; Transform coordinate */
/* into lamina coordiante system */
/* =======================================================*/
UNITS_SWITCH = CheckUnits();
co_coord = MatrixAllocIndirectDouble(p->no_dimen, p->nodes_per_elmt);
Direction_Matrix = MatrixAllocIndirectDouble(6,6);
for( i = 1; i <= 3; i++ ) {
for ( j = 1; j <= p->nodes_per_elmt; j++) {
co_coord[i-1][j-1] = p->coord[i-1][j-1].value; /* set initial co_coord */
}
}
Lamina_Sys_Implicit(p, Direction_Matrix, co_coord);
#ifdef DEBUG
dMatrixPrint("co_coord", co_coord, 3,4);
dMatrixPrint("Direction_Matrix", Direction_Matrix, 6,6);
#endif
x21 = p->coord[0][1].value - p->coord[0][0].value;
y21 = p->coord[1][1].value - p->coord[1][0].value;
z21 = p->coord[2][1].value - p->coord[2][0].value;
x31 = p->coord[0][2].value - p->coord[0][0].value;
y31 = p->coord[1][2].value - p->coord[1][0].value;
z31 = p->coord[2][2].value - p->coord[2][0].value;
x42 = p->coord[0][3].value - p->coord[0][1].value;
y42 = p->coord[1][3].value - p->coord[1][1].value;
z42 = p->coord[2][3].value - p->coord[2][1].value;
co_x31 = co_coord[0][2] - co_coord[0][0];
co_y31 = co_coord[1][2] - co_coord[1][0];
co_z31 = co_coord[2][2] - co_coord[2][0];
co_x42 = co_coord[0][3] - co_coord[0][1];
co_y42 = co_coord[1][3] - co_coord[1][1];
co_z42 = co_coord[2][3] - co_coord[2][1];
Jac = (x31*y42-x42*y31);
no_integ_pts = p->integ_ptr->thickness_pts;
dof = p->dof_per_node;
size = p->size_of_stiff;
/*****************************************************************/
#ifdef DEBUG
printf(" ***** In elmt_shell_4nodes_implicit(): \n");
printf(" : enter main switch \n");
#endif
switch(isw) {
case PROPTY:
#ifdef DEBUG
printf(" In elmt_shell_4nodes_implicit(): \n");
printf(" : enter case of PROPTY\n");
#endif
/* material properties: elastic */
E.value = p->work_material[0].value;
fy.value = p->work_material[2].value;
ET.value = p->work_material[3].value;
nu = p->work_material[4].value;
density.value = p->work_material[5].value;
lambda.value = nu*E.value/(1.0+nu)/(1.0-2.0*nu);
/* (1) check possion_ratio value */
if( nu == 0.0 || nu > 0.5 ) {
printf("WARNING >> ... In shell_Belytschko element() - nu value = %9.4f,reset to 0.3 !\n",
nu);
nu = 0.3; /* default poi_ratio value */
}
/* (2) calculate G value */
G.value = p->work_material[1].value = E.value/(1.0 - 2.0*nu) ;
if(E.value/((1.0 - 2.0*nu)) != p->work_material[1].value) {
printf(" elmt_shell(): WARNING: G is not equal to E/(1-2nu), check G for homogeneous material \n");
printf(" : ignore this message for non-homogeneous materials \n");
}
if( UNITS_SWITCH == ON ) {
E.dimen = p->work_material[0].dimen;
fy.dimen = p->work_material[2].dimen;
ET.dimen = p->work_material[3].dimen;
density.dimen = p->work_material[5].dimen;
lambda.dimen = E.dimen;
G.dimen = E.dimen;
}
Ixx = p->work_section[0].value;
Iyy = p->work_section[1].value;
Izz = p->work_section[2].value;
Ixy = p->work_section[3].value;
Ixz = p->work_section[4].value;
Iyz = p->work_section[5].value;
bf = p->work_section[7].value;
tf = p->work_section[8].value;
depth = p->work_section[9].value;
A = p->work_section[10].value;
h = p->work_section[11].value; /* thickness of the shell */
EA = E.value*A;
EIzz = E.value*Izz;
#ifdef DEBUG
printf(" In elmt_shell_4nodes_implicit(): \n");
printf(" : leaving case of PROPTY\n");
#endif
break;
case STIFF: /* form element stiffness */
#ifdef DEBUG
printf("*** In elmt_shell() : start case STIFF\n");
printf(" : Density = %14.4e\n", density.value);
printf(" : shell_thickness = %8.2f\n", h);
printf(" : Jac = %8.2f\n", Jac);
printf(" : E = %lf\n", E.value);
printf(" : nu = %lf\n", nu);
printf(" : no_integ_pts in z_direction = %d\n", no_integ_pts);
#endif
/* Integration loop over thickness */
stiff = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
T_matrix = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
T_Transpose = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
integ_coord = dVectorAlloc(no_integ_pts + 1);
weight = dVectorAlloc(no_integ_pts + 1);
size = p->size_of_stiff;
gauss(integ_coord,weight,no_integ_pts);
/* Integration loop over surface */
for (ii = 1; ii <= no_integ_pts; ii++) {
Shell_Stiff_Plane_4node(stiff, p, co_coord, integ_coord[ii], ii, E.value, nu);
for ( i = 1; i <= size; i++) {
for ( j = 1; j <= size; j++) {
p->stiff->uMatrix.daa[i-1][j-1] += stiff[i-1][j-1]*h*0.5*weight[ii];
}
}
} /* end of gaussian integration through thickness */
#ifdef DEBUG
printf(" End of integration through the thickness \n");
#endif
#ifdef DEBUG
dMatrixPrint(" stiff before rotation", p->stiff->uMatrix.daa, p->size_of_stiff, p->size_of_stiff);
#endif
free((char *) integ_coord);
free((char *) weight);
/***************************************************/
/* Transform stiffness matrix from local lamina */
/* coordinate system to global coordinate system */
/***************************************************/
/* store the direction matrix : T_matrix */
for (i = 1; i <= p->size_of_stiff; i++) {
for(j = 1; j <= p->nodes_per_elmt; j++) {
for( k = 1; k <= p->dof_per_node; k++) {
n1 = p->dof_per_node*(j-1);
n2 = p->dof_per_node*j;
n = n1 + k;
if(i > n1 && i <= n2) {
ii = i - n1;
T_matrix[i-1][n-1] = Direction_Matrix[ii-1][k-1];
}
else
T_matrix[i-1][n-1] = 0.0;
}
}
}
stiff = dMatrixMultRep(stiff, p->stiff->uMatrix.daa, size, size, T_matrix, size, size);
for(i = 1; i <= size; i++){
for(j = 1; j <= size; j++) {
T_Transpose[i-1][j-1] = T_matrix[j-1][i-1];
}
}
p->stiff->uMatrix.daa = dMatrixMultRep(p->stiff->uMatrix.daa,
T_Transpose, size, size, stiff, size, size);
#ifdef DEBUG
/* check the symmetry of the stiffness matrix */
for(i = 1; i <= size; i++){
for(j = 1; j <= size; j++) {
diff = p->stiff->uMatrix.daa[i-1][j-1] -p->stiff->uMatrix.daa[j-1][i-1];
if(diff > 1E-7 &&
(ABS(p->stiff->uMatrix.daa[i-1][j-1])) > 1E-7 &&
(ABS(diff/p->stiff->uMatrix.daa[i-1][j-1])) > 1E-1) {
printf("K[%d][%d] = %lf \n",
i, j, p->stiff->uMatrix.daa[i-1][j-1]);
printf("K[%d][%d] = %lf \n",
j, i, p->stiff->uMatrix.daa[j-1][i-1]);
printf("diff = %le \n",
p->stiff->uMatrix.daa[i-1][j-1] -
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -