Wikidot-ing!

December 29, 2007

Remember the C-source code for calculating the elastic stress fields of a circular cavity in a square plate under uniaxial stress? I have found another neat way of shipping code with comments (and equations and figures, if need be) thanks to Abi: here is the wikidot page of code, schematic and equations. Have fun!

The wonderful guys at WordPress have given a wraparound so that posting source code is easier; it can’t get any better than this. Look out these pages for more phase field codes! For now, as a test, here is the code which will calculate the elastic stress fields around a circular hole in a square plate which is under an uniaxial stress.


/************
This program calculates the elastic stress fields of a circular hole
in a square plate under an uniaxial stress. The stress is applied
along the x-axis. For the expressions used in calculating the stress
fields, see p. 109 of Elasticity by Barber -- Equations 8.74, 8.75, and 8.76.
************/

#include 
#include
#include

int main(void){

FILE *fp1, *fp2, *fp3;

double S; /* Applied uniaxial stress */
double a; /* Radius of the circular hole */
double r; /* Distance along the x- or y-axis */

int i;
double s11, s12, s22; /* The 11, 12 and 22 stress fields respectively */

S = 1.0;
a = 2.5;

fp1 = fopen("sigma11_x","w");
fp2 = fopen("sigma22_x","w");
fp3 = fopen("sigma12_x","w");

for(i=0;i<3000; ++i){
r = 0.01*i;
if (r < a){
s11 = s12 = s22 = 0.0;
}
else{
s11 = 0.5*S*(1. - a*a/(r*r)) + 0.5*S*(3.*a*a*a*a/(r*r*r*r) - 4.*a*a/(r*r) + 1.);
s22 = 0.5*S*(1. + a*a/(r*r)) - 0.5*S*(3.*a*a*a*a/(r*r*r*r) + 1.);
s12 = 0.0;
}
fprintf(fp1,"%le %le\n",r,s11);
fprintf(fp2,"%le %le\n",r,s22);
fprintf(fp3,"%le %le\n",r,s12);

}

fclose(fp1);
fclose(fp2);
fclose(fp3);

fp1 = fopen("sigma11_y","w");
fp2 = fopen("sigma22_y","w");
fp3 = fopen("sigma12_y","w");

for(i=0;i<3000; ++i){
r = 0.01*i;
if (r < a){
s11 = s12 = s22 = 0.0;
}
else{
s22 = 0.5*S*(1. - a*a/(r*r)) - 0.5*S*(3.*a*a*a*a/(r*r*r*r) - 4.*a*a/(r*r) + 1.);
s11 = 0.5*S*(1. + a*a/(r*r)) + 0.5*S*(3.*a*a*a*a/(r*r*r*r) + 1.);
s12 = 0.0;
}
fprintf(fp1,"%le %le\n",r,s11);
fprintf(fp2,"%le %le\n",r,s22);
fprintf(fp3,"%le %le\n",r,s12);

}

fclose(fp1);
fclose(fp2);
fclose(fp3);

return 0;
}

At the first go, looks like everything works fine except for (a) the #include things–for which, the stuff in triangular brackets, it seems to consider as html commands, and (b) the greater then signs, which it seem to have some problems showing as greater then signs. Code within the sourcecode wraparound are supposed to be formatted automatically. The language option I chose is cpp, which I assume is C++. Could that be the problem?

Test: <>

 #include

Hmm…I guess the wraparound requires a bit of fixing!!!

Test 2:

$features = file_get_contents( 'http://wordpress.com/features/' );
preg_match_all( '|<h3>(.*?)</h3>|is', $features, $why_wp_rocks );
foreach ( $why_wp_rocks[1] as $slick_feature )
$hotness[] = $slick_feature;
var_dump( $hotness );

If you are writing a code for solving for the equations of mechanical equilibrium in an elastically inhomogeneous system, then there are several test cases that you might wish to run to make sure that your code is working fine. One of them is the stress fields of a circular hole in a plate which is under an applied uniaxial stress. Typically, the problem is solved in polar coordinates and the stresses at any given point (r,\theta) are given by the following expressions; see Elasticity (2nd Edition) by J R Barber, p. 109, Equations 8.74 — 8.76, for example — by the way, a sample chapter from the book is available for download here (pdf).

\sigma_{rr} = \frac{S}{2} \left( 1 - \frac{a^{2}}{r^{2}}\right) + \frac{S \cos(2\theta)}{2}\left(\frac{3 a^{4}}{r^{4}} - \frac{4 a^{2}}{r^{2}} + 1 \right) ;

\sigma_{r \theta} =  \frac{S \sin(2\theta)}{2}\left(\frac{3 a^{4}}{r^{4}} - \frac{2 a^{2}}{r^{2}} - 1 \right) ;

\sigma_{\theta \theta} = \frac{S}{2} \left( 1 + \frac{a^{2}}{r^{2}}\right) - \frac{S \cos(2\theta)}{2}\left(\frac{3 a^{4}}{r^{4}} + 1 \right) ;

where S is applied stress, and a is the radius of the cavity.

Here is a code written in C which will generate the data files for the stresses along the x– and y-axes from the center of the circular cavity of unit radius under an applied (tensile) uniaxial stress of unity. Note that you should have write permissions in the directory for the output data files to be written. Compile the code (gcc –lm circ_hole_under_uniaxial_stress.c) and execute the resultant binary a.out.

You may have to tweak the code a bit if you have different applied stress and/or cavity radius. The data files generated using the code can then be plotted using gnuplot; and, such plots of the stress fields are shown here: Stress fields along x-axis(pdf) and Stress fields along y-axis (pdf). The following aspects of the plots are worth noting:

  1. The \sigma_{12} components of the stress fields are zero both along x– and y-axes;
  2. When the applied stress is along the x-direction, the \sigma_{11} component of the stress at the cavity along the y-axis jumps to thrice the applied stress; and,
  3. Far away from the cavity, the stress fields in the plate are just the applied fields (as one would expect from St. Venant principle, which was used to obtain the stress fields in the first place).

Have fun!