论文-基于MATLAB的汽车灯光优化分析.docx

北华航天工业学院毕业论文毕业设计报告(论文)报告(论文)题目： 基于MATLAB的汽车灯光 优化分析 作者所在系部： 机电工程学院 作者所在专业： 车辆工程 作者所在班级： B13142 作 者 姓 名 ： 宋应红 作 者 学 号 ： 2010322206 指导教师姓名： 白亚双 完 成 时 间 ： 2017年6月 北华航天工业学院教务处制摘 要汽车是现代人类生活中的重要工具，汽车在生活中的使用变得频繁。那么，它的安全性是人们非常重视的。汽车的灯光系统是汽车的重要组成部分，它保证行车安全，提高行驶速度，还可以在汽车转弯、超车、会车和停车时作为表达驾驶员意图的工具。汽车车灯光源是核心部分【1】关文达.汽车构造.第三版.北京：机械工业出版社，2010【2】 舒平、姚国平.汽车电器设备与维修.北京：北京理工大学出版社，2009【3】赵福堂.汽车电器宇电子设备.北京：北京理工大学出版社，2009。本文基于MATLAB软件对车灯光源功率和照亮区域进行优化分析。关键词：灯光优化 MATLAB 光源功率IAbstractAutomotive are an important tool in modern human life, and automotive are used more frequently in life. So, its safety is very important.Automotive lighting system is an important part of the car, it guarantee the driving safety, improve the speed, also can turn in the car, parking is passing, passing, and used as a tool of expression of driver intention. Automotive lamps are the core. Based on MATLAB software, this paper analyzes the power and illumination area of lamp light source.Key words:The lights to optimize MATLAB lamp powerII目 录摘 要I目 录III第1章 绪论11.1课题背景11.2 MATLAB在汽车领域的应用11.3 MATLAB在汽车灯光中的优点和发展前景1第2章 问题分析32.1 模型假设32.2 模型符号说明3第3章 点光源直射光形成的锥面方程确定43.1 旋转抛物面的方程确定43.2 点光源直射光锥面方程的确定43.3 光源直射功率有效值计算63.4直射光总功率计算11第4章 测试屏上直射光的亮区模型建立134.1 本章导读134.2亮区模型建立134.3 区域D1中任意一点处的光照强度计算154.4区域D21中任意一点出的光照强度确定164.5 区域D22中任意的一点的光照强度确定184.6 区域D23、D24中任意的一点的光照强度确定194.7 直射光区域内任意一点出的光照强度确定21总结27致 谢28参考文献29IV第1章 绪论1.1课题背景汽车头部的车灯的形状是一个旋转抛物面，汽车灯的水平正前方与旋转抛物面的对称轴重合，抛物面的开口半径为38mm，深度22.5mm。经过车灯的焦点，在与对称轴相垂直的水平方向，对称地放置长度为4mm的线光源，线光源均匀分布。在焦点F正前方25米处的A点放置一测试屛，屏与FA垂直。【4】韩瑞萍，董翔宇，郑毅，谭欣欣.大连：车灯线光源的优化设计,大连大学学报，20041.2 MATLAB在汽车领域的应用汽车产品在生产和制造的过程中，任何一道工序都需要离精细的加工和打磨，每一个微小的失误都有可能导致结构出现错误而无法验收，甚至造成重新返工，这样不仅延长生产的周期，还带来了经济的损失。MATLAB就可以在造出实车之前，通过虚拟仿真技术，看到遇到的问题，并找到合适的解决方案。1.3 MATLAB在汽车灯光中的优点和发展前景优点：（1）提高产品开发质量。在产品的开发研究中，通过虚拟样机来寻找最优设计，追求全性能最佳。（2）减短产品的开发周期。传统的产品开发，要经历方案设计图样设计样机制造工艺设计（寻找最合适的样机）测试性能投产，出现问题就需要往复循环，而且样机的制造，工艺测试都要花费很久的时间。而应用MATLAB的仿真技术，可以同时对多个方案进行性能模拟预测，找到最佳方案。在图样阶段，对结构参数进行验证，及早解决一些工艺设计，加工制造中的一些问题。（3）降低开发费用。用仿真代替样机或实体模型试验，可以节省大量的开发费用。发展前景:MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案,解决了实际应用的费时，【5】张志涌，杨祖樱.MATLAB教程.北京：北京航空航天学院出版社，2016经济的节约问题。在科研方面、在地产方面、在医学方面等等。MATLAB虚拟仿真的产品更安全、更直观、更有吸引力。越来越多的领域将引用仿真技术。35第2章 问题分析由于光源与反射面的距离与光源宽度相比较小，因此不能把光源作为一个点光源处理。为了方便处理，将灯丝当作长度一定、发光均匀的线光源，就可以把灯丝看作由若干个微小段组成，微小段相对于反射面的距离很小，可以看作点光源，用点光源积分处理。灯光光源发出的光线由旋转抛物面反射后的光线一般互相不平行，所以需要按照光的反射定律，求出反射光线。功率直射的有效值：点光源照射到物体上，以点光源为圆心作出一个球面，则光源射到物体上的光线在球面上的投影曲面面积与球面面积之比。是功率直射的有效值。功率比值：光源射到两物体上的光功率有效值之比，称为光源射到两物体上的光功率之比，即两物体在小球面上投影的曲面面积之比【6】石祥顺.物理光学与应用光学第二版.西安：西安电子科技大学出版社，2010。2.1 模型假设1.线光源长度分布是均匀的；2.将线光源视作无数个亮度相等的点光源的组合；3.只考虑一次反射的情况；4.线光源的功率额定值设定为1个单位；5.点光源各方向的光强相同；6.反射点前后功率损失率为常数。2.2 模型符号说明L：线光源长度（单位：mm），文中取4mm；P0：线光源的功率值，计算中取1；P1,P2:线光源直射光线和反射光线功率；R：直射光线总功率与反射光线总功率的比值；P（0，w，16）：线光源上的一点，其中;C:测试屏上的点；第3章 点光源直射光形成的锥面方程确定3.1 旋转抛物面的方程确定图3-1 旋转抛物面如图3-2所示，M、N为车灯边沿所在平面与抛物面、平面yOz的交点，因为MON是抛物面的母线，设母线方程为： （3-1）【7】同济大学数学系.高等数学.第六版.北京：高等教育出版社，2007：77根据M的位置关系，把点M的坐标为（0，38，22.5）带入 可得p=，并可求得焦点F的坐标为（0，0，16）。由高等数学知识可求得，母线绕Z轴旋转所形成的旋转抛物面方程为： （3-2）由于线光源经过焦点并且与y轴平行，所以有线光源上任意一点P的坐标是（0，w，16），并且。3.2 点光源直射光锥面方程的确定从（0，w,16）点的点光源发出的光线经过抛物面开口能形成一个锥面，该锥面的顶点坐标P，准线是旋转抛物面开口形成的圆，即： （3.3）点光源P直射光形成的锥面如图图3-2 点光源直射光锥面在锥面上的任意一点A（x，y，z），母线PA的参数方程为： （3.4）设PA和准线的交点是，则与其对应的参数为，所以； （3.5）还满足方程： （3.6）由上面解得：把 带入准线方程，就可求得准线方程是：， （3.7）3.3 光源直射功率有效值计算点光源P射到测试屏上的直射光功率和点光源的光线总功率的比值称为功率直射有效值。记成。计算方法:以球心是点光源P做一个半径为的小球面，直射到测试屏的光线在小球面上的投影面积为，则小球面的面积是,所以有： （3.8）为了方便计算，小球面的半径取1。锥面方程的标准化： （3.9)平移变换为： (3.10)在上面平移变换之下，变成： (3.11)令： （3.12）这是一个二次型方程，用正交变换的方法化成标准型。【8】牛莉.线性代数.第二版.北京：中国水利水电出版社，2009编写程序求A的三个特征值：输入程序：clear all;close all; % Initialsyms wA=1 0 0;0 1 w/6.5;0 w/6.5 (w2-382)/(6.52);D,V=eig(A)运行输出结果：D=1,0,00,-(2*(w2-1444)/(13*w)-(13*(4*w2-304*w+5945)*(4*w2+304*w+5945)(1/2)/338-(2*w2)/169+5607/338)/(2*w),(13*(4*w2-304*w+5945)*(4*w2+304*w+5945)(1/2)/338+(2*w2)/169-5607/338)/(2*w)-(2*(w2-1444)/(13*w)0,1,1V=1,0,00,(2*w2)/169-(4*w2-304*w+5945)*(4*w2+304*w+5945)(1/2)/338-5607/338,00,0,(4*w2-304*w+5945)*(4*w2+304*w+5945)(1/2)/338+(2*w2)/169-5607/338整理结果如下： 该准面在正交变换的标准型可设为： （3.13）锥面方程为： （3.14）令：则锥面方程可以表示为： (3.15) 接下来确定圆锥面截下单位球面部分的面积。 设单位小球面的方程为：，所求曲面在xOy 平面内的投影区域是： （3.16） 于是： （3.17）因为含有的形式，可以用极坐标代换:积分区域被变换为： （3.18）所以： （3.19）其中，根据题目的条件和问题的特殊性，简化问题计算。回到二次型矩阵的特征值计算结果，分析得出：考虑的极值分布，设：（3.20）可见是偶函数，经过配方为： (3.21)求导的： (3.22)近似的处理特征值程序如下： clear all; close all; format short % Initial syms w fw=(4*w2+(16*w4-44586*w2+35343025)(1/2)-5607)/388-1; dfw=diff(fw,w)dfw=(2*w)/97-(-64*w3+89172*w)/(776*(16*w4-44586*w2+35343025)(1/2) solve(dfw)ans =0 subs(fw,w,2)ans =35164937(1/2)/388 - 5979/388 er=fw/(4*w2+(16*w4-44586*w2+35343025)(1/2)-5607)/338); dferw=diff(er,w); solve(dferw)ans=0 subs(er,w,-2)ans=(35164937(1/2)/388 - 5979/388)/(35164937(1/2)/338 - 5591/338)整理程序： clear all; close all; format short % Initial syms w fw=(4*w2+(16*w4-44586*w2+35343025)(1/2)-5607)/388-1; dfw=diff(fw,w)solve(dfw)subs(fw,w,2) er=fw/(4*w2+(16*w4-44586*w2+35343025)(1/2)-5607)/338); dferw=diff(er,w); solve(dferw) subs(er,w,-2)整理输出结果：dfw=(2*w)/97-(-64*w3+89172*w)/(776*(16*w4-44586*w2+35343025)(1/2)Ans=0ans=35164937(1/2)/388-5979/388ans=0ans=(35164937(1/2)/388-5979/388)/(35164937(1/2)/338-5591/338) 此解小于0唯一驻点w=0,的绝对值不超过0.0030，所以取=1，此时的相对误差：经过计算，得到，在问题中，取，从而，计算时，。经过近似后，可以计算，程序是： format short syms a t r1 r r1=1/(1+1/(a2)(1/2); Sw=int(int(r/sqrt(1-r2),r,0,r1),t,0,2*pi) 运行程序结果：Sw=-2*pi*(1/(a2+1)(1/2)-1)整理结果得： （3.23）因此，在点光源P的直射有效值功率的问题中计算中，可近似的取： （3.24）所以，点光源P的功率是P（w）,到达测试屏的直射光线部分的功率是： （3.25）反射光线部分的功率为： （3.26）3.4直射光总功率计算建立直射光的功率模型图：图3-3 直射光总功率如图所示，在灯丝上取一点P（0，w，16）,并沿灯丝微元dw点（0，w+dw,16）,则的长度是dw,将发光灯丝当作（0，w,16）处的点光源，它的发光功率为： （3.27）由点光源直射光功率公式，在点P（0，w，16）的直射光线功率微元为： （3.28）所以，全段灯丝的直射光线总功率为： （3.29）其中有： （3.30）反射光线的总功率为： (3.31)直射光总功率与反射光总功率的比值为： (3.32)求解全段灯丝的直射光总功率公式定积分，由于L=4，都是常数，所以： (3.33)对上面积分，利用MATLAB数值积分函数quad 计算，程序如下：clc,clear,close all format shortfw=(w)(388)(1/2)./(-(4*w.2-5995)+(4*w.2+5607).2+19762).(1/2).(1/2);I1=quad(fw,0,2);P1=1/2-1/4*I1P2=1-P1R=P1/P2结果如下：P1=0.4099P2=0.5901R=0.6945整理结果为：P1=0.4099，P2=0.5901,R=0.6945。由此可得出，直射光线总功率与反射光线总功率的比值是0.6945。第4章 测试屏上直射光的亮区模型建立4.1 本章导读本章需要求解灯光亮区问题，所以需要了解线光源功率的度量问题。根据所学光学知识，光源的功率P与发光强度及线光源产生光通量成正相关，所以可以用发光强度I来衡量光源的功率。4.2亮区模型建立分析测试屏上直射光线的区域，P（0,w,16）处点光源发出的直射光线形成的锥面方程是： (4.1)该锥面与测试屏z=2016的交线是：将z值代入上式，得到方程：该交线在测试屏上是以点（）为圆心，为半径的圆。因此，P（0，w，16）处点光源发出的直射光在测试屏上的区域是一个圆形。随着点P（0，w，16）在直线上，从y=-2到y=2运动，测试屏上亮区的圆心在直线上，从到运动。如图：图4-1 测试屏上直射光区域该区域边界由直线段AB，圆弧BNC，直线段CD，圆弧DMA组成，各线段的方程分别为：直线段AB：直线的CD：圆弧BNC：圆弧DMA：从灯丝发出的直射光线在测试屏上的亮区，并不是灯丝任意一点发出的直射光都能到达亮区的任意一点，亮区中的点分为两个不一样的部分。如图：图4-2直射光区域图中D1：灯丝发出的直射光线能照射到区域中的任意一点，该区域是两个圆形区域的相交部分，两圆在z=25016的平面上，分别有方程表示： (4.2) (4.3)则边界曲线交点是：4.3 区域D1中任意一点处的光照强度计算设亮区中的一点T在D1中，则点T与点P(0,w,16)的距离为：直线PT与平面z=25016之间夹角的余弦值为：因此，点光源P（0，w,16）直射到点T的光照强度为： (4.4)式中，P（w）是点光源的功率。利用微元法，可求出线光源上的所有点在点T的光照强度是： (4.5)直射光线区域D1中的点光照度的积分求解程序：clc,clear,close all;format shortssyms x y wI1=int(1/(x2+(y-w)2+250002)(3/2),w)Iup=subs(I1,w,2);Ilow=subs(I1,w,-2);Iv=factor(simplify(Iup-Ilow)结果如下：I1=(w-y)/(x2+625000000)*(w-y)2+x2+625000000)(1/2)Iv=y*(y-2)2+x2+625000000)(1/2)-y*(y+2)2+x2+625000000)(1/2)+2*(y-2)2+x2+625000000)(1/2)+2*(y+2)2+x2+625000000)(1/2),1/(y-2)2+x2+625000000)(1/2),1/(y+2)2+x2+625000000)(1/2),1/(x2+625000000) 整理结果：(4.6)4.4区域D21中任意一点出的光照强度确定区域D21是平面z=25016中包含在一个圆形区域内部：但在下列圆形区域外：所以 对于区域D21有：与区域D1相似，点光源P（0，w，16）直射到点T的光照强度是为： (4.7)其中，P（w）是点光源功率。 确定线光源能够直射到D21中的光线有哪些；当w 从-2到2时，由点P（0，w，16）出点光源发出的直射光从直射光区域图最上面的圆形区域平移到最下面的圆形区域。对于D21中点T,当w=-2时，有点P（0，w，16）处点光源发出的直射光线能到达该点；当w从-2变化到某个w0值时，由点P（0，w，16）处点光源形成的直射光线亮区的边界圆周正好通过点T；当ww0时，点T就在由P（0，w，16）处点光源形成的直射光亮区的边界圆周之外，因此，用微元法，D21中的点T的光照强度为： (4.8)w0的确定方法，找出过点T的边界圆周的圆心。在z=25016的平面内，以点为圆心，作半径为的圆，确定圆与直线x=0的交点y1,解方程组：利用P（0，w，16）点光源与测试屏上形成的亮区圆心关系：从而求出,求得：直射光线区域D21中的点的光照度积分求解程序：clc,clear,close allformat shortsyms x y wI1=int(1/(x2+(y-w)2+250002)(3/2),w);w0=6.5/(22.5-25016)*(y+(38/6.5*25000)2-x2)(1/2);Iup=subs(I1,w,2);Ilow=subs(I1,w,w0);Ev=factor(simplify(Iup-Ilow)Iv1=subs(subs(char(Ev),x,X0(ID(i),JD(i),y,Y0(ID(i),JD(i)4.5 区域D22中任意的一点的光照强度确定区域D22是平面z=25016中包含在一个圆形区域内部：但在下列圆形区域外：所以 对于区域D22有：与区域D22相似，通过D21中的光照强度的相似方法求出点光源P（0，w，1直射到点T的光照强度是为： (4.9)其中：直射区域D22点的光照强度积分求解程序如下：clc,clear,close allformat shortsyms x y wI1=int(1/(x2+(y-w)2+250002)(3/2),w);w0=6.5/(22.5-25016)*(y-(38/6.5*25000)2-x2)(1/2);Iup=subs(I1,w,2);Ilow=subs(I1,w,w0);Ev=factor(simplify(Iup-Ilow)Iv1=subs(subs(char(Ev),x,X0(ID(i),JD(i),y,Y0(ID(i),JD(i)4.6 区域D23、D24中任意的一点的光照强度确定区域D23，D24为平面z=25016内的圆：分别与直线段AB：直线段CD：围成，因此对于区域D23和D24的描述为：对于区域D23、D24中的任意一点T,当w=-2时，由P（0，w，25016）处点光源发出的直射光线不能到达该点，当w从-2变化到时,点光源形成的直射光区域的边界正好到达T；当w从开始继续变化时，点TT正好在点光源形成的直射光亮区内；而当w达到某个（大于等于）时，点光源形成的直射光区域的边界正好到达T，当ww1时，点T就在点光源形成的直射光亮区的边界之外，因此，利用微元法，D23、D24中点T的光照强度为： (4.10)其中，的确定方法与D21中的方法类似，分别与在z=25016平面内以点T为圆心，半径为的圆与直线x=0的两个相交点对应。他们的值为：直射区域D22点的光照强度积分求解程序如下：clc,clear,close allformat shortsyms x y wI1=int(1/(x2+(y-w)2+250002)(3/2),w);w0=6.5/(22.5-25016)*(y+(38/6.5*25000)2-x2)(1/2);w1=6.6/(22.56-25016)*(y-(38/6.5*25000)2-x2)(1/2);Iup=subs(I1,w,w0);Ilow=subs(I1,w,w1);Iv=factor(simplify(Iup-Ilow)Iv1=subs(subs(char(Iv),x,X0(ID(i),JD(i),y,Y0(ID(i),JD(i)4.7 直射光区域内任意一点出的光照强度确定综上所述，直射光区域内任意一点的光照强度为：因此，直射光区域内的每一点的光照强度是关于坐标点（）的二元非负函数，我们可以通过二元函数的空间曲面来展示光照强度。（也可以用梯度线，梯度向量来展示）。直射光区域内任意一点的光照强度求解程序如下：%直射光的光照度计算clc,clear,close allformat short%功率P0=1;%线光源长度L=4;%整个亮区的矩形范围x0=-38*25000/6.5:1000:38*25000/6.5;y0=-(38*25000/6.5+2*(25016-22.5)/6.5):1000:(38*25000/6.5+2*(25016-22.5)/6.5);%区域D1中的点（x0,y0,25016）X0,Y0=meshgrid(x0,y0);E0=zeros(size(X0);ID,JD=find(X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2 & (X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2);For i=1:length(ID)E0(ID(i),JD(i)=25000*P0/L/(X0(ID(i),JD(i)2+250002)*(2-Y0(ID(i),JD(i)/(X0(ID(i),JD(i)2+(Y0(ID(i),JD(i)-2)2+250002)(1/2)+(2+Y0(ID(i),JD(i)/(X0(ID(i),JD(i)2+(Y0(ID(i),JD(i)+2)2+250002)(1/2);end%区域D21中的点（x0,y0,25016）ID,JD=find(X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2);For i=1:length(ID)E0(ID(i),JD(i)=(-590451689604761308672*(X0(ID(i),JD(i)2+Y0(ID(i),JD(i)2+4*Y0(ID(i),JD(i)+625000004)(1/2)*Y0(ID(i),JD(i)+304451652452455*(X0(ID(i),JD(i)2+Y0(ID(i),JD(i)2+4*Y0(ID(i),JD(i)+625000004)(1/2)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+2*(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2-359527985187027253421989208118379520*Y0(ID(i),JD(i)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+217781166658024909840441558307262072833108122128325)(1/2)+(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2-359527985187027253421989208118379520*Y0(ID(i),JD(i)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+217781166658024909840441558307262072833108122128325)(1/2)*Y0(ID(i),JD(i)/(X0(ID(i),JD(i)2+625000000)/(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2-359527985187027253421989208118379520*Y0(ID(i),JD(i)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+217781166658024909840441558307262072833108122128325)(1/2)/(X0(ID(i),JD(i)2+Y0(ID(i),JD(i)2+4*Y0(ID(i),JD(i)+625000004)(1/2);end%区域D22中的点（x0,y0,25016）ID,JD=find(X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2) & (X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2) & (X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2) &(X0=sqrt(38/6.5*25000)2-(2*(25016-22.5)/6.5)2) & (X0=(38*25000/6.5)2) & (X0.2+(Y0-2*(25016-22.5)/6.5).2=(38*25000/6.5)2) &(X0=-6/1.1*25000);for i=1:length(ID)E0(ID(i),JD(i)=-(590451689604761308672*(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2-359527985187027253421989208118379520*Y0(ID(i),JD(i)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+217781166658024909840441558307262072833108122128325)(1/2)*Y0(ID(i),JD(i)+304451652452455*(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2-359527985187027253421989208118379520*Y0(ID(i),JD(i)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+217781166658024909840441558307262072833108122128325)(1/2)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)-590451689604761308672*(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2+359527985187027253421989208118379520*Y0(ID(i),JD(i)*(4874578512396693-262144*X0(ID(i),JD(i)2)(1/2)+217781166658024909840441558307262072833108122128325)(1/2)*Y0(ID(i),JD(i)+304451652452455*(348449119428701635706449676008113452089344*X0(ID(i),JD(i)2+348633197757117393660201734482670062403584*Y0(ID(i),JD(i)2+3595279851870272534219892081