📄 elmt_shell_8n.c
字号:
free((char *) dp_force->units_name);
free((char *) dp_force);
free((char *) dp_length->units_name);
free((char *) dp_length);
free((char *) dp_moment->units_name);
free((char *) dp_moment);
break;
case OFF:
break;
default:
break;
}
#ifdef DEBUG
printf("*** In elmt_shell() : print equiv_nodal_load\n\n");
MatrixPrintIndirectDouble(p->equiv_nodal_load);
printf("*** In elmt_shell() : end case EQUIV_NODAL_LOAD\n");
#endif
break;
case STRESS_UPDATE:
/* update stress for given displacement */
/* and displacement incremental */
break;
case STRESS:
case STRESS_LOAD:
case LOAD_MATRIX:
/* ====================================== */
/* Form internal nodal load vector due to */
/* stress at previous load step */
/* ====================================== */
#ifdef DEBUG
printf("****** In elmt_shell_8nodes_implicit() : \n");
printf(" enter cases: STRESS_UPDATE, LOAD_MATRIX, STRESS_LOAD \n");
printf(" elemt_state = %d \n", p->elmt_state);
#endif
nodal_load = MatrixAllocIndirectDouble(p->size_of_stiff, 1);
nodal_load_temp = MatrixAllocIndirectDouble(p->size_of_stiff, 1);
h = p->work_section[11].value; /* thickness of the shell */
EE = E.value;
/* =================================================*/
/* nodal load due to the inital stress or stress at */
/* previous step */
/* =================================================*/
/* Integration loop over thickness */
integ_coord = dVectorAlloc(no_integ_pts+1);
weight = dVectorAlloc(no_integ_pts+1);
gauss(integ_coord,weight,no_integ_pts);
for(ii = 1; ii <= no_integ_pts; ii++) {
nodal_load_temp
= Shell_Nodal_Load_8Node(nodal_load_temp, p, co_coord,
integ_coord[ii], ii, E.value, nu, isw);
for(i = 1; i<= p->size_of_stiff; i++) {
nodal_load[i-1][0] += nodal_load_temp[i-1][0]*weight[ii];
}
} /* gaussian integration ends */
free((char *) integ_coord);
free((char *) weight);
size = p->size_of_stiff;
for(i= 1; i<= p->size_of_stiff; i++) {
p->nodal_loads[i-1].value = nodal_load[i-1][0];
}
/***************************************************/
/* Transform nodal_load vector from local lamina */
/* coordinate system to global coordinate system */
/***************************************************/
/* =====================================*/
/* Calculate the inverse of T_matrix */
/* T_matrix inverse [T]^-1 = */
/* = T_matrix trnspose : [T]^t */
/* =====================================*/
/* store the direction matrix : T_matrix */
if( isw == LOAD_MATRIX ) {
e1_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
e2_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
e3_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
Lamina_Sys_8node(p, e1_ptr, e2_ptr, e3_ptr);
Direction_Matrix = MatrixAllocIndirectDouble(6,6);
for ( i = 1; i <= 6; i++) {
if(i <= 3) {
Direction_Matrix[0][i-1] = e1_ptr[i-1][0];
Direction_Matrix[1][i-1] = e2_ptr[i-1][0];
Direction_Matrix[2][i-1] = e3_ptr[i-1][0];
Direction_Matrix[3][i-1] = 0.0;
Direction_Matrix[4][i-1] = 0.0;
Direction_Matrix[5][i-1] = 0.0;
}
else {
j = i - 3;
Direction_Matrix[0][i-1] = 0.0;
Direction_Matrix[1][i-1] = 0.0;
Direction_Matrix[2][i-1] = 0.0;
Direction_Matrix[3][i-1] = e1_ptr[j-1][0];
Direction_Matrix[4][i-1] = e2_ptr[j-1][0];
Direction_Matrix[5][i-1] = e3_ptr[j-1][0];
}
}
T_matrix = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
T_Transpose = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
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;
}
}
}
for( i = 1; i <= p->size_of_stiff; i++) {
for( j = 1; j <= p->size_of_stiff; j++)
T_Transpose[i-1][j-1] = T_matrix[j-1][i-1];
}
size = p->size_of_stiff;
nodal_load_temp = dMatrixMultRep(nodal_load_temp,
T_Transpose, size, size, nodal_load, size,1);
for(i= 1; i<= p->size_of_stiff; i++) {
p->nodal_loads[i-1].value = nodal_load_temp[i-1][0];
}
MatrixFreeIndirectDouble(e1_ptr,3);
MatrixFreeIndirectDouble(e2_ptr,3);
MatrixFreeIndirectDouble(e3_ptr,3);
MatrixFreeIndirectDouble(Direction_Matrix,6);
MatrixFreeIndirectDouble(T_matrix,p->size_of_stiff);
MatrixFreeIndirectDouble(T_Transpose,p->size_of_stiff);
}
MatrixFreeIndirectDouble(nodal_load,p->size_of_stiff);
MatrixFreeIndirectDouble(nodal_load_temp,p->size_of_stiff);
/* ------------NODAL LOAD UNITS ------------------------*/
switch( UNITS_SWITCH ) {
case ON:
if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
dp_force = DefaultUnits("N");
dp_length = DefaultUnits("m");
}
else {
dp_force = DefaultUnits("lbf");
dp_length = DefaultUnits("in");
}
/* node no 1 */
UnitsCopy( p->nodal_loads[0].dimen, dp_force );
UnitsCopy( p->nodal_loads[1].dimen, dp_force );
UnitsCopy( p->nodal_loads[2].dimen, dp_force );
UnitsMultRep( p->nodal_loads[3].dimen, dp_force, dp_length );
UnitsCopy( p->nodal_loads[4].dimen, p->nodal_loads[3].dimen );
/* node no > 1 */
for(i = 2; i <= p->nodes_per_elmt; i++) {
for(j = 1; j <= p->dof_per_node; j++) {
k = p->dof_per_node*(i-1)+j;
if(j <= 3)
UnitsCopy( p->nodal_loads[k-1].dimen, p->nodal_loads[0].dimen );
else
UnitsCopy( p->nodal_loads[k-1].dimen, p->nodal_loads[3].dimen );
}
}
free((char *) dp_force->units_name);
free((char *) dp_force);
free((char *) dp_length->units_name);
free((char *) dp_length);
break;
case OFF:
break;
default:
break;
}
#ifdef DEBUG
printf("*** In elmt_shell_8nodes_implicit() : end case LOAD_MATRIX, STRESS_UPDATE, STRESS_LOAD \n");
#endif
break;
case MASS_MATRIX: /* form mass matrix */
#ifdef DEBUG
printf("*** In elmt_shell_8nodes_implicit() : start case MASS\n");
printf(" : Density = %14.4e\n", density.value);
printf(" : shell_thickness = %8.2f\n", h);
#endif
/*========================== */
/* CALCULATING MASS MATRIX */
/*========================== */
Shell_8Node_Mass(p, p->stiff, co_coord, density.value, h);
#ifdef DEBUG
printf("*** In elmt_shell() : end case MASS\n");
#endif
break;
default:
break;
}
MatrixFreeIndirectDouble(co_coord, p->no_dimen);
#ifdef DEBUG
printf("*** leaving elmt_shell() \n");
#endif
return(p);
}
/* ============================================*/
/* Equivalent Load due to distributed pressure */
/* for 8 node shell element */
/* ============================================*/
ARRAY *sld108(p,task)
ARRAY *p;
int task;
{
#ifdef DEBUG
printf("**** enter sld108(): \n");
#endif
/* add details */
#ifdef DEBUG
printf("**** leaving sld108(): \n");
#endif
return(p);
}
void Stress_Update_8Node(p, co_coord, B1_matrix, integ_pt)
ARRAY *p;
double **co_coord;
double **B1_matrix;
int integ_pt;
{
double H, R, fy, E, Et, nu;
double **displ_incr,**m1,**strain_incr;
double **back_stress_incr, **strain_pl;
double **stress, **stress_dev, **stress_incr;
double temp, G, A_2, A, Lambda, mean_stress;
double eff_pl_strainTemp, eff_pl_strain;
double Beta1 = 0.05, temp1, effect_stress;
int iNo_iter_step, dof, size;
int i, j, k, ii, jj, kk, n;
DIMENSIONS *dimen;
/* Given i -state displ, stress mater_matrix, [Cep], */
/* and strain, calculate i+1 state displ_incr */
/* these i+1 state stress can be used for */
/* i+1 state [Cep] matrix */
#ifdef DEBUG
printf("\n enter Stress_Update_8Node(): \n");
#endif
kk = integ_pt;
dof = p->dof_per_node;
size = p->size_of_stiff;
E = p->work_material[0].value;
nu = p->work_material[4].value;
fy = p->work_material[2].value;
G = E/(1.0+nu)/2.0;
m1 = MatrixAllocIndirectDouble(p->dof_per_node, p->dof_per_node);
m1 = MATER_MAT_SHELL(m1, E, nu);
/* Step 1: calculate strain_incremental */
/* from displacement incremental */
/* at given integration point */
displ_incr = MatrixAllocIndirectDouble(p->size_of_stiff, 1);
strain_incr = MatrixAllocIndirectDouble(dof,1);
strain_pl = MatrixAllocIndirectDouble(dof,1);
stress = MatrixAllocIndirectDouble(dof,1); /* stress or stress incremental */
stress_dev = MatrixAllocIndirectDouble(dof,1); /* deviatric component of stress */
stress_incr = MatrixAllocIndirectDouble(dof,1); /* stress incremental */
back_stress_incr = MatrixAllocIndirectDouble(dof,1); /* stress incremental */
#ifdef DEBUG
printf(" In Stress_Update_8Node(): materials_name = %s \n",p->material_name);
#endif
/* ================= */
/* Elastic material */
/* ================= */
if(p->material_name != NULL && !strcmp(p->material_name,"ELASTIC")) {
/* Transfer p->displ into vector form */
/* -----------------------------------*/
/* There is no needed to calculate */
/* incrementally for the case of */
/* linear elasticity */
/* -----------------------------------*/
for(i = 1; i <= p->dof_per_node; i++) {
for(j = 1; j <= p->nodes_per_elmt; j++) {
k = p->dof_per_node*(j-1)+i;
displ_incr[k-1][0] = p->displ->uMatrix.daa[i-1][j-1];
}
}
strain_incr = dMatrixMultRep(strain_incr, B1_matrix, dof,
size, displ_incr, size, 1);
stress = dMatrixMultRep(stress, m1, dof, dof,
strain_incr, dof, 1);
#ifdef DEBUG
dMatrixPrint("displ_incr", displ_incr, size, 1);
dMatrixPrint("strain_incr", strain_incr, dof, 1);
dMatrixPrint("stress", stress, dof, 1);
#endif
for (i = 1; i <= dof; i++){
p->stress->uMatrix.daa[i-1][kk-1] += stress[i-1][0];
}
SaveRespondBuffer(p, kk);
}
/* ====================================*/
/* Elastic Plastic & Elastic Perfectly */
/* Plastic Materials */
/* ====================================*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -