این برنامه جهت محاسبه تنش و کرنش و همچنین تحلیل تخریب در کامپوزیت‌ پلیمری لایه های درنرم افزار MATLAB 7.12 نوشته شده است. جهت کد نویسی از تئوری کلاسیک صفحات کامپوزیتی استفاده شده است. ابتدا خواص ماده، ترتیب لایه چینی، تعداد لایه‌ها، زاویه لایه‌ها و ضخامت هر لایه کامپوزیت به‌عنوان ورودی وارد می‌شود. سپس ماتریس‌های B،A و D کامپوزیت محاسبه می‌گردند. بارگذاری روی لایه چینی می‌تواند دارای 3 مولفه نیرو و 3 مولفه تنش باشد که به‌عنوان ورودی وارد می‌شود. با توجه به بارگذاری وارد شده، در ابتدا کرنش‌های لایه ای (Laminate) در جهات بارگذاری محاسبه می‌شوند. کرنش در هر لایه کامپوزیت با استفاده از کرنش کلی و انحنا به دست می‌آید. با یک دوران، کرنش در هر لایه از جهات فرعی به جهات اصلی (مادی) انتقال داده می‌شود. با استفاده از روابط تنش و کرنش، تنش‌ها در هر لایه محاسبه می‌شوند.

محاسبه تنش، کرنش و تخریب در کامپوزیت (متلب)

پس از محاسبه تنش‌ها و کرنش‌ها در هر لایه کامپوزیت، می‌توان از تئوری‌های تخریب برای لایه چینی استفاده نمود. در این کد از تئوری‌های تنش ماکزیمم، کرنش ماکزیمم، سای هیل (Tsai-Hill) و سای وو (Tsai-Wu) استفاده شده است. برای هر کدام از این تئوری‌ها یک ضریب ایمنی توسط کد ارائه می‌شود. به‌طوری‌که اگر تمام بارهای وارده در این مقدار ضرب شوند، تخریب آغاز می‌شود. بنابراین می‌توان محاسبه نمود که چه میزان نیرو یا ممان وارده باعث آغاز تخریب خواهند شد. طبیعی است که اگر ضریب ایمنی از 1 کمتر باشد، تخریب رخ‌داده است.

اثرات حرارت و رطوبت نیز در محاسبه تنش‌ها لحاظ شده‌اند. بدین معنی که تغییرات در درجه حرارت و رطوبت، کرنش‌ها و تنش‌های جدیدی را ایجاد می‌کنند (تنش‌های پسماند) که در کنار تنش‌ها و کرنش‌های مکانیکی (تحت بارگذاری) روی روند تخریب نقش دارند. در شرایطی که تنش‌های پسماند موجود باشند، ضرایب انبساط حرارتی و رطوبتی و همچنین تغییرات در درجه حرارت و رطوبت نیز به‌عنوان ورودی درنرم افزار وارد می‌گردند.

شرح قسمت‌های مختلف کد

در این بخش قسمت‌های مختلف کد به‌صورت جداگانه بررسی می‌شود.

دریافت خواص مواد در این قسمت خواص کامپوزیت گرافیت اپوکسی به‌عنوان پیش فرض بهنرم افزار داده شده است.

clc
clear all
disp(‘Analysis of Composite Materials’)
% Material Constants (for Graphite/Epoxy)
Ex=181; %GPa
Ey=10.3; %GPa
NUx=.28;
NUy=(Ey/Ex)*NUx;
Es=7.17; %GPa

محاسبه ماتریس‌های نرمی و سفتی مادی

% S & Q Matrix (On axis compiance & Stiffness)
S_on=zeros(3,3);
S_on(1,1)=1/Ex;
S_on(1,2)=-NUy/Ey;
S_on(2,1)=S_on(1,2);
S_on(2,2)=1/Ey; S_on(3,3)=1/Es; Q_on=inv(S_on); disp(”)

ورود خواص لایه چینی

زاویه لایه‌ها به‌صورت یک بردار بهنرم افزار داده می‌شود. همچنین تعداد لایه‌های موجود نیز به‌صورت یک بردار با طول برابر با بردار لایه چینی وارد می‌گردد.

laminate=[0 90 0];
% Angles
number=[1 1 1];
% # of each ply
t=5e-3;
%Thickness of each ply
h=sum(number)*t;
% Total thickness

محاسبه ماتریس‌های B،A و D

در این کد برای محاسبه ماتریس‌های B،A و D از روش محاسباتی ارائه شده در کتاب‌های مرجعی نظیر کاو وسای استفاده شده است. این روش به دلیل وجود محاسبات ماتریسی کمتر دارای سرعت بالاتری برای تعداد لایه‌های زیاد می‌باشد. جزئیات روابط در مراجع ذکر شده موجود است. لازم به ذکر است که برای نمایش بهتر ماتریس‌ها (و همچنین تنش‌ها) توسط ضرایبی تغییر واحدهای مناسبی در کد انجام شده است.

بیشتر بخوانید: مدل سازی مواد مركب (کامپوزیت) در Abaqus

ورود ضرایب حرارتی و رطوبتی، تغییر در دما و رطوبت

a1=.02e-6; a2=22.5e-6; % Thermal expansion coefficients
a3=0; alfaon=[a1;a2;a3];
b1=0; b2=.6; %Moisture Coefficients
b3=0; betaon=[b1;b2;b3];
teta=0; % Temprature Change delta=0; % Moisture Change NT=[0;0;0]; %Stress resultant due to Temprature
MT=[0;0;0]; %Moment resultant due to Temprature
NB=[0;0;0]; %Stress resultant due to Moisture
MB=[0;0;0]; %Moment resultant due to Temprature

محاسبه نیروهای حرارتی و رطوبتی

for i=1:length(laminate)
m(i)=cos(laminate(i)*pi/180);
n(i)=sin(laminate(i)*pi/180);
alfa1(i)=a1*m(i)^2+a2*n(i)^2;
alfa2(i)=a2*m(i)^2+a1*n(i)^2;
alfa3(i)=2*(a1-a2)*m(i)*n(i);
beta1(i)=b1*m(i)^2+b2*n(i)^2;

%Tranforming thermal coefficients to off axis
beta2(i)=b2*m(i)^2+b1*n(i)^2;
beta3(i)=2*(b1-b2)*m(i)*n(i);

%Tranforming Moisture coefficients to off axis
T=[m(i)^2 n(i)^2 2*m(i)*n(i);
n(i)^2 m(i)^2 -2*m(i)*n(i);
m(i)*n(i) m(i)*n(i) m(i)^2-n(i)^2];
%Trans. Matrix
Rm=[1 0 0;0 1 0;0 0 2];
Qoff=inv(T)*Q_on*Rm*T*inv(Rm);
%off-axis Stiffness alfa=[alfa1(i);
alfa2(i);alfa3(i)];

%off-axis thermal coefficients
NT=NT+Qoff*1e9*alfa*ttt*teta;
MT=MT+.5*Qoff*1e9*alfa*teta*(zzz(i+1)^2-zzz(i)^2);
beta=[beta1(i);beta2(i);beta3(i)];

%off-axis thermal coefficients
NB=NB+Qoff*1e9*beta*ttt*teta;
MB=MB+.5*Qoff*1e9*beta*delta*(zzz(i+1)^2-zzz(i)^2);
end

ورود بارگذاری مکانیکی و جمع آن با بارهای حرارتی و رطوبتی تبدیل واحد برای نیروها و ممانها انجام شده است.

% Input Loadings, input N1,N2,N6 in N/m
N1=1; %N/m
N2=0; %N/m
N6=0; %N/m
N1=N1*1e-9; %GN/m
N2=N2*1e-9; %GN/m
N6=N6*1e-9; %GN/m

%Input M1,M2,M6 in N
M1=0;
M2=0;
M6=0;
M1=M1*1e-6;
M2=M2*1e-6;
M6=M6*1e-6;
N = [N1; N2; N6];
M = [M1; M2; M6];

%Adding the effect of 5emperature & Moisture
N=N+NT*1e-6+NB*1e-6;
M=M+MT*1e-6+MB*1e-6;

محاسبه تنش‌ها و کرنش‌ها در هر لایه

eps0_k = Compliance * [N; M];
eps0 = eps0_k(1:3);

% Midplane Strains
k= eps0_k(4:6);

% Midplane curvatures
for i=1:length(laminate)
zz(i)=.5*(z(i)+z(i+1));
%midplane of each ply
end

%Calculting strain in each ply for i=1:length(laminate) ep1(i)=eps0(1)+k(1)*zz(i);
ep2(i)=eps0(2)+k(2)*zz(i);
ep6(i)=eps0(3)+k(3)*zz(i);
end

%Calculating stress in each ply for i=1:length(laminate) m(i)=cos(laminate(i)*pi/180);
n(i)=sin(laminate(i)*pi/180);

T_eps=[m(i)^2 n(i)^2 m(i)*n(i);
n(i)^2 m(i)^2 -m(i)*n(i);-2*m(i)*n(i)

2*m(i)*n(i) m(i)^2-n(i)^2]; eps_on=T_eps*[ep1(i);ep2(i);ep6(i)];
Sigma_on=Q_on*(eps_on-alfaon*teta-betaon*delta);
Sigma_X(i)=Sigma_on(1); %GPa

Sigma_Y(i)=Sigma_on(2); %GPa

%On axis Stresses
Sigma_S(i)=Sigma_on(3); %GPa
e1(i)=eps_on(1);
e2(i)=eps_on(2);
e6(i)=eps_on(3);
end

disp(‘Stress in each layer (Pa): ‘)

for i=1:length(laminate)
disp(‘Layer’);
disp(i);
disp(‘Sigma_X’)
1e9*Sigma_X(i)
disp(‘Sigma_Y’)
1e9*Sigma_Y(i)
disp(‘Sigma_S’)
1e9*Sigma_S(i)
end

ورود خواص استحکامی ماده بهنرم افزار

در این مرحله 4 معیار تنش و کرنش ماکزیمم، تسای هیل و تسای وو نمایش داده شده و از کاربر درخواست می شود معیار موردنظر خود را انتخاب کند.

disp(‘ 1- Maximum On axis Stress Critera’)
disp(‘ 2- Tsi_Wu Criterian ‘)
disp(‘ 3- Tsi-Hill Criteria ‘)
disp(‘ 4- Maximum Strain ‘)
disp(‘ ‘)

n=input(‘Enter the number of Analysis: 1,2,3 or 4? ‘);
disp(‘ ‘)

X_T=input(‘Enter long. tensile strength in MPa: X_T? ‘)*1e-3;

X_C=input(‘Enter long. compressive strength in MPa: X_C? ‘)*1e-3; %MPa to GPa

Y_T=input(‘Enter trans. tensile strength in MPa: Y_T? ‘)*1e-3;

Y_C=input(‘Enter trans. compressive strength in MPa: Y_T? ‘)*1e-3;

S=input(‘Enter shear strength S? ‘)*1e-3;

بررسی تخریب با استفاده از هر کدام از تئوری‌ها، ضریب ایمنی محاسبه‌شده و اعلام می‌گردد.

حل مثال نمونه: مسئله 5.3 کتاب برای لایه چینی سه لایه s [09/0] گرافیت اپوکسی، مقادیر به دست آمده از کد با نرم افزار مقایسه می‌شود. برای بارگذاری Nx=1N/m:

تنش‌ها در لایه 0 درجه:

ضریب ایمنی لایه اول با استفاده از سای-وو

ضریب ایمنی لایه اول با استفاده از کرنش ماکزیمم