哈工大 材料力学 MATLAB 梁 上机 大作业
哈工大材料力学大作业--matlab编程
H a r b i n I n s t i t u t e o f T e c h n o l o g y材料力学上机作业课程名称：材料力学设计题目：应力状态分析院系：机电学院班级：分析者：学号：指导教师：***设计时间：2013年6月18日哈尔滨工业大学材料力学上机课设计说明书一， 设计题目题目7 应力状态分析 输入：1. 平面应力状态输入：x y xy σστ（,,）；某截面方位角α2. 空间应力状态输入：,x y z xy yz zx σσστττ（,,,,）输出： 1. 输出主应力123σσσ（,,）2. 最大切应力（13max 132σσττ-==）3.如为平面应力状态则需要输出方位角α斜截面上的应力ααστ、及主方向角*σα4. 画出应力圆示意图二， 程序计算设计过程1. 平面应力状态分析对于任意平面应力状态，有max min σσ=2x y σσ+±主应力为：1max 23min ,0,σσσσσ===并且由 2tan 2xyx yστασσ=-可求得主应力方向角13σσαα、。
对于任意一个方位角α，有：=cos 2sin 222sin 2cos 22x yx yxy x yxy αασσσσσατασστατα++++-=-+从而，输入任意角α，即可求得该截面的应力状态ααστ、并且ααστ、都是关于α的函数，上式即为应力圆的参数方程，参数为α。
将α从0到pi 取一系列的值，则可以求出一系列的ααστ、，在坐标系中找到对应点，连接即可作出应力圆。
2. 三向应力状态分析解特征方程 321230I I I σσσ-+-=即可求出主应力123σσσ、、 其中：123||||||||x y z xyx y zy z xz xy y yz z zx x x yx zx xyy zy xzyz z I I I σσσστστσττστστσστττστττσ=-+⎛⎫⎛⎫⎛⎫=++⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎛⎫⎪= ⎪ ⎪⎝⎭再由 13max 132σσττ-== 可求得最大切应力。
材料力学上机大作业(哈工大)
H a r b i n I n s t i t u t e o f T e c h n o l o g y材料力学上机报告课程名称：材料力学设计题目：二向应力状态分析院系：XXXXXX班级：XXXXXX设计者：XXXXXX学号：XXXXXX设计时间：2013.06.18哈尔滨工业大学二向应力状态分析一：课题要求1.输入：任意一点的应力状态：（σx、σy、τxy）；某截面方位角α2.输出：输入点的主应力（σ1、σ2、σ3），方位角α斜截面上的应力σα、τα。
及主方向角α3.画出应力圆示意图。
4.程序运行时为界面显示形式。
二：程序框图三：所编程序x=str2double(get(handles.edit1,'string'));y=str2double(get(handles.edit2,'string'));xy=str2double(get(handles.edit3,'string'));M=str2double(get(handles.edit4,'string')); %将窗口输入值分别赋给x,y,xy,M b=sqrt((x/2-y/2)^2+xy^2);x1=(x+y)/2+b;x3=(x+y)/2-b;x2=0;if x1<0x2=x1;x1=0;endt=(x1-x3)/2;M=M*pi/180;b1=(x+y)/2+(x-y)*cos(2*M)/2-xy*sin(2*M);b2=(x-y)*sin(2*M)/2+xy*cos(2*M);b3=90*atan((-2*xy)/(x+y))/pi;%计算输出的主切应力大小、方向和截面上的应力并赋值set(handles.edit5,'string',x1);set(handles.edit6,'string',x2);set(handles.edit7,'string',x3);set(handles.edit9,'string',t);set(handles.edit10,'string',b3);set(handles.edit11,'string',b1);set(handles.edit12,'string',b2);%在输出窗口显示主切应力大小、方向和截面上应力b4=sqrt(b.^2+t.^2);v1=(x+y)/2-b4:0.001:(x+y)/2+b4;b11=sqrt(b4.^2-(v1-(x+y)/2).^2);b12=-sqrt(b4.^2-(v1-(x+y)/2).^2);%绘制应力圆上的点axes(handles.axes1); %选择应力圆的输出地址plot(v1,b11,v1,b12);grid on%绘制应力圆以上程序为在matlab中使用GUI编程时的主代码，界面代码请见m文件。
哈工大材料力学上机实验报告
材料力学I上机实验设计报告院系：机电学院班级： 1308***姓名：***学号： 11308*****指导教师：张桂莲时间：2015年6月一、问题描述1、应力状态分析对于空间或者是平面应力状态的相关计算，如果采用人工计算的方式比较繁琐而且容易出错，对于这种简单的重复计算，编制相应的程序则可以大大提高计算准确度和人工计算强度。
对于平面应力状态，输入量应为（,,x y xy σστ），以及某截面的方位角α，其输出数据应为该单元体所受主应力（123,,σσσ），所受最大剪应力（13max 132σσττ-==），以及方位角为α的斜截面上的应力（,ααστ）以及主方向角σα，同时还要画出其应力圆示意图，以直观的显示其应力状态。
对于空间应力状态，输入量则应该为各应力（,,,,,x y z xy yz xz σσστττ），其输出数据应该为该单元体所受主应力（123,,σσσ），所受最大剪应力（13max 132σσττ-==），同时还要画出其应力圆示意图，以直观的显示其应力状态。
这样，应力状态分析的基本任务就可以完成。
2、常用截面图形几何性质的分析在生活中，有各种各样的几何形状，但是对于工程实际中经常用到的构件，其截面的几何形状则非常有限。
对于不同的截面，其形心位置、对于形心轴的惯性矩也就有所不同，这样在进行如弯曲、扭转等的应力分析时就会到来不便，因此编制相应的程序来计算相关截面的几何性质也就具有了实际应用价值和可行性。
在这部分程序中，截面几何形状分为三角形、矩形、椭圆形、梯形、圆形、扇形等多种形式，对于不同的截面形状，输入量也就不同。
例如，对于扇形应输入直径和圆心角（,d α）；对于梯形则应输入上底、下底和高（,,a b h ）；对于椭圆形，则要输入长轴长和短轴长（,a b ）等等，在此不一一列举，具体输入数据请参看程序运行。
不过对于不同的截面，其输出的量都是相同的，即截面形心的位置、面积、对于形心轴的惯性矩（,,,,C C C C y z y z S I I ），这些输出量就是这些截面的基本几何性质参数，有了这些参数之后则可以对其进行进一步的计算和接下来的分析等问题。
【原创】哈工大材料力学上机——应力状态分析
Harbin Institute of Technology材料力学上机大作业课程名称：材料力学题目：应力状态分析院系：机电学院班级：姓名：学号：指导教师：设计时间： 20130621哈尔滨工业大学一、题目要求：应力状态分析输入：1.平面应力状态输入：（xy y x τσσ、、）；某截面方位角α2.空间应力状态输入：（zx yz xy z y x τττσσσ、、、、、）输出：1.输出主应力（321σσσ、、）。
2.最大剪应力（23113max σσττ-==）。
3.如为平面应力状态则需输出方位角α斜截面上的应力αατσ、及主方位角*σα。
4.画出应力圆示意图。
算法：1.二向应力状态给定xy y x τσσ、、，某截面方位角α，由公式ατασστατασσσσσ2cos 2sin 22sin 2cos 22'''xy y x y x xy y x yx x +--=+-++=即可求得方位角α斜截面上的应力αατσ、， 再由公式y x xyσστασ-=22tan ，推出2/)2arctan(y x xy σστασ-=，于是求出主方位角σα的两个值。
最后由公式22min max )2(2xy y x y x τσσσσσσ+-±+=⎭⎬⎫得到min max ,σσ，并将0,min max 与σσ比较大小，最后得到321,,σσσ 再由公式23113max σσττ-==计算出max τ程序：由matlAB 程序GUI 功能实现可视化窗口，使界面更加友好。
已知正交界面的正应力和剪应力，求其应力状态，包括主应力大小、最大剪切应力、以及具体界面位置应力和剪切力。
第一个按钮用来绘制莫尔圆；第二个按钮用来计算具体界面位置应力和剪切力。
function varargout = gui_moore_13_4(varargin)% GUI_MOORE_13_4 M-file for gui_moore_13_4.fig% GUI_MOORE_13_4, by itself, creates a new GUI_MOORE_13_4 or raises the existing% singleton*.%% H = GUI_MOORE_13_4 returns the handle to a new GUI_MOORE_13_4 or the handle to% the existing singleton*.%% GUI_MOORE_13_4('CALLBACK',hObject,eventData,handles,...) calls the local% function named CALLBACK in GUI_MOORE_13_4.M with the given input arguments.%% GUI_MOORE_13_4('Property','Value',...) creates a newGUI_MOORE_13_4 or raises the% existing singleton*. Starting from the left, property value pairs are% applied to the GUI before gui_moore_13_4_OpeningFunction gets called. An% unrecognized property name or invalid value makes property application% stop. All inputs are passed to gui_moore_13_4_OpeningFcn via varargin.%% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one% instance to run (singleton)".%% See also: GUIDE, GUIDATA, GUIHANDLES% Edit the above text to modify the response to help gui_moore_13_4% Last Modified by GUIDE v2.5 07-Sep-2007 18:01:33% Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @gui_moore_13_4_OpeningFcn, ...'gui_OutputFcn', @gui_moore_13_4_OutputFcn, ...'gui_LayoutFcn', [] , ...'gui_Callback', []);if nargin && ischar(varargin{1})gui_State.gui_Callback = str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); elsegui_mainfcn(gui_State, varargin{:});end% End initialization code - DO NOT EDIT% --- Executes just before gui_moore_13_4 is made visible.function gui_moore_13_4_OpeningFcn(hObject, eventdata, handles, varargin)% This function has no output args, see OutputFcn.% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% varargin command line arguments to gui_moore_13_4 (see VARARGIN) % Choose default command line output for gui_moore_13_4set(handles.sx_edit,'String', 30);set(handles.sy_edit,'String', 0);set(handles.txy_edit,'String', 8);set(handles.rad_edit,'String', 1.0472); %pi/3handles.output = hObject;handles.output = hObject;% Update handles structureguidata(hObject, handles);% UIWAIT makes gui_moore_13_4 wait for user response (see UIRESUME) % uiwait(handles.figure1);% --- Outputs from this function are returned to the command line. function varargout = gui_moore_13_4_OutputFcn(hObject, eventdata, handles)% varargout cell array for returning output args (see VARARGOUT);% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Get default command line output from handles structurevarargout{1} = handles.output;function txy_edit_Callback(hObject, eventdata, handles)% hObject handle to txy_edit (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of txy_edit as text% str2double(get(hObject,'String')) returns contents of txy_edit as a double% --- Executes during object creation, after setting all properties. function txy_edit_CreateFcn(hObject, eventdata, handles)% hObject handle to txy_edit (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunction sy_edit_Callback(hObject, eventdata, handles)% hObject handle to sy_edit (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of sy_edit as text% str2double(get(hObject,'String')) returns contents of sy_edit as a double% --- Executes during object creation, after setting all properties. function sy_edit_CreateFcn(hObject, eventdata, handles)% hObject handle to sy_edit (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunction sx_edit_Callback(hObject, eventdata, handles)% hObject handle to sx_edit (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of sx_edit as text% str2double(get(hObject,'String')) returns contents of sx_edit as a double% --- Executes during object creation, after setting all properties. function sx_edit_CreateFcn(hObject, eventdata, handles)% hObject handle to sx_edit (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end% --- Executes on button press in moore_pushbutton.function moore_pushbutton_Callback(hObject, eventdata, handles)% hObject handle to moore_pushbutton (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)。
哈工大材料力学大作业--matlab编程
H a r b i n I n s t i t u t e o f T e c h n o l o g y材料力学上机作业课程名称：材料力学设计题目：应力状态分析院系：机电学院班级：分析者：学号：指导教师：***设计时间：2013年6月18日哈尔滨工业大学材料力学上机课设计说明书一， 设计题目题目7 应力状态分析 输入：1. 平面应力状态输入：x y xy σστ（,,）；某截面方位角α2. 空间应力状态输入：,x y z xy yz zx σσστττ（,,,,）输出： 1. 输出主应力123σσσ（,,）2. 最大切应力（13max 132σσττ-==）3.如为平面应力状态则需要输出方位角α斜截面上的应力ααστ、及主方向角*σα4. 画出应力圆示意图二， 程序计算设计过程1. 平面应力状态分析对于任意平面应力状态，有max min σσ=2x y σσ+±主应力为：1max 23min ,0,σσσσσ===并且由 2tan 2xyx yστασσ=-可求得主应力方向角13σσαα、。
对于任意一个方位角α，有：=cos 2sin 222sin 2cos 22x yx yxy x yxy αασσσσσατασστατα++++-=-+从而，输入任意角α，即可求得该截面的应力状态ααστ、并且ααστ、都是关于α的函数，上式即为应力圆的参数方程，参数为α。
将α从0到pi 取一系列的值，则可以求出一系列的ααστ、，在坐标系中找到对应点，连接即可作出应力圆。
2. 三向应力状态分析解特征方程 321230I I I σσσ-+-=即可求出主应力123σσσ、、 其中：123||||||||x y z xyx y zy z xz xy y yz z zx x x yx zx xyy zy xzyz z I I I σσσστστσττστστσστττστττσ=-+⎛⎫⎛⎫⎛⎫=++⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎛⎫⎪= ⎪ ⎪⎝⎭再由 13max 132σσττ-== 可求得最大切应力。
材料力学上机大作业(matlab编)
一、可实现课题在如图所示的悬臂梁中，杆件为圆杆。
杆长为L，直径为D，材料弹性模量为E。
输入集中力F大小，作用点a，弯矩M，作用点b，即可求得悬臂梁的挠度曲线图。
二、程序代码clear alldisp('请给定材料信息'); %输入材料信息L=input('圆杆长度L(/M)=');D=input('圆杆直径D(/M)=');E=input('弹性模量E(/GPa)=');I=double(D^4*3.14/32);disp('请给定受力情况'); %输入受力情况F=input('切向集中力大小F(/N)=');a=input('切向集中力作用位置(/M)=');M=input('弯矩大小M(/N*M)=');b=input('弯矩作用位置(/M)=');x1=0:0.01:a; %F引入的挠度vx1=(-F*x1.^2*3*a+F*x1.^3)*(1/(6*E*10^9*I));x2=a:0.01:L;vx2=(-F*a.^2*3*x2+F*a.^3)*(1/(6*E*10^9*I));v11=[vx1,vx2];x11=[x1,x2];x3=0:0.01:b; %M引入的挠度vx3=(-M*x3.^2)*(1/(2*E*10^9*I));x4=b:0.01:L;vx4=(-M*b*x4+M*0.5*b.^2)*(1/(E*10^9*I));x22=[x3,x4];v22=[vx3,vx4];v33=v22+v11; %叠加plot(x11,v33),xlabel('x /M'),ylabel('v(x) /M')title('挠曲线图')grid on;三、使用方法运行代码输入圆杆长度（单位：m）输入圆杆直径（单位：m）输入弹性模量（单位：GPa）输入集中力大小（单位：N）（向下为正，若无请输入0）输入集中力作用位置（单位：m）（若无请输入0）输入弯矩大小（单位：N*m）（逆时针为正，若无请输入0）输入弯矩作用位置（单位：m）（若无请输入0）输出挠曲线图四、运行实例【实例1】圆杆同时受集中力与弯矩作用，输入、输出见下图。
哈工大材料力学计算机计算大作业_范例
B
ql 3 48 EI z
500 23 1.5215mrad 3 0.034 0.054 3 9 0.04 0.06 48 200 10 12 12 500 2 4 1.1411mm 3 0.034 0.0543 9 0.04 0.06 128 200 10 12 12
max
1000 32 1 5.77 mm ， x 3 1732.1mm 9 3 100000 3
电算解：
故：与理论解一致。
第 十四 页 共 二十一 页
材料力学大作业
（8）理论解：
增补条件： 代入约定条件算得： A M 3000 N m ， a 1 .6 m ， 电算解：
3
电算解：
故：与理论解一致。
第 十七 页 共 二十一 页
材料力学大作业
（11）理论解：
代入约定数据解之得： A B 11 .25mrad ； C
5 1000 34 10.55mm 384 100000
电算解：
故：与理论解一致。
第 十八 页 共 二十一 页
然后，我们需要输入抗弯刚度。题目没有给怎么办？不怕，我们有“抗弯刚度计算器”！点击对 应区段抗弯刚度边上的按钮即可开始计算！
第 三 页 共 二十一 页
材料力学大作业
在计算完毕之后，点击“置入并返回” ，则我们又来到了基本信息输入界面。由于这一根简单 梁是等截面、等抗弯刚度的，所以我们不用去管“区段二”和“区段三”的相关文本框（实际上在 你选定对应的区段之前，这些文本框都是无法激活的） 。我们需要点击“更新基本信息”按键，以 核查、写入数据。核查无误，我们点击“受力情况输入”按钮。
2016新编哈工大材料力学上机大作业四——绘制梁的剪力弯矩图
H a r b i n I n s t i t u t e o f T e c h n o l o g y材料力学上机作业课程名称：材料力学设计题目：绘制梁的剪力弯矩图院系：班级：分析者：学号：指导教师：设计时间：哈尔滨工业大学材料力学上机课设计说明书一，设计题目题目4 绘制梁的剪力弯矩图输入：1.梁的总长度l2.支撑条件及量的各区段长度输入（左、右固定端悬臂梁；简支梁；左、右、双外伸梁）3.各载荷大小、作用位置及方向（q i、a i、b i；p j、c j、m k、d k）输出：1.结构构型图（图示）2.剪力、弯矩（图示）3.输出剪力、弯矩的最大值及截面位置。
二，方向规定本程序规定集中作用力及均布载荷以向下为正，右固定端悬臂梁与左外伸梁集中力偶以顺时针方向为正，其他情况集中力偶以逆时针方向为正。
当取出梁的一段为研究对象时，梁左端面力以向下为正，力偶以顺时针方向为正，梁右端面力以向上为正，力偶以逆时针方向为正。
三，程序设计过程1，制作程序框架，显示提示内容，提示操作者需要输入的作用条件及各作用位置；2，编写程序使计算机读入操作者输入的作用条件；3，草稿拟写各种情况下为达到题目要求所需使用的计算公式设量的长度为l,集中力大小为p，作用位置为c，集中力大小为q，作用起始位置a，终止位置为b，集中力偶大小为m，作用位置d。
（1）左固定端悬臂梁：在任意位置x处，取x以右部分为研究对象a，终止位置为b，集中力作用位置为c，大小为P，集中力偶作用位置为d，大小为m。
①若c<x,a<b<x,d>x,则Fs y=0,M(x)=m；②若c>x,a<b<x,d<x,则Fs y=-P,M(x)=Px-Pc；③若c<x,a<x<b,d<x,则Fs y=-q(b-x),M(x)=-2q(b-x)²；④若c<x, x<a<b,d>x,则Fs y=-q(b-a),M(x)=-q(b-a)(2ba+-x)；⑤若c>x,a<b<x,d>x,则Fs y=-P,M(x)= Px-Pc+m；(第①、②两种情况合成)⑥若c<x,a<x<b,d>x,则Fs y=-q(b-x),M(x)=m-2q(b-x)²；(第①、③两种情况合成)⑦若c<x,x<a<b,d>x,则Fs y=-q(b-a),M(x)=m-q(b-a)(2ba+-x)；(第①、④两种情况合成)⑧若c>x,a<x<b,d<x,则Fs y=-P-q(b-x),M(x)=Px-Pc-2q(b-x)²；(第②、③两种情况合成)⑨若c>x,x<a<b,d<x,则Fs y=-P-q(b-a),M(x)=Px-Pc-q(b-a)*(2ba+-x)；(第②、④两种情况合成)⑩若c>x,x<a<b,d>x,则Fs y=-P-q(b-a), M(x)=m+Px-Pc-q(b-a)*(2ba+-x)；(第①、②、④两种情况合成)⑪c>x,a<x<b,d>x, 则Fs y=-P-q(b-x), M(x)m+Px-Pc-q(b-a)*(2ba+-x); (第①、②、③两种情况合成)将上述公式编入程序即可计算出在固定端悬臂梁情况下任意位置处的剪力和弯矩，采用散点法作出梁的剪力弯矩图。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本程序只支持静定结构的梁（左端悬臂梁、简支梁）函数输入格式：beamsolver(L,EI,supports,loads,maxdx);参量的输入格式：L=10.0EI=2e8supports={{'f',0},{'v',10.0}} ----左端悬臂梁supports={{'p',2.0},{'r',8.0}} ----简支梁loads={{'f',[2.0,1000]},{'m',[4.0,500]},{'d',[7.0,9.0,3.0,20,100]}}maxdx=0.01输出为V,M,vy，x的一维行向量和三张坐标图----图1:剪力V图,图2:弯矩M图,图3:挠度vy图.范例：悬臂梁：纯受集中力：beamsolver(10.0,2e8,{{'f',0},{'v',10.0}},{{'f',[2.0,1000]}},0.01)纯受集中矩：beamsolver(10.0,2e8,{{'f',0},{'v',10.0}},{{'m',[4.0,500]}},0.01)纯受分布力（格式一）：beamsolver(10.0,2e8,{{'f',0},{'v',10.0}},{{'d',[7.0,9.0,3.0,20,100]}},0.01) 纯受分布力（格式二）：beamsolver(10.0,2e8,{{'f',0},{'v',10.0}},{{'d',[(7.0:0.01:9.0);(0:0.01:2.0).^2]} },0.01)受混合力：beamsolver(10.0,2e8,{{'f',0},{'v',10.0}},{{'f',[2.0,1000]},{'m',[4.0,500]},{'d', [7.0,9.0,3.0,20,100]}},0.01)简支梁：纯受集中力：beamsolver(10.0,2e8,{{'p',2.0},{'r',8.0}},{{'f',[5.0,1000]}},0.01)纯受集中矩：beamsolver(10.0,2e8,{{'p',2.0},{'r',8.0}},{{'m',[4.0,500]}},0.01)纯受分布力（格式一）：beamsolver(10.0,2e8,{{'p',2.0},{'r',8.0}},{{'d',[7.0,9.0,3.0,20,100]}},0.01) 纯受分布力（格式二）：beamsolver(10.0,2e8,{{'p',2.0},{'r',8.0}},{{'d',[(7.0:0.01:9.0);(0:0.01:2.0).^2] }},0.01)受混合力：beamsolver(10.0,2e8,{{'p',2.0},{'r',8.0}},{{'f',[5.0,1000]},{'m',[4.0,500]},{'d' ,[7.0,9.0,3.0,20,100]}},0.01)以上范例输出的剪力V图、弯矩M图都经过笔算检验完全正确。
悬臂梁的挠度经检验为正确的，简支梁的挠度有误（整个程序的唯一错误）。
function [x,V,M,vy]=beamsolver(L,EI,supports,loads,maxdx)%材料力学上机课作业（满分）%本DOC文件由M文件转换而成，可能有格式错误%班级;1036***(专业班级1036***)%支持两种梁的模型:%1.左端悬臂梁%2.一端pin一端roller的简支梁%输入格式：beamsolver(L,EI,supports,loads,maxdx);%例如：L=10.0,EI=2e8,supports={{'f',0},{'v',10.0}},loads={{'f',[2.0,1000]},{'m',[4.0,50 0]},{'d',[7.0,9.0,3.0,20,100]}},maxdx=0.01%输出为三张坐标图,图1:剪力V图,图2:弯矩M图,图3:挠度vy图.%fixed-freeif (supports{1}{1}=='f' && supports{2}{1}=='v')n=1;M=linspace(0,0,L/maxdx);V=linspace(0,0,L/maxdx);while n<=length(loads)%P------------------------fixed-free--------------------------------------if loads{n}{1}=='f'x=linspace(maxdx,L,L/maxdx); %将x分为小量n1=loads{n}{2}(1)/maxdx; %确定x=L1处对应的下标V1=linspace(loads{n}{2}(2),loads{n}{2}(2),n1); %第一段剪力赋值V2=zeros(1,L/maxdx-n1); %第二段剪力赋值M1=-loads{n}{2}(2)*(loads{n}{2}(1)-x(1:n1)); %第一段弯矩赋值M2=zeros(1,L/maxdx-n1); %第二段弯矩赋值M=M+[M1,M2]; %全梁的弯矩V=V+[V1,V2]; %全梁的剪力end%M-----------------------fixed-free----------------------------------------- if loads{n}{1}=='m'x=linspace(maxdx,L,L/maxdx); %将x分为小量n1=loads{n}{2}(1)/maxdx; %确定x=L1处对应的下标M1=linspace(loads{n}{2}(2),loads{n}{2}(2),n1); %第一段弯矩赋值M2=zeros(1,L/maxdx-n1); %第二段弯矩赋值M=M+[M1,M2]; %全梁的弯矩end%q-----------------------fixed-free---------------------------------------if loads{n}{1}=='d'x=linspace(maxdx,L,L/maxdx); %将x分为小量n1=loads{n}{2}(1)/maxdx; %确定x=L1处对应的下标n2=loads{n}{2}(2)/maxdx; %确定x=L2处对应的下标k=length(loads{n}{2})-4;A=size(loads{n}{2});if A(1)==1 %如果是第一种q的输入情况ks=1;q2=linspace(loads{n}{2}(4),loads{n}{2}(4),n2-n1+1); while ks<=kk0=loads{n}{2}(4+ks);q2=q2+k0*(x(n1:n2)-loads{n}{2}(3));ks=ks+1;endq1=zeros(1,n1-1);q3=zeros(1,L/maxdx-n2);q0=[q1,q2,q3];q=fliplr(q0); %求出qVq0=cumsum(q)*maxdx;Vq=fliplr(Vq0); %求出剪力VMq0=-cumsum(Vq0)*maxdx;Mq=fliplr(Mq0); %求出弯矩Mendif A(1)==2 %如果是q的第二种输入情况q2=loads{n}{2}(2,:); %表示出q2n1=loads{n}{2}(1,1)/maxdx; %定出q开始的坐标B=length(loads{n}{2}); %求q的长度n2=loads{n}{2}(1,B)/maxdx; %定出q结束的坐标q1=zeros(1,n1-1);q3=zeros(1,L/maxdx-n2);q0=[q1,q2,q3]; %求出q0q=fliplr(q0);Vq0=cumsum(q)*maxdx;Vq=fliplr(Vq0); %求出剪力VMq0=-cumsum(Vq0)*maxdx;Mq=fliplr(Mq0); %求出弯矩MendM=M+Mq; %全梁的弯矩V=V+Vq; %全梁的剪力endn=n+1;endV=VM=MA=cumsum(M)*maxdx/EI; %对弯矩积分求转角vy=cumsum(A)*maxdx %对转角积分求挠度subplot(3,1,1),plot(x,V),grid; %绘图subplot(3,1,2),plot(x,M),grid;subplot(3,1,3),plot(x,vy),grid;end%分割线---------------------------------------------------------------------------------%pin-roller-----------------------------------------------------------------------------if (supports{1}{1}=='p' || supports{1}{1}=='r')&&(supports{2}{1}=='p' || supports{2}{1}=='r')n=1;np=supports{1}{2}/maxdx;nr=supports{2}{2}/maxdx;M=linspace(0,0,L/maxdx);V=linspace(0,0,L/maxdx);while n<=length(loads)%P------------------------pin-roller--------------------------------------if loads{n}{1}=='f'x=linspace(maxdx,L,L/maxdx); %将x分为小量n1=loads{n}{2}(1)/maxdx; %确定x=L1处对应的下标Nr=loads{n}{2}(2)*(loads{n}{2}(1)-supports{1}{2})/(supports{2}{2}-supports{1}{2} );if n1<=np %力作用在在支架左边V1=zeros(1,n1-1);V2=-linspace(loads{n}{2}(2),loads{n}{2}(2),np-n1);V3=-linspace(Nr,Nr,nr-np);V4=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3,V4];M1=zeros(1,n1-1);M2=-loads{n}{2}(2)*(x(n1:np-1)-loads{n}{2}(1));M30=Nr*x(1:nr-np);M3=fliplr(M30);M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4];endif (n1>np && n1<=nr) %力作用在两个支架之间V1=zeros(1,np-1);V2=linspace(loads{n}{2}(2)-Nr,loads{n}{2}(2)-Nr,n1-np);V3=-linspace(Nr,Nr,nr-n1);V4=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3,V4];M1=zeros(1,np-1);M2=(loads{n}{2}(2)-Nr)*(x(np:n1-1)-supports{1}{2});M3=Nr*(supports{2}{2}-x(n1:nr-1));M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4];endif n1>nr %力作用在支架右边V1=zeros(1,np-1);V2=linspace(loads{n}{2}(2)-Nr,loads{n}{2}(2)-Nr,nr-np);V3=linspace(loads{n}{2}(2),loads{n}{2}(2),n1-nr);V4=zeros(1,L/maxdx-n1+1);V=V+[V1,V2,V3,V4];M1=zeros(1,np-1);M2=(loads{n}{2}(2)-Nr)*(x(np:nr-1)-supports{1}{2});M3=-loads{n}{2}(2)*(loads{n}{2}(1)-x(nr:n1-1));M4=zeros(1,L/maxdx-n1+1);M=M+[M1,M2,M3,M4];endend%M-----------------------pin-roller----------------------------------------- if loads{n}{1}=='m'x=linspace(maxdx,L,L/maxdx); %将x分为小量n1=loads{n}{2}(1)/maxdx; %确定x=L1处对应的下标if n1<=np %弯矩作用在支架左边Np=loads{n}{2}(2)/(supports{2}{2}-supports{1}{2});Nr=Np;V1=zeros(1,np-1);V2=linspace(Np,Np,nr-np);V3=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3];M1=zeros(1,n1-1);M2=-linspace(loads{n}{2}(2),loads{n}{2}(2),np-n1);M3=-Nr*(supports{2}{2}-x(np:nr-1));M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4]; %全梁的弯矩endif (n1>np && n1<=nr) %弯矩作用在两支架中间Np=loads{n}{2}(2)/(supports{2}{2}-supports{1}{2});Nr=Np;V1=zeros(1,np-1);V2=linspace(Np,Np,nr-np);V3=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3];M1=zeros(1,np-1);M2=Np*(x(np:n1-1)-supports{1}{2});M3=-Nr*(supports{2}{2}-x(n1:nr-1));M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4]; %全梁的弯矩endif n1>nr %弯矩作用在支架右边Np=loads{n}{2}(2)/(supports{2}{2}-supports{1}{2});Nr=NpV1=zeros(1,np-1);V2=linspace(Np,Np,nr-np);V3=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3];M1=zeros(1,np-1);M2=Np*(x(np:nr-1)-supports{1}{2});M3=linspace(loads{n}{2}(2),loads{n}{2}(2),n1-nr);M4=zeros(1,L/maxdx-n1+1);M=M+[M1,M2,M3,M4]; %全梁的弯矩endend%q-----------------------pin-roller--------------------------------------- if loads{n}{1}=='d'x=linspace(maxdx,L,L/maxdx); %将x分为小量n1=loads{n}{2}(1)/maxdx; %确定x=L1处对应的下标n2=loads{n}{2}(2)/maxdx; %确定x=L2处对应的下标k=length(loads{n}{2})-4;A=size(loads{n}{2});if A(1)==1 %如果是第一种q的输入情况ks=1;q2=linspace(loads{n}{2}(4),loads{n}{2}(4),n2-n1+1);while ks<=kk0=loads{n}{2}(4+ks);q2=q2+k0*(x(n1:n2)-loads{n}{2}(3));ks=ks+1;endq1=zeros(1,n1-1);q3=zeros(1,L/maxdx-n2);q0=[q1,q2,q3];F=q0*maxdx; %等效为很多个力作用在杆上nF=length(F); %找出F的最后一点值得坐标X=1;while X<=nFif F(X)~=0Nr=F(X)*(X-np)/(nr-np);if X<=np %力作用在在支架左边V1=zeros(1,X-1);V2=-linspace(F(X),F(X),np-X);V3=-linspace(Nr,Nr,nr-np);V4=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3,V4];M1=zeros(1,X-1);M2=-F(X)*(x(X:np-1)-(X*maxdx));M30=Nr*x(1:nr-np);M3=fliplr(M30);M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4];endif (X>np && X<=nr) %力作用在两个支架之间V1=zeros(1,np-1);V2=linspace(F(X)-Nr,F(X)-Nr,X-np);V3=-linspace(Nr,Nr,nr-X);V4=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3,V4];M1=zeros(1,np-1);M2=(F(X)-Nr)*(x(np:X-1)-supports{1}{2}); M3=Nr*(supports{2}{2}-x(X:nr-1));M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4];endif X>nr %力作用在支架右边V1=zeros(1,np-1);V2=linspace(F(X)-Nr,F(X)-Nr,nr-np);V3=linspace(F(X),F(X),X-nr);V4=zeros(1,L/maxdx-X+1);V=V+[V1,V2,V3,V4];M1=zeros(1,np-1);M2=(F(X)-Nr)*(x(np:nr-1)-supports{1}{2}); M3=-F(X)*(X*maxdx-x(nr:X-1));M4=zeros(1,L/maxdx-X+1);M=M+[M1,M2,M3,M4];endendX=X+1;endendif A(1)==2 %如果是q的第二种输入情况q2=loads{n}{2}(2,:); %表示出q2n1=loads{n}{2}(1,1)/maxdx; %定出q开始的坐标B=length(loads{n}{2}); %求q的长度n2=loads{n}{2}(1,B)/maxdx; %定出q结束的坐标q1=zeros(1,n1-1);q3=zeros(1,L/maxdx-n2);q0=[q1,q2,q3]; %求出qF=q0*maxdx; %等效为很多个力作用在杆上nF=length(F); %找出F的最后一点值得坐标X=1;while X<=nFif F(X)~=0Nr=F(X)*(X-np)/(nr-np);if X<=np %力作用在在支架左边V1=zeros(1,X-1);V2=-linspace(F(X),F(X),np-X);V3=-linspace(Nr,Nr,nr-np);V4=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3,V4];M1=zeros(1,X-1);M2=-F(X)*(x(X:np-1)-(X*maxdx));M30=Nr*x(1:nr-np);M3=fliplr(M30);M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4];endif (X>np && X<=nr) %力作用在两个支架之间V1=zeros(1,np-1);V2=linspace(F(X)-Nr,F(X)-Nr,X-np);V3=-linspace(Nr,Nr,nr-X);V4=zeros(1,L/maxdx-nr+1);V=V+[V1,V2,V3,V4];M1=zeros(1,np-1);M2=(F(X)-Nr)*(x(np:X-1)-supports{1}{2}); M3=Nr*(supports{2}{2}-x(X:nr-1));M4=zeros(1,L/maxdx-nr+1);M=M+[M1,M2,M3,M4];endif X>nr %力作用在支架右边V1=zeros(1,np-1);V2=linspace(F(X)-Nr,F(X)-Nr,nr-np);V3=linspace(F(X),F(X),X-nr);V4=zeros(1,L/maxdx-X+1);V=V+[V1,V2,V3,V4];M1=zeros(1,np-1);M2=(F(X)-Nr)*(x(np:nr-1)-supports{1}{2}); M3=-F(X)*(X*maxdx-x(nr:X-1));M4=zeros(1,L/maxdx-X+1);M=M+[M1,M2,M3,M4];endendX=X+1;endendendn=n+1;endV=VM=MA=cumsum(M)*maxdx/EI; %对弯矩积分求转角 vy=cumsum(A)*maxdx %对转角积分求挠度subplot(3,1,1),plot(x,V),grid; %绘图 subplot(3,1,2),plot(x,M),grid;subplot(3,1,3),plot(x,vy),grid;end。